Added blowing snow probability estimation.

Calculations/parameterizations.py

@@ -8,9 +8,27 @@ Collection of parameterizations for INP concentrations and INAS densities.
 """
 
 import numpy as np
+import numpy.typing as npt
 
 
-def inas_parameterizations(T, material):
+def inas_parameterizations(T: npt.ArrayLike, material: str) -> npt.ArrayLike:
+    """
+    Ice Nucleation Active Site (INAS) density parameterizations for different
+    materials.
+
+    Parameters
+    ----------
+    T : npt.ArrayLike
+        Temperature in K.
+    material : str
+        Aerosol type.
+
+    Returns
+    -------
+    npt.ArrayLike
+        INAS density in cm-2.
+
+    """
     match material:
         case 'dust':
             a = 1e-4
@@ -63,3 +81,40 @@ def inas_parameterizations(T, material):
             b0 = 141
             b1 = -0.495
             return 1e-4 * np.exp(b1 * (T + 273.15)**1 + b0)
+
+
+def blowing_snow_probability(
+    wind_speed: npt.ArrayLike, temperature: npt.ArrayLike, snow_age: npt.ArrayLike
+) -> npt.ArrayLike:
+    """
+    Blowing snow probability as a function of the wind speed, temperature and snow age.
+    Parameterization developed by _[1] and valid for wind speeds between 3 and 20 m/s,
+    and for air temperatures above -35 °C.
+
+    Parameters
+    ----------
+    wind_speed : npt.ArrayLike
+        Wind speed in m/s.
+    temperature : npt.ArrayLike
+        Temperature in °C.
+    snow_age : npt.ArrayLike
+        Snow age in hours.
+
+    Returns
+    -------
+    npt.ArrayLike
+        Probability of blowing snow.
+
+    References
+    ----------
+    [1] Li, L. and J. W. Pomeroy (1997). “Probability of occurrence of blowing snow”.
+        In: Journal of Geophysical Research: Atmospheres 102.D18, pp. 21955–21964.
+        issn: 0148-0227. doi: 10.1029/97jd01522.
+    """
+    from scipy.special import erf
+
+    T = temperature
+    I = np.log(snow_age)
+    mean_wind_speed = 0.365 * T + 0.00706 * T**2 + 0.9 * I + 11.2
+    delta = 0.145 * T + 0.00196 * T**2 + 4.3
+    return -1 / 2 * erf((mean_wind_speed - wind_speed) / (np.sqrt(2) * delta)) + 0.5