diff --git a/PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py b/PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py
index f17186751..f0c6814c5 100644
--- a/PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py
+++ b/PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py
@@ -11,14 +11,14 @@ def __init__(self, _):
pass
@staticmethod
- def alpha_kinetic(alpha_equilibrium, saturation, D_ratio_heavy_to_light):
+ def alpha_kinetic(alpha_equilibrium, relative_humidity, D_ratio_heavy_to_light):
"""eq. (11)
Parameters
----------
alpha_equilibrium
Equilibrium fractionation factor.
- saturation
+ relative_humidity
Over liquid water or ice.
D_ratio_heavy_to_light
Diffusivity ratio for heavy to light isotope.
@@ -27,6 +27,6 @@ def alpha_kinetic(alpha_equilibrium, saturation, D_ratio_heavy_to_light):
----------
alpha_kinetic
Kinetic fractionation factor for liquid water or ice."""
- return saturation / (
- alpha_equilibrium / D_ratio_heavy_to_light * (saturation - 1) + 1
+ return relative_humidity / (
+ alpha_equilibrium / D_ratio_heavy_to_light * (relative_humidity - 1) + 1
)
diff --git a/examples/PySDM_examples/Fisher_1991/fig_2.ipynb b/examples/PySDM_examples/Fisher_1991/fig_2.ipynb
index ceaab1041..9a6df827b 100644
--- a/examples/PySDM_examples/Fisher_1991/fig_2.ipynb
+++ b/examples/PySDM_examples/Fisher_1991/fig_2.ipynb
@@ -21,15 +21,13 @@
},
{
"cell_type": "code",
- "execution_count": 1,
"id": "31481014375c7d36",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:17.048483Z",
- "start_time": "2025-06-25T12:41:17.041726Z"
+ "end_time": "2025-12-20T08:39:20.178984Z",
+ "start_time": "2025-12-20T08:39:20.175130Z"
}
},
- "outputs": [],
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
@@ -37,19 +35,19 @@
" !pip --quiet install open-atmos-jupyter-utils\n",
" from open_atmos_jupyter_utils import pip_install_on_colab\n",
" pip_install_on_colab('PySDM-examples', 'PySDM')"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 1
},
{
"cell_type": "code",
- "execution_count": 2,
"id": "c36c2f5c-17c1-422e-9ab9-f7d612064199",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:21.768860Z",
- "start_time": "2025-06-25T12:41:17.056517Z"
+ "end_time": "2025-12-20T08:39:22.493785Z",
+ "start_time": "2025-12-20T08:39:20.182948Z"
}
},
- "outputs": [],
"source": [
"from matplotlib import pyplot\n",
"import numpy as np\n",
@@ -64,19 +62,19 @@
" vapour_mixing_ratio,\n",
" ice_saturation_curve_4\n",
")"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 2
},
{
"cell_type": "code",
- "execution_count": 3,
"id": "91bb290498429f53",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:21.875872Z",
- "start_time": "2025-06-25T12:41:21.845639Z"
+ "end_time": "2025-12-20T08:39:22.528300Z",
+ "start_time": "2025-12-20T08:39:22.497806Z"
}
},
- "outputs": [],
"source": [
"formulae= Formulae(\n",
" isotope_meteoric_water_line=\"Dansgaard1964\",\n",
@@ -94,39 +92,39 @@
"for isotope in isotopes:\n",
" alpha_eq[isotope] = getattr(formulae.isotope_equilibrium_fractionation_factors, f'alpha_i_{isotope}')\n",
" diffusivity_ratio[isotope] = getattr(formulae.isotope_diffusivity_ratios, f'ratio_{isotope}_heavy_to_light')"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 3
},
{
"cell_type": "code",
- "execution_count": 4,
"id": "754af83e4ab216bc",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:21.884154Z",
- "start_time": "2025-06-25T12:41:21.881306Z"
+ "end_time": "2025-12-20T08:39:22.533972Z",
+ "start_time": "2025-12-20T08:39:22.531976Z"
}
},
- "outputs": [],
"source": [
"def alpha_kin(iso, T):\n",
" return formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n",
" alpha_equilibrium = alpha_eq[iso](T),\n",
" D_ratio_heavy_to_light=diffusivity_ratio[iso](T),\n",
- " saturation = ice_saturation_curve_4(const=const, T=T)\n",
+ " relative_humidity = ice_saturation_curve_4(const=const, T=T)\n",
" )"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 4
},
{
"cell_type": "code",
- "execution_count": 5,
"id": "9ca4e3256065274b",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:21.894710Z",
- "start_time": "2025-06-25T12:41:21.891746Z"
+ "end_time": "2025-12-20T08:39:22.540134Z",
+ "start_time": "2025-12-20T08:39:22.537463Z"
}
},
- "outputs": [],
"source": [
"def d_delta_dT(T, delta):\n",
" y = yf(T=T)\n",
@@ -143,19 +141,19 @@
" / (alpha * (y + alpha * y_e))\n",
" )\n",
" return res"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 5
},
{
"cell_type": "code",
- "execution_count": 6,
"id": "71979cfe7348344d",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:21.901806Z",
- "start_time": "2025-06-25T12:41:21.899568Z"
+ "end_time": "2025-12-20T08:39:22.544815Z",
+ "start_time": "2025-12-20T08:39:22.542865Z"
}
},
- "outputs": [],
"source": [
"delta_18O_0 = -15 * PER_MILLE\n",
"delta_2H_0 = const.CRAIG_1961_SLOPE_COEFF * delta_18O_0\n",
@@ -163,19 +161,19 @@
"\n",
"y_e = 0\n",
"yf = partial(vapour_mixing_ratio, formulae)\n"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 6
},
{
"cell_type": "code",
- "execution_count": 7,
"id": "4e112378083e50f5",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:25.203838Z",
- "start_time": "2025-06-25T12:41:21.907105Z"
+ "end_time": "2025-12-20T08:39:27.719208Z",
+ "start_time": "2025-12-20T08:39:22.547553Z"
}
},
- "outputs": [],
"source": [
"result = solve_ivp(\n",
" fun=d_delta_dT,\n",
@@ -189,1297 +187,19 @@
" delta_2H=delta_2H,\n",
" delta_18O=delta_18O\n",
")"
- ]
+ ],
+ "outputs": [],
+ "execution_count": 7
},
{
"cell_type": "code",
- "execution_count": 8,
"id": "153a6e26bc2be38a",
"metadata": {
"ExecuteTime": {
- "end_time": "2025-06-25T12:41:25.731698Z",
- "start_time": "2025-06-25T12:41:25.217926Z"
+ "end_time": "2025-12-20T08:39:28.456103Z",
+ "start_time": "2025-12-20T08:39:27.743940Z"
}
},
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- "\n",
- "\n",
- "\n"
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "application/vnd.jupyter.widget-view+json": {
- "model_id": "d85508aab73e4fbc91cf47f982b6d8a7",
- "version_major": 2,
- "version_minor": 0
- },
- "text/plain": [
- "HBox(children=(HTML(value=\"./fig_1.pdf
\"), HTML(value=\"./fig_1.pdf
\"), HTML(value=\"./fig_2.pdf
\"), HTML(value=\"./fig_2.pdf
\"), HTML(value=\"\n",
- "\n",
- "\n"
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "application/vnd.jupyter.widget-view+json": {
- "model_id": "8c15bde98c9c4698979a4cb47ad75d48",
- "version_major": 2,
- "version_minor": 0
- },
- "text/plain": [
- "HBox(children=(HTML(value=\"./fig_8.pdf
\"), HTML(value=\""
+ ],
+ "image/svg+xml": "\n\n\n"
+ },
+ "metadata": {},
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
+ },
+ {
+ "data": {
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_8.pdf
\"), HTML(value=\"./fig_9.pdf
\"), HTML(value=\""
+ ],
+ "image/svg+xml": "\n\n\n"
+ },
+ "metadata": {},
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
+ },
+ {
+ "data": {
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./fig_9.pdf
\"), HTML(value=\""
+ ],
+ "image/svg+xml": "\n\n\n"
+ },
+ "metadata": {},
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
+ },
+ {
+ "data": {
+ "text/plain": [
+ "HBox(children=(HTML(value=\"./tmp9botejd8.pdf
\"), HTML(value…"
+ ],
+ "application/vnd.jupyter.widget-view+json": {
+ "version_major": 2,
+ "version_minor": 0,
+ "model_id": "8efa2723e9c3471a947d3a9f7ce74f52"
+ }
+ },
+ "metadata": {},
+ "output_type": "display_data",
+ "jetTransient": {
+ "display_id": null
+ }
+ }
+ ],
+ "execution_count": 5
+ },
+ {
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2025-12-20T08:38:41.676493Z",
+ "start_time": "2025-12-20T08:38:41.673218Z"
+ }
+ },
+ "cell_type": "code",
+ "source": "",
+ "id": "b03a73fd279fa0ce",
+ "outputs": [],
+ "execution_count": null
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py b/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py
index 84feeb7ed..c9dfd44a7 100644
--- a/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py
+++ b/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py
@@ -22,20 +22,20 @@ def test_units_alpha_kinetic():
# arrange
alpha_eq = 1 * physics.si.dimensionless
D_ratio = 1 * physics.si.dimensionless
- saturation_over_ice = 1 * physics.si.dimensionless
+ rh_over_ice = 1 * physics.si.dimensionless
# act
sut = JouzelAndMerlivat1984.alpha_kinetic(
alpha_equilibrium=alpha_eq,
D_ratio_heavy_to_light=D_ratio,
- saturation=saturation_over_ice,
+ relative_humidity=rh_over_ice,
)
# assert
assert sut.check("[]")
@staticmethod
- def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
+ def test_fig_9_from_jouzel_and_merlivat_1984(plot=PLOT):
"""[Jouzel & Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749)"""
# arrange
formulae = Formulae(
@@ -44,7 +44,7 @@ def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
isotope_diffusivity_ratios="Stewart1975",
)
temperatures = formulae.trivia.C2K(np.asarray([-30, -20, -10]))
- saturation = np.linspace(start=1, stop=1.35)
+ rh = np.linspace(start=1, stop=1.35)
alpha_s = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O
diffusivity_ratio_heavy_to_light = (
formulae.isotope_diffusivity_ratios.ratio_18O_heavy_to_light
@@ -56,7 +56,7 @@ def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
alpha_k = {
temperature: sut(
alpha_equilibrium=alpha_s[temperature],
- saturation=saturation,
+ relative_humidity=rh,
D_ratio_heavy_to_light=diffusivity_ratio_heavy_to_light(temperature),
)
for temperature in temperatures
@@ -68,12 +68,12 @@ def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
}
# plot
- pyplot.xlim(saturation[0], saturation[-1])
+ pyplot.xlim(rh[0], rh[-1])
pyplot.ylim(1.003, 1.022)
- pyplot.xlabel("S")
- pyplot.ylabel("alpha_k * alpha_s")
+ pyplot.xlabel("Si [1]")
+ pyplot.ylabel(r"$\alpha_\text{kin} \alpha_\text{eq}$ [1]")
for k, v in alpha_s_times_alpha_k.items():
- pyplot.plot(saturation, v, label=k)
+ pyplot.plot(rh, v, label=k)
pyplot.legend()
if plot:
pyplot.show()
@@ -88,10 +88,10 @@ def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
@staticmethod
@pytest.mark.parametrize(
- ("temperature_C", "saturation", "alpha"),
+ ("temperature_C", "relative_humidity", "alpha"),
((-10, 1, 1.021), (-10, 1.35, 1.0075), (-30, 1, 1.0174), (-30, 1.35, 1.004)),
)
- def test_fig9_values(temperature_C, saturation, alpha):
+ def test_fig9_values(temperature_C, relative_humidity, alpha):
# arrange
formulae = Formulae(
isotope_kinetic_fractionation_factors="JouzelAndMerlivat1984",
@@ -105,7 +105,7 @@ def test_fig9_values(temperature_C, saturation, alpha):
alpha_s = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O(T)
alpha_k = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(
alpha_equilibrium=alpha_s,
- saturation=saturation,
+ relative_humidity=relative_humidity,
D_ratio_heavy_to_light=diffusivity_ratio_18O(T),
)
@@ -123,13 +123,14 @@ def test_alpha_kinetic_jouzel_merlivat_vs_craig_gordon(
):
# arrange
T = Formulae().trivia.C2K(temperature_C)
- RH = np.linspace(0.3, 1)
+ rh = np.linspace(0.3, 1)
formulae = Formulae(
isotope_equilibrium_fractionation_factors="VanHook1968",
isotope_diffusivity_ratios="HellmannAndHarvey2020",
isotope_kinetic_fractionation_factors="JouzelAndMerlivat1984",
)
- Si = formulae.saturation_vapour_pressure.pvs_ice(T)
+ vapour_partial_pressure = rh * formulae.saturation_vapour_pressure.pvs_water(T)
+ Si = vapour_partial_pressure / formulae.saturation_vapour_pressure.pvs_ice(T)
alpha_eq = getattr(
formulae.isotope_equilibrium_fractionation_factors, f"alpha_l_{isotope}"
)(T)
@@ -138,14 +139,14 @@ def test_alpha_kinetic_jouzel_merlivat_vs_craig_gordon(
)(T)
alpha_kin_jm = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(
alpha_equilibrium=alpha_eq,
- saturation=Si,
+ relative_humidity=Si,
D_ratio_heavy_to_light=D_heavy_to_light,
)
formulae = Formulae(
isotope_kinetic_fractionation_factors="CraigGordon",
)
alpha_kin_cg = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(
- relative_humidity=RH,
+ relative_humidity=rh,
turbulence_parameter_n=1,
delta_diff=alpha_eq - 1,
theta=1,
@@ -155,7 +156,7 @@ def test_alpha_kinetic_jouzel_merlivat_vs_craig_gordon(
n = (alpha_kin_jm + 1) / (alpha_kin_cg + 1)
# plot
- pyplot.plot(1 - RH, n)
+ pyplot.plot(1 - rh, n)
pyplot.gca().set(
xlabel="1-RH",
ylabel="turbulence parameter n",