new formula for isotopic kinetic fractionation factor: Bolot 2013

PySDM/physics/drop_growth/howell_1949.py

@@ -14,7 +14,17 @@ class Howell1949(Fick):  # pylint: disable=too-few-public-methods
 
     @staticmethod
     def Fk(const, T, K, lv):
-        """thermodynamic term associated with heat conduction"""
+        """Thermodynamic term associated with heat conduction.
+
+        Parameters
+        ----------
+        T
+            Temperature.
+        K
+            Thermal diffusivity with heat ventilation factor.
+        lv
+            Latent heat of evaporation or sublimation.
+        """
         return const.rho_w * lv / T / K * (lv / T / const.Rv)
 
     @staticmethod

PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py

@@ -1,6 +1,8 @@
 """
-kinetic fractionation factor from [Jouzel & Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749)
-(as eq. 3e for n=1 in [Stewart 1975](https://doi.org/10.1029/JC080i009p01133))
+kinetic fractionation factor [Jouzel & Merlivat 1984](https://doi.org/10.1029/JD089iD07p11749),
+as eq. 3e for n=1 in [Stewart 1975](https://doi.org/10.1029/JC080i009p01133)
+and eq. 1 in [Bolot 2013](https://doi.org/10.5194/acp-13-7903-2013),
+where alpha_kinetic is multiplied by alpha equilibrium (eq. 1 defines effective alpha).
 """
 
 
@@ -9,13 +11,58 @@ def __init__(self, _):
         pass
 
     @staticmethod
