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 f17186751d..783388dc6a 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 @@ -30,3 +30,39 @@ def alpha_kinetic(alpha_equilibrium, saturation, D_ratio_heavy_to_light): return saturation / ( alpha_equilibrium / D_ratio_heavy_to_light * (saturation - 1) + 1 ) + + @staticmethod + def transfer_coefficient(const, rho_s, D, Fk): + """ + eq. (A4) in Jouzel & Merlivat 1975, + eq. (A6) in Bolot 2013. + + Parameters + ---------- + D + Diffusion coefficient for light isotope. + Fk + Term associated with heat transfer. + + Returns + ---------- + Thermal transfer coefficient between water vapour and condensate. + """ + return 1 / (1 + D * rho_s / const.rho_w * Fk) + + @staticmethod + def effective_saturation(transfer_coefficient, RH): + """ + + Parameters + ---------- + transfer_coefficient + Thermal transfer coefficient between water vapour and condensate (liquid water or ice). + RH + Relative humidity. + + Returns + ---------- + Effective vapour saturation over liquid water or ice. + """ + return 1 / (1 - transfer_coefficient * (1 - 1 / RH)) diff --git a/docs/bibliography.json b/docs/bibliography.json index 337dc2a2b6..fafca06568 100644 --- a/docs/bibliography.json +++ b/docs/bibliography.json @@ -881,7 +881,8 @@ "examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb", "examples/PySDM_examples/Jouzel_and_Merlivat_1984/thermodynamic_profiles.py", "tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py", - "examples/PySDM_examples/Fisher_1991/__init__.py" + "examples/PySDM_examples/Fisher_1991/__init__.py", + "tests/smoke_tests/no_env/jouzel_and_merlivat_1984/test_fig_8.py" ], "title": "Deuterium and oxygen 18 in precipitation: Modeling of the isotopic effects during snow formation", "label": "Jouzel & Merlivat 1984 (J. Geophys. Res. Atmos. 89)" diff --git a/examples/PySDM_examples/Fisher_1991/fig_2.ipynb b/examples/PySDM_examples/Fisher_1991/fig_2.ipynb index 6ebee29c8b..9e786c967f 100644 --- a/examples/PySDM_examples/Fisher_1991/fig_2.ipynb +++ b/examples/PySDM_examples/Fisher_1991/fig_2.ipynb @@ -36,7 +36,7 @@ ], "id": "31481014375c7d36", "outputs": [], - "execution_count": 1 + "execution_count": 75 }, { "cell_type": "code", @@ -63,7 +63,7 @@ ")" ], "outputs": [], - "execution_count": 2 + "execution_count": 76 }, { "metadata": { @@ -91,9 +91,9 @@ " 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')" ], - "id": "91bb290498429f53", + "id": "92e11666a7b83d4d", "outputs": [], - "execution_count": 3 + "execution_count": 77 }, { "metadata": { @@ -111,9 +111,9 @@ " saturation = ice_saturation_curve_4(const=const, T=T)\n", " )" ], - "id": "754af83e4ab216bc", + "id": "cd4693f604f74961", "outputs": [], - "execution_count": 4 + "execution_count": 78 }, { "metadata": { @@ -142,7 +142,7 @@ ], "id": "9ca4e3256065274b", "outputs": [], - "execution_count": 5 + "execution_count": 79 }, { "metadata": { @@ -160,9 +160,9 @@ "y_e = 0\n", "yf = partial(vapour_mixing_ratio, formulae)\n" ], - "id": "71979cfe7348344d", + "id": "7bce3bc5674ae8cb", "outputs": [], - "execution_count": 6 + "execution_count": 80 }, { "metadata": { @@ -180,15 +180,11 @@ " t_eval=temperature\n", ")\n", "assert result.success, result.message\n", - "delta_2H, delta_18O = result.y\n", - "d_excess = formulae.isotope_meteoric_water_line.excess_d(\n", - " delta_2H=delta_2H,\n", - " delta_18O=delta_18O\n", - ")" + "delta_2H, delta_18O = result.y\n" ], - "id": "4e112378083e50f5", + "id": "200252902b10ccb6", "outputs": [], - "execution_count": 7 + "execution_count": 81 }, { "metadata": { @@ -225,7 +221,7 @@ "pyplot.legend()\n", "show_plot(\"fig_1\")" ], - "id": "153a6e26bc2be38a", + "id": "262ff6a72a6d27ec", "outputs": [ { "data": { @@ -252,7 +248,7 @@ "output_type": "display_data" } ], - "execution_count": 8 + "execution_count": 82 }, { "metadata": { @@ -265,10 +261,13 @@ "source": [ "d_excess_from_fig2 = np.array((10, 8, 7, 6.8, 6, 5.3, 4, 2.5, 0)) * PER_MILLE\n", "T_from_fig2 = np.linspace(-50, -10, n_points)\n", - "\n", + "d_excess = formulae.isotope_meteoric_water_line.excess_d(\n", + " delta_2H=delta_2H,\n", + " delta_18O=delta_18O\n", + ")\n", "\n", "pyplot.plot(\n", - " K2C(temperature[::-1]),\n", + " K2C(temperature)[::-1],\n", " in_unit(d_excess, PER_MILLE)[::-1],\n", " label=\"d-excess, TODO #1601 (WORK IN PROGRESS - NO MATCH WITH THE PAPER YET)\",\n", ")\n", @@ -281,7 +280,7 @@ "pyplot.legend()\n", "show_plot(\"fig_2\")" ], - "id": "571d258dc7dea96f", + "id": "2cf7385894279338", "outputs": [ { "data": { diff --git a/examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb b/examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb index 446d44e577..e38ee8b6fb 100644 --- a/examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb +++ b/examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb @@ -14,16 +14,17 @@ "metadata": {}, "cell_type": "markdown", "source": [ - "# figs 8 and 9 from [Jouzel and Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749) \"_Deuterium and oxygen 18 in precipitation: Modeling of the isotopic effects during snow formation_\"\n", - "\n" + "# Comparison between saturation wrt. ice and effective saturation - based on [Jouzel and Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749) \"_Deuterium and oxygen 18 in precipitation: Modeling of the isotopic effects during snow formation_\"\n", + "\n", + "figs 8 and 9" ], "id": "841538e30be219a1" }, { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:40.862379Z", - "start_time": "2025-06-25T23:25:40.859115Z" + "end_time": "2025-07-10T10:48:11.678433Z", + "start_time": "2025-07-10T10:48:11.674981Z" } }, "cell_type": "code", @@ -44,8 +45,8 @@ "metadata": { "collapsed": true, "ExecuteTime": { - "end_time": "2025-06-25T23:25:42.736860Z", - "start_time": "2025-06-25T23:25:40.951437Z" + "end_time": "2025-07-10T10:48:13.278425Z", + "start_time": "2025-07-10T10:48:11.743948Z" } }, "source": [ @@ -62,16 +63,21 @@ { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:43.030796Z", - "start_time": "2025-06-25T23:25:42.742196Z" + "end_time": "2025-07-10T10:48:13.549139Z", + "start_time": "2025-07-10T10:48:13.285324Z" } }, "cell_type": "code", "source": [ - "formulae = Formulae(\n", + "ANY_NUMBER = 44\n", + "\n", + "formulae= Formulae(\n", + " latent_heat_vapourisation=\"Lowe2019\",\n", + " saturation_vapour_pressure=\"MurphyKoop2005\",\n", " isotope_diffusivity_ratios=\"Stewart1975\",\n", - " isotope_equilibrium_fractionation_factors=\"Majoube1970\",\n", - " isotope_kinetic_fractionation_factors=\"JouzelAndMerlivat1984\"\n", + " isotope_kinetic_fractionation_factors=\"JouzelAndMerlivat1984\",\n", + " isotope_equilibrium_fractionation_factors=\"MerlivatAndNief1967+Majoube1970\",\n", + " diffusion_thermics=\"Neglect\",\n", ")\n", "const = formulae.constants\n", "svp = formulae.saturation_vapour_pressure\n", @@ -81,46 +87,151 @@ "n_points = 100\n", "T_0_50 = C2K(np.linspace(0.0, -50, n_points))" ], - "id": "87e416d611926ad5", + "id": "6f7078d9d56a0a7c", "outputs": [], "execution_count": 3 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:44.982409Z", - "start_time": "2025-06-25T23:25:43.041753Z" + "end_time": "2025-07-10T10:48:13.561668Z", + "start_time": "2025-07-10T10:48:13.559532Z" + } + }, + "cell_type": "code", + "source": [ + "diffusivity_ratio = {}\n", + "alpha_eq = {}\n", + "alpha_kinetic = {}\n", + "isotopes = (\"2H\", \"18O\")\n", + "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')" + ], + "id": "2c7e390144e7499", + "outputs": [], + "execution_count": 4 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "## Saturation wrt ice - comparison with effective saturation", + "id": "bbb54a367a947a49" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-10T10:48:15.236361Z", + "start_time": "2025-07-10T10:48:13.569626Z" } }, "cell_type": "code", "source": [ + "pvs_ice = formulae.saturation_vapour_pressure.pvs_ice\n", + "pvs_water = formulae.saturation_vapour_pressure.pvs_water\n", + "pressure = ANY_NUMBER\n", + "\n", + "D_light = formulae.diffusion_thermics.D(T=T_0_50, p=pressure)\n", + "Fk = formulae.drop_growth.Fk(\n", + " T=T_0_50,\n", + " K=formulae.diffusion_thermics.K(T=T_0_50, p=pressure),\n", + " lv=formulae.latent_heat_vapourisation.lv(T_0_50)\n", + ")\n", + "\n", + "saturation_wrt_ice_at_RH100 = pvs_water(T_0_50) / pvs_ice(T_0_50)\n", + "\n", + "eff_saturation_wrt_ice_at_RH100 = formulae.isotope_kinetic_fractionation_factors.effective_saturation(\n", + " transfer_coefficient=formulae.isotope_kinetic_fractionation_factors.transfer_coefficient(\n", + " rho_s = pvs_ice(T_0_50) / const.Rv / T_0_50,\n", + " D=D_light,\n", + " Fk=Fk\n", + "),\n", + " RH=saturation_wrt_ice_at_RH100\n", + ")" + ], + "id": "47503ee93cc639b7", + "outputs": [], + "execution_count": 5 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-10T10:48:16.174530Z", + "start_time": "2025-07-10T10:48:15.248702Z" + } + }, + "cell_type": "code", + "source": [ + "isotope_considered = '2H'\n", + "alpha_kinetic = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n", + " alpha_equilibrium = alpha_eq[isotope_considered](T_0_50),\n", + " D_ratio_heavy_to_light=diffusivity_ratio[isotope_considered](T_0_50),\n", + " saturation = saturation_wrt_ice_at_RH100\n", + " )\n", + "\n", + "eff_alpha_kinetic = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n", + " alpha_equilibrium = alpha_eq[isotope_considered](T_0_50),\n", + " D_ratio_heavy_to_light=diffusivity_ratio[isotope_considered](T_0_50),\n", + " saturation = eff_saturation_wrt_ice_at_RH100\n", + " )" + ], + "id": "ad049bb8d20fa39b", + "outputs": [], + "execution_count": 6 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-10T10:48:16.674907Z", + "start_time": "2025-07-10T10:48:16.184371Z" + } + }, + "cell_type": "code", + "source": [ + "pyplot.plot(\n", + " K2C(T_0_50),\n", + " saturation_wrt_ice_at_RH100,\n", + " 'k',\n", + " label='saturation wrt ice at RH100',\n", + ")\n", "pyplot.plot(\n", " K2C(T_0_50),\n", - " svp.pvs_water(T_0_50) / svp.pvs_ice(T_0_50),\n", - " label='saturation wrt liquid phase',\n", + " pvs_ice(T_0_50),\n", + " 'b',\n", + " label='',\n", + ")\n", + "pyplot.plot(\n", + " K2C(T_0_50),\n", + " eff_saturation_wrt_ice_at_RH100,\n", + " 'k--',\n", + " label=\"effective saturation wrt ice at RH100\",\n", ")\n", "pyplot.plot(\n", " K2C(T_0_50)[n_points//5:],\n", " ice_saturation_curve_4(const,T_0_50[n_points//5:]),\n", + " 'g',\n", " label='curve 4',\n", ")\n", + "pyplot.legend()\n", "pyplot.gca().set(\n", - " ylabel='Si',\n", - " xlabel='Temperature (C)',\n", + " title=f\"Saturation wrt ice for {isotope_considered} at RH=100\",\n", + " xlabel=\"Temperature [°C]\",\n", + " ylabel=\"Si [ ]\",\n", " xlim=(K2C(T_0_50[0]), K2C(T_0_50[-1])),\n", + " ylim=(1, 1.7)\n", ")\n", "pyplot.grid()\n", "pyplot.legend()\n", "show_plot('fig_8')" ], - "id": "3f97b336a478a65a", + "id": "fd51ae9370fec8d7", "outputs": [ { "data": { "text/plain": [ "
" ], - "image/svg+xml": "\n\n\n \n \n \n \n 2025-06-26T01:25:44.970122\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-10T12:48:16.658508\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" }, "metadata": {}, "output_type": "display_data" @@ -133,45 +244,160 @@ "application/vnd.jupyter.widget-view+json": { "version_major": 2, "version_minor": 0, - "model_id": "8c15bde98c9c4698979a4cb47ad75d48" + "model_id": "a1dfa4c3b72242c39e736397d22a942d" } }, "metadata": {}, "output_type": "display_data" } ], - "execution_count": 4 + "execution_count": 7 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:45.731460Z", - "start_time": "2025-06-25T23:25:44.995396Z" + "end_time": "2025-07-10T10:48:16.991147Z", + "start_time": "2025-07-10T10:48:16.706807Z" + } + }, + "cell_type": "code", + "source": [ + "saturation_difference = saturation_wrt_ice_at_RH100 - eff_saturation_wrt_ice_at_RH100\n", + "pyplot.plot(K2C(T_0_50), abs(saturation_difference), 'k')\n", + "pyplot.grid()\n", + "pyplot.gca().set(\n", + " title=\"Absolute difference between environment and effective saturation wrt ice\",\n", + " xlabel=\"temperature [C]\",\n", + " ylabel=\"$|Si - S_{{eff}}|$ [ ]\",\n", + ")\n", + "show_plot(\"Si_minus_effective_Si\")" + ], + "id": "2573a83fb26b7f0f", + "outputs": [ + { + "data": { + "text/plain": [ + "
" + ], + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-10T12:48:16.975254\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\"./Si_minus_effective_Si.pdf" + ], + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-10T12:48:17.136976\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\"./alpha_kinetic_of_tem…" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "7e88c170759c47f68245fadf3f3ad6f5" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 9 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "## Fig. 9: Effective fractionation factors\n", + "$\\alpha_{eff} = \\alpha_s\\alpha_k$\n", + "for $^{18}$O of saturation over ice" + ], + "id": "552c5304753132f5" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-10T10:48:18.029386Z", + "start_time": "2025-07-10T10:48:17.159609Z" } }, "cell_type": "code", "source": [ "temperatures = C2K(np.array([-10., -20., -30.]))\n", + "isotope_considered = \"18O\"\n", + "\n", "alphas_eff = np.zeros((n_points, len(temperatures)))\n", "Si = np.linspace(1.0, 1.4, n_points)\n", "for i, T in enumerate(temperatures):\n", - " alpha_s = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O(T)\n", + " alpha_s = alpha_eq[isotope_considered](T)\n", " alpha_k = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n", " alpha_equilibrium=alpha_s,\n", " saturation=Si,\n", - " D_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_18O_heavy_to_light(T)\n", + " D_ratio_heavy_to_light=diffusivity_ratio[isotope_considered](T),\n", " )\n", " alphas_eff[:, i] = alpha_k * alpha_s" ], "id": "1a7e5ae0c7093428", "outputs": [], - "execution_count": 5 + "execution_count": 10 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:45.959914Z", - "start_time": "2025-06-25T23:25:45.739364Z" + "end_time": "2025-07-10T10:48:18.247014Z", + "start_time": "2025-07-10T10:48:18.054434Z" } }, "cell_type": "code", @@ -180,27 +406,27 @@ " pyplot.plot(\n", " Si,\n", " alphas_eff[:, i],\n", - " label=f'{K2C(T):.5g} (C)',\n", - " color='k',\n", + " 'g',\n", + " label=f'${K2C(T):.5g}$ [°C]',\n", " linewidth = 1 + i/2\n", " )\n", "pyplot.gca().set(\n", " ylabel='$\\\\alpha_s\\\\alpha_k$',\n", " xlabel='Si',\n", - " title = \"$\\\\alpha_s\\\\alpha_k$ for $^{18}$O of saturation over ice\"\n", + " title = f\"Effective fractionation factors for {isotope_considered}\"\n", ")\n", "pyplot.grid()\n", "pyplot.legend()\n", "show_plot('fig_9')" ], - "id": "8b1287b4696c7dda", + "id": "7c07f76691edf597", "outputs": [ { "data": { "text/plain": [ "
" ], - "image/svg+xml": "\n\n\n \n \n \n \n 2025-06-26T01:25:45.940460\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-10T12:48:18.226016\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" }, "metadata": {}, "output_type": "display_data" @@ -213,25 +439,25 @@ "application/vnd.jupyter.widget-view+json": { "version_major": 2, "version_minor": 0, - "model_id": "a478d96b93d34ed2936063de54732783" + "model_id": "df8d0f3528124da78a006fd0b6c88e9c" } }, "metadata": {}, "output_type": "display_data" } ], - "execution_count": 6 + "execution_count": 11 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-06-25T23:25:45.972541Z", - "start_time": "2025-06-25T23:25:45.970842Z" + "end_time": "2025-07-10T10:48:18.272016Z", + "start_time": "2025-07-10T10:48:18.269850Z" } }, "cell_type": "code", "source": "", - "id": "94ff2c9e78541ae", + "id": "d482baa7f1c16462", "outputs": [], "execution_count": null } diff --git a/tests/smoke_tests/no_env/jouzel_and_merlivat_1984/test_fig_8.py b/tests/smoke_tests/no_env/jouzel_and_merlivat_1984/test_fig_8.py new file mode 100644 index 0000000000..1064a89c5f --- /dev/null +++ b/tests/smoke_tests/no_env/jouzel_and_merlivat_1984/test_fig_8.py @@ -0,0 +1,96 @@ +""" +tests for fig. 8 from [Jouzel & Merlivat 1984 (J. Geophys. Res. Atmos. 89)](https://doi.org/10.1029/JD089iD07p11749) +""" # pylint: disable=line-too-long + +from pathlib import Path +import pytest +import numpy as np +from open_atmos_jupyter_utils import notebook_vars + +from PySDM_examples import Jouzel_and_Merlivat_1984 +from PySDM.physics.constants import T0 + +PLOT = False + + +@pytest.fixture(scope="session", name="notebook_variables") +def notebook_variables_fixture(): + return notebook_vars( + file=Path(Jouzel_and_Merlivat_1984.__file__).parent / "fig_8_9.ipynb", + plot=PLOT, + ) + + +class TestFig8: + @staticmethod + @pytest.mark.parametrize( + "if_effective, temp_C, value", + ( + (False, 0, 1), + (False, -10, 1.1), + (False, -20, 1.22), + (False, -30, 1.34), + (False, -45, 1.54), + (True, 0, 1), + (True, -10, 1.08), + (True, -20, 1.18), + (True, -30, 1.30), + (True, -45, 1.52), + ), + ) + def test_fig8_values_against_the_paper( + notebook_variables, if_effective, temp_C, value + ): + # arrange + temperature_C = notebook_variables["T_0_50"] - T0 + if if_effective: + Si = notebook_variables["eff_saturation_wrt_ice_at_RH100"] + else: + Si = notebook_variables["saturation_wrt_ice_at_RH100"] + + # act + sut = Si[np.argmin(np.abs(temperature_C - temp_C))] + + # assert + np.testing.assert_approx_equal(sut, value, significant=1) + + @staticmethod + @pytest.mark.parametrize("temperature_C, expected", ((-20, 0.05),)) + def test_fig_8_max_difference(notebook_variables, temperature_C, expected): + """ + Test maximum difference between saturation and effective saturation over ice. + Based on comment the on eq. (13) + in [Jouzel & Merlivat 1984 (J. Geophys. Res. Atmos. 89)](https://doi.org/10.1029/JD089iD07p11749) + """ # pylint: disable=line-too-long + # arrange + saturation_difference = notebook_variables["saturation_difference"] + temp = notebook_variables["T_0_50"] + temp_C = temp - T0 + temp_C_tolerance = 0.01 + + # act + max_difference = np.max(saturation_difference) + temp_C_range_max = temp_C[ + np.isclose(saturation_difference, max_difference, rtol=0.1) + ] + + # assert + assert max_difference <= expected + assert ( + np.min(temp_C_range_max) - temp_C_tolerance + <= temperature_C + <= np.max(temp_C_range_max) + temp_C_tolerance + ) + + @staticmethod + def test_alpha_less_then_eff_alpha(notebook_variables): + # arrange + alpha_eff = notebook_variables["eff_alpha_kinetic"] + alpha = notebook_variables["alpha_kinetic"] + + # act + sut = alpha_eff + + # assert + np.testing.assert_array_less(sut, 1) + np.testing.assert_array_less(alpha, alpha_eff) 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 84feeb7ed4..1a1855c362 100644 --- a/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py +++ b/tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py @@ -34,6 +34,23 @@ def test_units_alpha_kinetic(): # assert assert sut.check("[]") + @staticmethod + def test_units_transfer_coefficient(): + with DimensionalAnalysis(): + # arrange + const = Formulae().constants + rho_s = 1 * physics.si.g / physics.si.m**3 + D = 1 * physics.si.m**2 / physics.si.s + Fk = 1 * physics.si.s / physics.si.m**2 + + # act + sut = JouzelAndMerlivat1984.transfer_coefficient( + const=const, D=D, Fk=Fk, rho_s=rho_s + ) + + # assert + assert sut.check("[]") + @staticmethod def test_fig_9_from_jouzel_and_merlivat_1984(plot=False): """[Jouzel & Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749)""" @@ -170,3 +187,35 @@ def test_alpha_kinetic_jouzel_merlivat_vs_craig_gordon( # assert np.testing.assert_equal(n > 0.5, True) np.testing.assert_equal(n < 1, True) + + @staticmethod + @pytest.mark.parametrize("T", np.linspace(240, 300, 11)) + def test_effective_saturation(T): + # arrange + formulae = Formulae(drop_growth="Mason1971") + const = formulae.constants + saturation = np.linspace(0.8, 1.2, 21) + rho_s = formulae.saturation_vapour_pressure.pvs_ice(T) / const.Rv / T + Fk = formulae.drop_growth.Fk( + T=T, + K=const.K0, + lv=formulae.latent_heat_vapourisation.lv(T), + ) + A_li = JouzelAndMerlivat1984.transfer_coefficient( + const=const, rho_s=rho_s, D=const.D0, Fk=Fk + ) + + # act + sut = JouzelAndMerlivat1984.effective_saturation( + transfer_coefficient=A_li, + RH=saturation, + ) + idx_subsaturation = np.where(saturation < 1) + + # assert + np.testing.assert_array_less( + saturation[idx_subsaturation], sut[idx_subsaturation] + ) + np.testing.assert_array_less( + sut[~idx_subsaturation[0]], saturation[~idx_subsaturation[0]] + )