Merge #95

95: Add the double moment Seifert and Beheng (2001) parametrization r=sajjadazimi a=sajjadazimi



Co-authored-by: Sajjad Azimi <[email protected]>

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.9.2"
+version = "0.10.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

docs/bibliography.bib

@@ -274,3 +274,14 @@ @article{Wood2005
   doi = {10.1175/JAS3530.1},
   year = {2005}
 }
+
+@article{SeifertBeheng2001,
+  title={A double-moment parameterization for simulating autoconversion, accretion and selfcollection},
+  author={Seifert, Axel and Beheng, Klaus D},
+  journal={Atmospheric research},
+  volume={59},
+  pages={265--281},
+  doi = {https://doi.org/10.1016/S0169-8095(01)00126-0},
+  year={2001},
+  publisher={Elsevier}
+}

docs/src/API.md

@@ -40,8 +40,12 @@ Microphysics1M.snow_melt
 # 2-moment precipitation microphysics
 ```@docs
 Microphysics2M
-Microphysics2M.conv_q_liq_to_q_rai
+Microphysics2M.LiqRaiRates
+Microphysics2M.autoconversion
 Microphysics2M.accretion
+Microphysics2M.liquid_self_collection
+Microphysics2M.rain_self_collection
+Microphysics2M.conv_q_liq_to_q_rai
 ```
 
 # Aerosol model
@@ -90,4 +94,5 @@ CommonTypes.KK2000Type
 CommonTypes.B1994Type
 CommonTypes.TC1980Type
 CommonTypes.LD2004Type
+CommonTypes.SB2001Type
 ```

src/Common.jl

@@ -51,6 +51,13 @@ function G_func(param_set::APS, T::FT, ::TD.Ice) where {FT <: Real}
     )
 end
 
+"""
+    A Heaviside step function
+"""
+function heaviside(x::FT) where {FT <: Real}
+    return FT(x > 0)
+end
+
 """
     logistic_function(x, x_0, k)
 

src/CommonTypes.jl

@@ -102,4 +102,11 @@ The type for 2-moment precipitation formation by Liu and Daum (2004)
 """
 struct LD2004Type <: Abstract2MPrecipType end
 
+"""
+    SB2001Type
+
+The type for 2-moment precipitation formation by Seifert and Beheng (2001)
+"""
+struct SB2001Type <: Abstract2MPrecipType end
+
 end #module CommoniTypes.jl

src/Microphysics1M.jl

@@ -198,7 +198,6 @@ function unpack_params(
     )
 end
 
-
 """
     lambda(q, ρ, n0, m0, me, r0, χm, Δm)
 

src/Microphysics2M.jl

@@ -1,9 +1,9 @@
 """
-    Two-moment bulk microphysics autoconversion and accretion rates from:
-      - Khairoutdinov and Kogan 2000,
-      - Beheng 1994,
-      - Tripoli and Cotton 1980,
-      - Liu and Daum 2004.
+Double-moment bulk microphysics parametrizations including:
+ - autoconversion, accretion, and self-collection rates from Seifert and Beheng 2001,
+ - additional double-moment bulk microphysics autoconversion and accretion rates
+   from: Khairoutdinov and Kogan 2000, Beheng 1994, Tripoli and Cotton 1980, and 
+   Liu and Daum 2004.
 """
 module Microphysics2M
 
@@ -13,20 +13,227 @@ const CO = Common
 import ..CommonTypes
 const CT = CommonTypes
 
+import SpecialFunctions
+const SF = SpecialFunctions
+
+import Thermodynamics
+const TD = Thermodynamics
+
 import ..Parameters
 const CMP = Parameters
 const APS = CMP.AbstractCloudMicrophysicsParameters
 
-export conv_q_liq_to_q_rai
+export autoconversion
 export accretion
+export liquid_self_collection
+export rain_self_collection
+export conv_q_liq_to_q_rai
 
 """
-    A Heaviside step function
+A structure containing the rates of change of the specific humidities and number
+densities of liquid and rain water.
 """
-function heaviside(x::FT) where {FT <: Real}
-    return FT(x > 0)
+Base.@kwdef struct LiqRaiRates{FT <: Real}
+    "Rate of change of the liquid water specific humidity"
+    dq_liq_dt::FT
+    "Rate of change of the liquid water number density"
+    dN_liq_dt::FT
+    "Rate of change of the rain water specific humidity"
+    dq_rai_dt::FT
+    "Rate of change of the rain water number density"
+    dN_rai_dt::FT
 end
 
+# Double-moment bulk microphysics autoconversion, accretion, and self-collection
+# rates from Seifert and Beheng 2001
+
+"""
+    autoconversion(param_set, scheme, q_liq, q_rai, ρ, N_liq)
+
+ - `param_set` - abstract set with Earth parameters
+ - `scheme` - type for 2-moment rain autoconversion parameterization
+ - `q_liq` - cloud water specific humidity
+ - `q_rai` - rain water specific humidity
+ - `ρ` - air density
+ - `N_liq` - cloud droplet number density
+
+Returns a LiqRaiRates object containing `q_liq`, `N_liq`, `q_rai`, `N_rai` tendencies due to
+collisions between cloud droplets (autoconversion) for `scheme == SB2001Type`
+"""
+function autoconversion(
+    param_set::APS,
+    scheme::CT.SB2001Type,
+    q_liq::FT,
+    q_rai::FT,
+    ρ::FT,
+    N_liq::FT,
+) where {FT <: Real}
+
+    if q_liq < eps(FT)
+        return LiqRaiRates(
+            dq_liq_dt = FT(0),
+            dN_liq_dt = FT(0),
+            dq_rai_dt = FT(0),
+            dN_rai_dt = FT(0),
+        )
+    end
+
+    kc::FT = CMP.kc_SB2001(param_set)
+    ν::FT = CMP.ν_SB2001(param_set)
+    xstar::FT = CMP.xstar_SB2001(param_set)
+    A::FT = CMP.A_phi_au_SB2001(param_set)
+    a::FT = CMP.a_phi_au_SB2001(param_set)
+    b::FT = CMP.b_phi_au_SB2001(param_set)
+
+    L_liq = ρ * q_liq
+    x_liq = L_liq / N_liq
+    τ = FT(1) - q_liq / (q_liq + q_rai)
+    ϕ_au = A * τ^a * (FT(1) - τ^a)^b
+
+    dL_rai_dt =
+        kc / 20 / xstar * (ν + 2) * (ν + 4) / (ν + 1)^2 *
+        L_liq^2 *
+        x_liq^2 *
+        (1 + ϕ_au / (1 - τ)^2)
+    dN_rai_dt = dL_rai_dt / xstar
+    dL_liq_dt = -dL_rai_dt
+    dN_liq_dt = -2 * dN_rai_dt
+
+    return LiqRaiRates(
+        dq_liq_dt = dL_liq_dt / ρ,
+        dN_liq_dt = dN_liq_dt,
+        dq_rai_dt = dL_rai_dt / ρ,
+        dN_rai_dt = dN_rai_dt,
+    )
+end
+
+"""
+    accretion(param_set, scheme, q_liq, q_rai, ρ)
+
+ - `param_set` - abstract set with Earth parameters
+ - `scheme` - type for 2-moment accretion parameterization
+ - `q_liq` - cloud water specific humidity
+ - `q_rai` - rain water specific humidity
+ - `ρ` - air density
+ - `N_liq` - cloud droplet number density
+
+Returns a LiqRaiRates object containing `q_liq`, `N_liq`, `q_rai`, `N_rai` tendencies due to
+collisions between raindrops and cloud droplets (accretion) for `scheme == SB2001Type`
+"""
+function accretion(
+    param_set::APS,
+    scheme::CT.SB2001Type,
+    q_liq::FT,
+    q_rai::FT,
+    ρ::FT,
+    N_liq::FT,
+) where {FT <: Real}
+
+    if q_liq < eps(FT) || q_rai < eps(FT)
+        return LiqRaiRates(
+            dq_liq_dt = FT(0),
+            dN_liq_dt = FT(0),
+            dq_rai_dt = FT(0),
+            dN_rai_dt = FT(0),
+        )
+    end
+
+    kr::FT = CMP.kr_SB2001(param_set)
+    τ0::FT = CMP.τ_0_phi_ac_SB2001(param_set)
+    c::FT = CMP.c_phi_ac_SB2001(param_set)
+
+    L_liq = ρ * q_liq
+    L_rai = ρ * q_rai
+    x_liq = L_liq / N_liq
+    τ = FT(1) - q_liq / (q_liq + q_rai)
+    ϕ_ac = (τ / (τ + τ0))^c
+
+    dL_rai_dt = kr * L_liq * L_rai * ϕ_ac
+    dN_rai_dt = FT(0)
+    dL_liq_dt = -dL_rai_dt
+    dN_liq_dt = dL_liq_dt / x_liq
+
+    return LiqRaiRates(
+        dq_liq_dt = dL_liq_dt / ρ,
+        dN_liq_dt = dN_liq_dt,
+        dq_rai_dt = dL_rai_dt / ρ,
+        dN_rai_dt = dN_rai_dt,
+    )
+end
+
+"""
+    liquid_self_collection(param_set, scheme, q_liq, ρ, dN_liq_dt_au)
+
+ - `param_set` - abstract set with Earth parameters
+ - `scheme` - type for 2-moment liquid self-collection parameterization
+ - `q_liq` - cloud water specific humidity
+ - `ρ` - air density
+ - `dN_liq_dt_au` - rate of change of cloud droplets number density due to autoconversion
+
+Returns the cloud droplets number density tendency due to collisions of cloud droplets
+that produce larger cloud droplets (self-collection) for `scheme == SB2001Type`
+"""
+function liquid_self_collection(
+    param_set::APS,
+    scheme::CT.SB2001Type,
+    q_liq::FT,
+    ρ::FT,
+    dN_liq_dt_au::FT,
+) where {FT <: Real}
+
+    if q_liq < eps(FT)
+        return FT(0)
+    end
+
+    kc::FT = CMP.kc_SB2001(param_set)
+    ν::FT = CMP.ν_SB2001(param_set)
+
+    L_liq = ρ * q_liq
+
+    dN_liq_dt_sc = -kc * (ν + 2) / (ν + 1) * L_liq^2 - dN_liq_dt_au
+
+    return dN_liq_dt_sc
+end
+
+"""
+    rain_self_collection(param_set, scheme, q_rai, ρ, N_rai)
+
+ - `param_set` - abstract set with Earth parameters
+ - `scheme` - type for 2-moment rain self-collection parameterization
+ - `q_rai` - rain water specific humidity
+ - `ρ` - air density
+ - `N_rai` - raindrops number density
+
+Returns the raindrops number density tendency due to collisions of raindrops
+that produce larger raindrops (self-collection) for `scheme == SB2001Type`
+"""
+function rain_self_collection(
+    param_set::APS,
+    scheme::CT.SB2001Type,
+    q_rai::FT,
+    ρ::FT,
+    N_rai::FT,
+) where {FT <: Real}
+
+    if q_rai < eps(FT)
+        return FT(0)
+    end
+
+    kr::FT = CMP.kr_SB2001(param_set)
+
+    L_rai = ρ * q_rai
+
+    dN_rai_dt_sc = -kr * N_rai * L_rai
+
+    return dN_rai_dt_sc
+end
+
+# Additional double moment autoconversion and accretion parametrizations:
+# - Khairoutdinov and Kogan (2000)
+# - Beheng (1994)
+# - Tripoli and Cotton (1980)
+# - Liu and Daum (2004)
+
 """
     conv_q_liq_to_q_rai(param_set, scheme, q_liq, ρ; N_d, smooth_transition)
 
@@ -38,8 +245,8 @@ end
 
 Returns the q_rai tendency due to collisions between cloud droplets
 (autoconversion), parametrized following:
- - Khairoutdinov and Kogan (2000) for `scheme == KK200Type`
- - Beheng (1994) for `scheme == Beheng1994Type`
+ - Khairoutdinov and Kogan (2000) for `scheme == KK2000Type`
+ - Beheng (1994) for `scheme == B1994Type`
  - Tripoli and Cotton (1980) for `scheme == TC1980Type`
  - Liu and Daum (2004) for `scheme ==LD2004Type`
 
@@ -126,7 +333,7 @@ function conv_q_liq_to_q_rai(
         _k::FT = CMP.k_thrshld_stpnss(param_set)
         _output = CO.logistic_function(q_liq, q_liq_threshold, _k)
     else
-        _output = heaviside(q_liq - q_liq_threshold)
+        _output = CO.heaviside(q_liq - q_liq_threshold)
     end
     return D * q_liq^a * N_d^b * _output
 end
@@ -160,7 +367,7 @@ function conv_q_liq_to_q_rai(
             _k::FT = CMP.k_thrshld_stpnss(param_set)
             _output = CO.logistic_function(R_6, R_6C, _k)
         else
-            _output = heaviside(R_6 - R_6C)
+            _output = CO.heaviside(R_6 - R_6C)
         end
         return E * (q_liq * ρ)^3 / N_d / ρ * _output
     end
@@ -173,11 +380,11 @@ end
  - `scheme` - type for 2-moment rain accretion parameterization
  - `q_liq` - cloud water specific humidity
  - `q_rai` - rain water specific humidity
- - `ρ` - air density (for `KK200Type` and `Beheng1994Type`)
+ - `ρ` - air density (for `KK2000Type` and `Beheng1994Type`)
 
  Returns the accretion rate of rain, parametrized following
- - Khairoutdinov and Kogan (2000) for `scheme == KK200Type`
- - Beheng (1994) for `scheme == Beheng1994Type`
+ - Khairoutdinov and Kogan (2000) for `scheme == KK2000Type`
+ - Beheng (1994) for `scheme == B1994Type`
  - Tripoli and Cotton (1980) for `scheme == TC1980Type`
 """
 function accretion(

src/Parameters.jl

@@ -95,6 +95,15 @@ Base.@kwdef struct CloudMicrophysicsParameters{FT, TP} <:
     R_6C_coeff_LD2004::FT
     E_0_LD2004::FT
     k_thrshld_stpnss::FT
+    kc_SB2001::FT
+    kr_SB2001::FT
+    xstar_SB2001::FT
+    ν_SB2001::FT
+    A_phi_au_SB2001::FT
+    a_phi_au_SB2001::FT
+    b_phi_au_SB2001::FT
+    τ_0_phi_ac_SB2001::FT
+    c_phi_ac_SB2001::FT
 end
 
 Base.eltype(::CloudMicrophysicsParameters{FT}) where {FT} = FT
@@ -139,5 +148,4 @@ cv_v(ps::CMPS) = TD.Parameters.cv_v(thermodynamics_params(ps))
 cv_l(ps::CMPS) = TD.Parameters.cv_l(thermodynamics_params(ps))
 cv_i(ps::CMPS) = TD.Parameters.cv_i(thermodynamics_params(ps))
 
-
 end # module