-    def alpha_kinetic(
-        alpha_equilibrium, saturation_over_ice, diffusivity_ratio_heavy_to_light
-    ):
-        """eq. (11)"""
-        return saturation_over_ice / (
-            alpha_equilibrium
-            / diffusivity_ratio_heavy_to_light
-            * (saturation_over_ice - 1)
-            + 1
+    def alpha_kinetic(alpha_equilibrium, saturation, D_ratio_heavy_to_light):
+        """eq. (11)
+
+        Parameters
+        ----------
+        alpha_equilibrium
+            Equilibrium fractionation factor.
+        saturation
+            Over liquid water or ice.
+        D_ratio_heavy_to_light
+            Diffusivity ratio for heavy to light isotope.
+
+        Returns
+        ----------
+        alpha_kinetic
+            Kinetic fractionation factor for liquid water or ice."""
+        return saturation / (
+            alpha_equilibrium / D_ratio_heavy_to_light * (saturation - 1) + 1
         )
+
+    @staticmethod
+    def transfer_coefficient(D, Fk):
+        """
+        eq. (A4) in Jouzel & Merlivat 1984,
+        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 * 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

@@ -529,6 +529,7 @@
   "https://doi.org/10.5194/acp-13-7903-2013": {
     "usages": [
       "examples/PySDM_examples/Bolot_et_al_2013/__init__.py",
+      "PySDM/physics/isotope_kinetic_fractionation_factors/jouzel_and_merlivat_1984.py",
       "tests/unit_tests/physics/test_isotope_equilibrium_fractionation_factors.py"
     ],
     "label": "Bolot et al. 2013 (Atmos. Chem. Phys. 13)",

examples/PySDM_examples/Fisher_1991/fig_2.ipynb

@@ -108,7 +108,7 @@
     "    return formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n",
     "        alpha_equilibrium = alpha_eq[iso](T),\n",
     "        diffusivity_ratio_heavy_to_light=diffusivity_ratio[iso](T),\n",
-    "        saturation_over_ice = ice_saturation_curve_4(const=const, T=T)\n",
+    "        saturation = ice_saturation_curve_4(const=const, T=T)\n",
     "    )"
    ],
    "id": "754af83e4ab216bc",

examples/PySDM_examples/Jouzel_and_Merlivat_1984/fig_8_9.ipynb

@@ -158,7 +158,7 @@
     "    alpha_s = formulae.isotope_equilibrium_fractionation_factors.alpha_i_18O(T)\n",
     "    alpha_k = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(\n",
     "        alpha_equilibrium=alpha_s,\n",
-    "        saturation_over_ice=Si,\n",
+    "        saturation=Si,\n",
     "        diffusivity_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_18O_heavy_to_light(T)\n",
     "    )\n",
     "    alphas_eff[:, i] = alpha_k * alpha_s"

tests/unit_tests/physics/test_isotope_kinetic_fractionation_factors.py

@@ -11,10 +11,12 @@
 from PySDM.physics.dimensional_analysis import DimensionalAnalysis
 from PySDM.physics.isotope_kinetic_fractionation_factors import JouzelAndMerlivat1984
 
+PLOT = False
+
 
 class TestIsotopeKineticFractionationFactors:
     @staticmethod
-    def test_units():
+    def test_units_alpha_kinetic():
         """checks that alphas are dimensionless"""
         with DimensionalAnalysis():
             # arrange
@@ -25,8 +27,36 @@ def test_units():
             # act
             sut = JouzelAndMerlivat1984.alpha_kinetic(
                 alpha_equilibrium=alpha_eq,
-                diffusivity_ratio_heavy_to_light=D_ratio,
-                saturation_over_ice=saturation_over_ice,
+                D_ratio_heavy_to_light=D_ratio,
+                saturation=saturation_over_ice,
+            )
+
+            # assert
+            assert sut.check("[]")
+
+    @staticmethod
+    def test_units_transfer_coefficient():
+        with DimensionalAnalysis():
+            # arrange
+            D = 1 * physics.si.m**2 / physics.si.s
+            Fk = 1 * physics.si.s / physics.si.m**2
+
+            # act
+            sut = JouzelAndMerlivat1984.transfer_coefficient(D=D, Fk=Fk)
+
+            # assert
+            assert sut.check("[]")
+
+    @staticmethod
+    def test_units_eff_saturation():
+        with DimensionalAnalysis():
+            # arrange
+            transfer_coeff = 1 * physics.si.dimensionless
+            relative_humidity = 1 * physics.si.dimensionless
+
+            # act
+            sut = JouzelAndMerlivat1984.effective_saturation(
+                transfer_coefficient=transfer_coeff, RH=relative_humidity
             )
 
             # assert
@@ -54,10 +84,8 @@ def test_fig_9_from_jouzel_and_merlivat_1984(plot=False):
         alpha_k = {
             temperature: sut(
                 alpha_equilibrium=alpha_s[temperature],
-                saturation_over_ice=saturation,
-                diffusivity_ratio_heavy_to_light=diffusivity_ratio_heavy_to_light(
-                    temperature
-                ),
+                saturation=saturation,
+                D_ratio_heavy_to_light=diffusivity_ratio_heavy_to_light(temperature),
             )
             for temperature in temperatures
         }
@@ -105,12 +133,68 @@ 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_over_ice=saturation,
-            diffusivity_ratio_heavy_to_light=diffusivity_ratio_18O(T),
+            saturation=saturation,
+            D_ratio_heavy_to_light=diffusivity_ratio_18O(T),
         )
 
         # act
         sut = alpha_s * alpha_k
 
         # assert
         np.testing.assert_approx_equal(actual=sut, desired=alpha, significant=3)
+
+    @staticmethod
+    @pytest.mark.parametrize("isotope", ("2H", "18O", "17O"))
+    @pytest.mark.parametrize("temperature_C", (-30, -20, -1))
+    def test_alpha_kinetic_jouzel_merlivat_vs_craig_gordon(
+        isotope, temperature_C, plot=PLOT
+    ):
+        # arrange
+        T = Formulae().trivia.C2K(temperature_C)
+        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)
+        alpha_eq = getattr(
+            formulae.isotope_equilibrium_fractionation_factors, f"alpha_l_{isotope}"
+        )(T)
+        D_heavy_to_light = getattr(
+            formulae.isotope_diffusivity_ratios, f"ratio_{isotope}_heavy_to_light"
+        )(T)
+        alpha_kin_jm = formulae.isotope_kinetic_fractionation_factors.alpha_kinetic(
+            alpha_equilibrium=alpha_eq,
+            saturation=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,
+            turbulence_parameter_n=1,
+            delta_diff=alpha_eq - 1,
+            theta=1,
+        )
+
+        # act
+        n = (alpha_kin_jm + 1) / (alpha_kin_cg + 1)
+
+        # plot
+        pyplot.plot(1 - RH, n)
+        pyplot.gca().set(
+            xlabel="1-RH",
+            ylabel="turbulence parameter n",
+        )
+        pyplot.grid()
+
+        if plot:
+            pyplot.show()
+        else:
+            pyplot.clf()
+
+        # assert
+        np.testing.assert_equal(n > 0.5, True)
+        np.testing.assert_equal(n < 1, True)