new formula for isotopic kinetic fractionation factor: Bolot 2013 (#1611)

Author: AgnieszkaZaba

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,22 @@ 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
         )

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

@@ -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,8 @@ 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
@@ -54,10 +56,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 +105,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)