add effective saturation wrt ice/liquid from Jouzel&Merlivat 1984

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))

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)"

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": {

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)

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]]
+        )