Merge #144

144: Implementing Chen 2022 Terminal Velocity Formulas into 2M scheme r=trontrytel a=apoorva-thanvantri

## The purpose of this pull request is to add the new terminal velocity formulas for raindrops (and potentially ice particles) from Chen 2022 to the 2 Moment Cloud Microphysics Scheme.  

## To-do
#### implement formulas for ice
## Content
#### wrote documentation page! (might need some edits)
#### initial implementation of formulas - still have questions/concerns on size distribution; need to edit
#### made plots comparing terminal velocity formulas between chen and seifert and Ogura(1 moment formula)
#### made unit tests


Co-authored-by: apoorva-thanvantri <apoorva.thanvantri@gmail.com>
Co-authored-by: Anna Jaruga <ajaruga@caltech.edu>

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.11.1"
+version = "0.12.0"
 
 [deps]
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"

docs/bibliography.bib

@@ -35,6 +35,16 @@ @article{Beheng1994
   doi = {10.1016/0169-8095(94)90020-5}
 }
 
+@article{Chen2022,
+  title={Accurate parameterization of precipitation particles' fall speeds for bulk cloud microphysics schemes},
+  author={Chen, Jen-Ping and Hsieh, Ting-Wei and Lin, Chiu-Yi and Yu, Cheng-Ku},
+  journal={Atmospheric Research},
+  volume={293},
+  number={106171},
+  doi = {https://doi.org/10.1016/j.atmosres.2022.106171},
+  year = {2022}
+}
+
 @article{Desai2019,
   title = {Aerosol-Mediated Glaciation of Mixed-Phase Clouds: Steady-State Laboratory Measurements},
   author = {Desai, Neel and Chandrakar, KK and Kinney, G and Cantrell, W and Shaw, RA},
@@ -493,3 +503,5 @@ @article{SeifertBeheng2006
   year={2006},
   publisher={Springer}
 }
+
+

docs/src/API.md

@@ -115,6 +115,7 @@ CommonTypes.B1994Type
 CommonTypes.TC1980Type
 CommonTypes.LD2004Type
 CommonTypes.SB2006Type
+CommonTypes.Chen2022Type
 CommonTypes.AbstractAerosolType
 CommonTypes.ArizonaTestDustType
 CommonTypes.DesertDustType

docs/src/Microphysics2M.md

@@ -20,15 +20,15 @@ The [SeifertBeheng2006](@cite) parametrization provides process rates for autoco
 The piece-wise polynomial collection Kernel, used for the derivation of the parametrization, is given by:
 ```math
 \begin{align}
-    K(x,y) = 
+    K(x,y) =
     \begin{cases}
     k_{cc}(x^2+y^2), \quad & x\wedge y < x^*\\
     k_{cr}(x+y), \quad & x\oplus y \geq x^*,\\
     k_{rr}(x+y)\exp[-\kappa_{rr} (x^{1/3} + y^{1/3})], \quad & x\wedge y \geq x^*,
     \end{cases}
 \end{align}
 ```
-where ``x`` and ``y`` are drop masses and ``x^*`` is the mass threshold chosen to separate the cloud and rain portions of the mass distribution. For ``K`` in ``m^3 s^{-1}`` the constants are 
+where ``x`` and ``y`` are drop masses and ``x^*`` is the mass threshold chosen to separate the cloud and rain portions of the mass distribution. For ``K`` in ``m^3 s^{-1}`` the constants are
 
 |   symbol      | default value                                       |
 |---------------|-----------------------------------------------------|
@@ -79,7 +79,7 @@ where:
   - ``\rho_0 = 1.225 \, kg \cdot m^{-3}`` is the air density at surface conditions,
   - ``k_{cc}`` is the cloud-cloud collection kernel constant,
   - ``\nu`` is the cloud droplet gamma distribution parameter,
-  - ``x^*`` is the drop mass separating the cloud and rain categories 
+  - ``x^*`` is the drop mass separating the cloud and rain categories
   - ``\overline{x}_c = (q_{liq} \rho) / N_{liq}`` is the cloud droplet mean mass with ``N_{liq}`` denoting the cloud droplet number density. Here, to ensure numerical stability, we limit ``\overline{x}_c`` by the upper bound of ``x^*``.
 
 The function ``\phi_{acnv}(\tau)`` is used to correct the autoconversion rate for the undeveloped cloud droplet spectrum and the early stage rain evolution assumptions. This is a universal function which is obtained by fitting to numerical results of the SCE:
@@ -113,7 +113,7 @@ and the rate of change of liquid water specific humidity and cloud droplets numb
   \left. \frac{\partial N_{liq}}{\partial t} \right|_{acnv} = -2 \left. \frac{\partial N_{rai}}{\partial t} \right|_{acnv}.
 \end{align}
 ```
-!!! note 
+!!! note
     The Seifert and Beheng parametrization is formulated for the rate of change of liquid water content ``L = \rho q``. Here, we assume constant ``\rho`` and divide the rates by ``\rho`` to derive the equations for the rate of change of specific humidities.
 
 ### Accretion
@@ -231,7 +231,7 @@ Raindrops breakup is modeled by assuming that in a precipitation event coalescen
 where ``\Delta \overline{D}_r = \overline{D}_r - \overline{D}_{eq}`` with ``\overline{D}_r`` denoting the mean volume raindrop diameter and ``\overline{D}_{eq}`` being the equilibrium mean diameter. The function ``\Phi_{br}(\Delta \overline{D}_r)`` is given by
 ```math
   \begin{align}
-    \Phi_{br}(\Delta \overline{D}_r) = 
+    \Phi_{br}(\Delta \overline{D}_r) =
     \begin{cases}
     -1, \quad & \overline{D}_r < \overline{D}_{threshold},\\
     k_{br} \Delta \overline{D}_r, \quad & \overline{D}_{threshold} < \overline{D}_r < \overline{D}_{eq},\\
@@ -260,7 +260,13 @@ For the two moment scheme which is based on number density and mass, it is strai
   \overline{v}_{r,\, k} = \frac{1}{M_r^k} \int_0^\infty x^k f_r(x) v(x) dx,
 \end{equation}
 ```
-where the superscript ``k`` indicates the moment number, ``k=0`` for number density and ``k=1`` for mass. The individual terminal velocity of particles is approximated by ``v_(x) = (\rho_0/\rho)^{1/2} [a_R - b_R exp(-c_R D_r)]`` where ``a_R``, ``b_R`` and ``c_R`` are three free parameters and ``D_r`` is the particle diameter. Evaluating the integral results in the following equation for terminal velocity:
+where the superscript ``k`` indicates the moment number, ``k=0`` for number density and ``k=1`` for mass.
+The individual terminal velocity of particles is approximated by
+```math
+v(x) = \left(\rho_0/\rho\right)^{\frac{1}{2}} [a_R - b_R exp(-c_R D_r)]
+```
+where ``a_R``, ``b_R`` and ``c_R`` are three free parameters and ``D_r`` is the particle diameter.
+Evaluating the integral results in the following equation for terminal velocity:
 ```math
 \begin{equation}
   \overline{v}_{r,\, k} = \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}}\left[a_R - b_R \left(1+\frac{c_R}{\lambda_r}\right)^{-(3k+1)}\right],
@@ -338,7 +344,7 @@ The two-moment parametrization of evaporation suffers from the similar numerical
   \overline{x}_r = max \left(\overline{x}_{r,\, min} , min \left(\overline{x}_{r,\, max} , \frac{\rho q_{rai} \lambda_r}{N_0}\right)\right).
 \end{equation}
 ```
-This mean mass is used for computing the evaporation rate. 
+This mean mass is used for computing the evaporation rate.
 
 The default free parameter values are:
 
@@ -349,7 +355,62 @@ The default free parameter values are:
 |``\alpha_r``                | ``159 \, m \cdot s^{-1} \cdot kg^{-\beta_r}``   |
 |``\alpha_r``                | ``0.266``                                       |
 
-## Other double-moment autoconversion and accretion schemes
+## Additional 2-moment microphysics options
+
+### Terminal Velocity
+
+[Chen2022](@cite) provides a terminal velocity parameterisation
+based on an empirical fit to a high accuracy model.
+It consideres the deformation effects of large rain drops,
+as well as size-specific air density dependence.
+The fall speed of individual raindrops $v(D)$ is parameterized as:
+```math
+\begin{equation}
+    v(D) = (\phi)^{\kappa} \Sigma_{i = 1,3} \; a_i D^{b_i} e^{-c_i*D}
+\end{equation}
+```
+where D is the diameter of the particle,
+$a_i$, $b_i$, and $c_i$ are free parameers that account for deformation at larger sizes and air density dependance,
+$\phi$ is the aspect ratio (assumed to be 1 for spherical droplets), and
+$\kappa$ is 0 (corresponding to a spherical raindrop).
+$a_i$, $b_i$, and $c_i$ are listed in the table below.
+The formula is applicable when $D > 0.1 mm$,
+$q$ refers to $q = e^{0.115231 * \rho_a}$, where  $\rho_a$ refers to air density.
+The units are: [v] = m/s, [D]=mm, [$a_i$] = $mm^{-b_i} m/s$, [$b_i$] is dimensionless, [$c_i$] = 1/mm.
+
+|  $i$  |            $a_i$                        |         $b_i$                      |   $c_i$           |
+|-------|-----------------------------------------|------------------------------------|-------------------|
+|   1   |   $ 0.044612 \; q$                      |   $2.2955 \; -0.038465 \; \rho_a$  |   $0$             |
+|   2   |   $-0.263166 \; q$                      |   $2.2955 \; -0.038465 \; \rho_a$  |   $0.184325$      |
+|   3   |   $4.7178 \; q \; (\rho_a)^{-0.47335}$  |   $1.1451 \; -0.038465 \; \rho_a$  |   $0.184325$      |
+
+Assuming the same size distribution as in [SeifertBeheng2006](@cite),
+the number- and mass-weighted mean terminal velocities are:
+```math
+\begin{equation}
+  \overline{v_k} = \frac{\int_0^\infty v(D) \, D^k \, n(D) \, dD}
+             {\int_0^\infty D^k \, n(D) \, dD}
+             = (\phi)^{\kappa} \Sigma_{i} \frac{a_i \lambda^\delta \Gamma(b_i + \delta)}
+             {(\lambda + c_i)^{b_i + \delta} \; \Gamma(\delta)}
+\end{equation}
+```
+where $\Gamma$ is the gamma function, $\delta = k + 1$, and
+$\lambda$ is the size distribution parameter.
+$\overline{v_k}$ corresponds to the number-weighted mean terminal velocity when k = 0,
+and mass-weighted mean terminal velocity when k = 3.
+
+Below, we compare the individual terminal velocity formulas for [Chen2022](@cite) and [SeifertBeheng2006](@cite).
+We also compare bulk number weighted [ND] and mass weighted [MD]
+terminal velocities for both formulas integrated over the size distribution from [SeifertBeheng2006](@cite).
+We also show the mass weighted terminal velocity from the 1-moment scheme.
+
+```@example
+include("TerminalVelocityComparisons.jl")
+```
+![](terminal_velocity_bulk_comparisons.svg)
+![](terminal_velocity_individual_raindrop.svg)
+
+### Accretion and Autoconversion
 
 The other autoconversion and accretion rates in the `Microphysics2M.jl` module are implemented after Table 1 from [Wood2005](@cite)
   and are based on the works of

docs/src/TerminalVelocityComparisons.jl

@@ -0,0 +1,172 @@
+import Plots as PL
+
+FT = Float64
+
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+const CMT = CM.CommonTypes
+const CM1 = CM.Microphysics1M
+const CM2 = CM.Microphysics2M
+const CMP = CM.Parameters
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+const param_set = cloud_microphysics_parameters(toml_dict)
+
+const rain = CMT.RainType()
+const SB2006 = CMT.SB2006Type()
+const Chen2022 = CMT.Chen2022Type()
+
+"""
+    rain_terminal_velocity_individual_C(ρ, scheme, D_r)
+
+ - `ρ`` - air density
+ - `scheme` - type for parameterization
+ - `D_r` - diameter of the raindrops
+
+ Returns the fall velocity of a raindrop using multiple gamma-function terms
+ to increase accuracy for `scheme == Chen2022Type`
+"""
+function rain_terminal_velocity_individual_Chen(
+    param_set,
+    ρ::FT,
+    scheme::CMT.Chen2022Type,
+    D_r::FT, #in m
+) where {FT <: Real}
+
+    # coefficients from Table B1 from Chen et. al. 2022
+    ρ0::FT = CMP.q_coeff_rain_Ch2022(param_set)
+    a1::FT = CMP.a1_coeff_rain_Ch2022(param_set)
+    a2::FT = CMP.a2_coeff_rain_Ch2022(param_set)
+    a3::FT = CMP.a3_coeff_rain_Ch2022(param_set)
+    a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(param_set)
+    b1::FT = CMP.b1_coeff_rain_Ch2022(param_set)
+    b2::FT = CMP.b2_coeff_rain_Ch2022(param_set)
+    b3::FT = CMP.b3_coeff_rain_Ch2022(param_set)
+    b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(param_set)
+    c1::FT = CMP.c1_coeff_rain_Ch2022(param_set)
+    c2::FT = CMP.c2_coeff_rain_Ch2022(param_set)
+    c3::FT = CMP.c3_coeff_rain_Ch2022(param_set)
+
+    q = exp(ρ0 * ρ)
+    ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow)
+    bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ)
+    ci = (c1, c2, c3)
+
+    D_r = D_r * 1000 #D_r is in mm in the paper --> multiply D_r by 1000
+
+    v = 0
+    for i in 1:3
+        v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i]))
+    end
+    return v
+end
+
+function rain_terminal_velocity_individual_SB(
+    param_set,
+    ρ::FT,
+    D_r::FT,
+) where {FT <: Real}
+
+    a_r::FT = CMP.aR_tv_SB2006(param_set)
+    b_r::FT = CMP.bR_tv_SB2006(param_set)
+    c_r::FT = CMP.cR_tv_SB2006(param_set)
+    ρ_air_ground::FT = CMP.ρ0_SB2006(param_set)
+
+    v = (ρ_air_ground / ρ)^(1 / 2) * (a_r - b_r * exp(-c_r * D_r))
+    return v
+end
+
+q_liq_range = range(1e-8, stop = 1e-3, length = 1000)
+q_rai_range = range(1e-8, stop = 1e-3, length = 1000)
+D_r_range = range(50e-6, stop = 9e-3, length = 1000)
+N_d_range = range(1e7, stop = 1e9, length = 1000)
+q_liq = 5e-4
+# q_rai = 5e-4
+ρ_air = 1.0 # kg m^-3
+N_rai = 1e4 #this value matters for the individual terminal velocity
+
+PL.plot(
+    q_rai_range * 1e3,
+    [
+        CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[1]
+        for q_rai in q_rai_range
+    ],
+    linewidth = 3,
+    xlabel = "q_rain [g/kg]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-SB2006 [ND]",
+)
+PL.plot!(
+    q_rai_range * 1e3,
+    [
+        CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[1]
+        for q_rai in q_rai_range
+    ],
+    linewidth = 3,
+    xlabel = "q_rain [g/kg]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-Chen2022 [ND]",
+)
+PL.plot!(
+    q_rai_range * 1e3,
+    [
+        CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[2]
+        for q_rai in q_rai_range
+    ],
+    linewidth = 3,
+    xlabel = "q_rain [g/kg]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-SB2006 [M]",
+)
+PL.plot!(
+    q_rai_range * 1e3,
+    [
+        CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[2]
+        for q_rai in q_rai_range
+    ],
+    linewidth = 3,
+    xlabel = "q_rain [g/kg]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-Chen2022 [M]",
+)
+PL.plot!(
+    q_rai_range * 1e3,
+    [
+        CM1.terminal_velocity(param_set, rain, ρ_air, q_rai) for
+        q_rai in q_rai_range
+    ],
+    linewidth = 3,
+    xlabel = "q_rain [g/kg]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-Ogura-1-moment",
+)
+
+PL.title!("Terminal Velocity vs. Specific Humidity of Rain", titlefontsize = 9)
+PL.savefig("terminal_velocity_bulk_comparisons.svg")
+
+PL.plot(
+    D_r_range,
+    [
+        rain_terminal_velocity_individual_SB(param_set, ρ_air, D_r) for
+        D_r in D_r_range
+    ],
+    linewidth = 3,
+    xlabel = "D [m]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-SB2006",
+)
+PL.plot!(
+    D_r_range,
+    [
+        rain_terminal_velocity_individual_Chen(param_set, ρ_air, Chen2022, D_r)
+        for D_r in D_r_range
+    ],
+    linewidth = 3,
+    xlabel = "D [m]",
+    ylabel = "terminal velocity [m/s]",
+    label = "Rain-Chen2022",
+)
+PL.title!("Terminal Velocity vs. Diameter", titlefontsize = 13)
+PL.savefig("terminal_velocity_individual_raindrop.svg")

src/CommonTypes.jl

@@ -18,6 +18,7 @@ export B1994Type
 export TC1980Type
 export LD2004Type
 export SB2006Type
+export Chen2022Type
 export AbstractAerosolType
 export ArizonaTestDustType
 export DesertDustType
@@ -113,6 +114,13 @@ The type for 2-moment precipitation formation by Seifert and Beheng (2006)
 """
 struct SB2006Type <: Abstract2MPrecipType end
 
+"""
+    Chen2022Type
+
+The type for 2-moment precipitation terminal velocity by Chen et al 2022
+"""
+struct Chen2022Type <: Abstract2MPrecipType end
+
 """
     AbstractAerosolType