0M: updraft precipitation parameters and evaporation of the precipitation flux

Project.toml

@@ -25,7 +25,7 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
 EmulatorModelsExt = ["DataFrames", "MLJ"]
 
 [compat]
-ClimaParams = "1.2"
+ClimaParams = "1.3.2"
 DataFrames = "1.6"
 DocStringExtensions = "0.8, 0.9"
 FastGaussQuadrature = "1"

docs/bibliography.bib

@@ -326,6 +326,39 @@ @article{Kaul2015
   year = {2015}
 }
 
+@book{Kessler1969,
+  title = {On the Distribution and Continuity of Water Substance in Atmospheric Circulations},
+  author = {Kessler, Edwin},
+  series = {Meteorological Monographs},
+  volume = {10},
+  number = {32},
+  publisher = {American Meteorological Society},
+  address = {Boston, MA},
+  year = {1969},
+  doi = {10.1007/978-1-935704-36-2_1}
+}
+
+@article{Tiedtke1993,
+  title = {Representation of Clouds in Large-Scale Models},
+  author = {Tiedtke, M.},
+  journal = {Monthly Weather Review},
+  volume = {121},
+  number = {11},
+  pages = {3040--3061},
+  year = {1993},
+  doi = {10.1175/1520-0493(1993)121<3040:ROCILS>2.0.CO;2}
+}
+
+@book{IFS2023,
+  title = {IFS Documentation CY48R1 - Part IV: Physical Processes},
+  author = {{ECMWF}},
+  series = {IFS Documentation},
+  number = {4},
+  publisher = {ECMWF},
+  year = {2023},
+  url = {https://www.ecmwf.int/en/publications/ifs-documentation}
+}
+
 @article{Kessler1995,
   title = {On the continuity and distribution of water substance in atmospheric circulations},
   author = {Kessler, Edwin},

docs/src/API.md

@@ -26,6 +26,7 @@ MicrophysicsNonEq.gamma_helper
 Microphysics0M
 Microphysics0M.remove_precipitation
 Microphysics0M.∂remove_precipitation_∂q_tot
+Microphysics0M.evaporate_precipitation
 ```
 
 ## Parameters

docs/src/Microphysics0M.md

@@ -53,3 +53,43 @@ We provide the derivative of the precipitation removal rate with respect to the
 The derivative is equal to ``-1/\tau_{precip}`` if the cloud condensate
   (liquid water plus ice) exceeds the threshold, and 0 otherwise.
 
+
+## Updraft precipitation parameters
+
+`Microphysics0MParams` also stores a second set of `Parameters0M` in
+  `precip_updraft`, for condensate in convective updrafts
+  (e.g., EDMF updrafts in ClimaAtmos).
+It uses the same thresholds as `precip` and the relaxation timescale
+  `precipitation_timescale_updraft_ratio` ``\times`` `precipitation_timescale`
+  from `ClimaParams.jl` (the default ratio of 1 gives the updrafts the same
+  timescale as `precip`).
+
+## Evaporation of precipitation
+
+The precipitation removed by the 0-moment scheme can evaporate
+  (or sublimate) as it falls through sub-saturated air.
+Following [Kessler1969](@cite), as used in the IFS
+  ([Tiedtke1993](@cite), [IFS2023](@cite)), the evaporation rate is
+```math
+\begin{equation}
+  E = a_p \, k_E \, \max(0, RH_c \, q_{vap}^{sat} - q_{vap})
+    \left[ \frac{\sqrt{p/p_s} \, F / a_p}{\alpha_2} \right]^{\alpha_3}
+\end{equation}
+```
+where:
+  - ``F`` is the grid-mean precipitation flux (positive downward),
+  - ``a_p`` is the fraction of the grid box covered by precipitation,
+  - ``k_E`` is the evaporation coefficient (``k_E = 0`` disables evaporation),
+  - ``RH_c`` is the critical relative humidity above which precipitation does not evaporate,
+  - ``q_{vap}^{sat}`` is the saturation specific humidity
+    (over liquid for rain and over ice for snow),
+  - ``p/p_s`` is the ratio of the air pressure to the surface pressure,
+  - ``\alpha_2``, ``\alpha_3`` are the flux scale and exponent.
+The parameters are stored in `Evaporation0M` and defined in `ClimaParams.jl`.
+The default ``k_E = 0`` turns evaporation off; the IFS uses
+  ``k_E = 5.44 \times 10^{-4} \, s^{-1}``.
+
+The rate is not limited.
+Time-stepping schemes should limit it so that no more precipitation
+  evaporates than falls into a grid cell, and so that the air is not
+  moistened beyond ``RH_c``.

src/Microphysics0M.jl

@@ -3,7 +3,8 @@
 
 Zero-moment bulk microphysics scheme that removes cloud condensate 
 above a threshold, equivalent to instantaneous conversion to precipitation 
-with infinite terminal velocity.
+with infinite terminal velocity. The resulting precipitation flux can
+evaporate in sub-saturated air below cloud.
 """
 module Microphysics0M
 
@@ -13,6 +14,7 @@ import CloudMicrophysics.ThermodynamicsInterface as TDI
 
 export remove_precipitation
 export ∂remove_precipitation_∂q_tot
+export evaporate_precipitation
 
 """
     remove_precipitation(params_0M::Parameters0M, q_lcl, q_icl)
@@ -74,4 +76,33 @@ The derivative is −1/τ_precip whenever precipitation is actively being remove
     q_vap_sat,
 ) = ifelse(q_lcl + q_icl > S_0 * q_vap_sat, -1 / τ_precip, zero(q_lcl))
 
+"""
+    evaporate_precipitation(params::Evaporation0M, deficit, p_over_ps, F)
+
+Returns the evaporation (or sublimation) rate of the precipitation flux
+in sub-saturated air (kg/kg/s, positive for a vapor source),
+following [Kessler1969](@cite) as used in the IFS ([Tiedtke1993](@cite)):
+
+```math
+E = a_p k_E \\max(0, \\mathrm{deficit})
+    \\left[\\frac{\\sqrt{p/p_s} \\, F / a_p}{α_2}\\right]^{α_3}
+```
+
+The rate is not limited by the available precipitation or by the vapor
+needed to saturate the air; time-stepping schemes should apply those limits.
+
+# Arguments
+- `params`: 0-moment evaporation parameters (contains `k_E`, `a_p`, `α₂`, `α₃`)
+- `deficit`: humidity deficit `RH_c q_vap_sat - q_vap` (kg/kg), with `RH_c` from
+  `params` and `q_vap_sat` over liquid for rain and over ice for snow
+- `p_over_ps`: ratio of air pressure to surface pressure (≥ 0) (-)
+- `F`: grid-mean precipitation flux, positive downward (≥ 0) (kg/m2/s)
+"""
+evaporate_precipitation(
+    (; k_E, a_p, α₂, α₃)::CMP.Evaporation0M,
+    deficit,
+    p_over_ps,
+    F,
+) = a_p * k_E * max(0, deficit) * (sqrt(p_over_ps) * F / (a_p * α₂))^α₃
+
 end #module Microphysics0M.jl

src/parameters/Microphysics0M.jl

@@ -1,4 +1,4 @@
-export Parameters0M
+export Parameters0M, Evaporation0M
 
 """
     Parameters0M{FT}
@@ -26,3 +26,38 @@ function Parameters0M(td::CP.ParamDict)
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     return Parameters0M(; parameters...)
 end
+
+"""
+    Evaporation0M{FT}
+
+Parameters for the evaporation (and sublimation) of the zero-moment
+precipitation flux in sub-saturated air, following
+[Kessler1969](@cite) as used in the IFS ([Tiedtke1993](@cite), [IFS2023](@cite)).
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+@kwdef struct Evaporation0M{FT} <: ParametersType
+    "evaporation coefficient (`0` disables evaporation) [1/s]"
+    k_E::FT
+    "critical relative humidity above which precipitation does not evaporate [-]"
+    RH_c::FT
+    "fraction of the grid box covered by precipitation [-]"
+    a_p::FT
+    "precipitation flux scale [kg/m2/s]"
+    α₂::FT
+    "precipitation flux exponent [-]"
+    α₃::FT
+end
+
+function Evaporation0M(td::CP.ParamDict)
+    name_map = (;
+        :precipitation_evaporation_coefficient => :k_E,
+        :precipitation_evaporation_rh_crit => :RH_c,
+        :precipitation_evaporation_area_fraction => :a_p,
+        :precipitation_evaporation_flux_scale => :α₂,
+        :precipitation_evaporation_exponent => :α₃,
+    )
+    parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
+    return Evaporation0M(; parameters...)
+end

src/parameters/Microphysics0MParams.jl

@@ -1,7 +1,7 @@
 export Microphysics0MParams
 
 """
-    Microphysics0MParams{P}
+    Microphysics0MParams{P, PU, E}
 
 Unified parameter container for 0-moment microphysics.
 
@@ -12,14 +12,38 @@ $(DocStringExtensions.FIELDS)
 
     Microphysics0MParams(FT)
     Microphysics0MParams(toml_dict::ClimaParams.ParamDict)
+    Microphysics0MParams(; precip, precip_updraft = precip, evaporation = nothing)
 
 where
 - `FT::Type`: A floating point type, e.g. `Float32`
 - `toml_dict::ParamDict`: A ClimaParams parameter TOML dictionary (`ParamDict`)
+
+The TOML constructor builds `precip_updraft` with the same thresholds as `precip`
+and the timescale `precipitation_timescale_updraft_ratio * precipitation_timescale`.
 """
-@kwdef struct Microphysics0MParams{P} <: ParametersType
+@kwdef struct Microphysics0MParams{P, PU, E} <: ParametersType
     "Precipitation removal parameters"
     precip::P
+    "Precipitation removal parameters for convective updrafts"
+    precip_updraft::PU = precip
+    "Precipitation evaporation parameters (`nothing` if not provided)"
+    evaporation::E = nothing
+end
+function Microphysics0MParams(toml_dict::CP.ParamDict)
+    precip = Parameters0M(toml_dict)
+    (; precipitation_timescale_updraft_ratio) = CP.get_parameter_values(
+        toml_dict,
+        "precipitation_timescale_updraft_ratio",
+        "CloudMicrophysics",
+    )
+    precip_updraft = Parameters0M(;
+        precip.qc_0,
+        precip.S_0,
+        τ_precip = precipitation_timescale_updraft_ratio * precip.τ_precip,
+    )
+    return Microphysics0MParams(;
+        precip,
+        precip_updraft,
+        evaporation = Evaporation0M(toml_dict),
+    )
 end
-Microphysics0MParams(toml_dict::CP.ParamDict) =
-    Microphysics0MParams(; precip = Parameters0M(toml_dict))

test/common_types_tests.jl

@@ -32,6 +32,7 @@ function test_common_types_broadcasts(FT)
         CMP.Chen2022VelTypeLargeIce{FT},
         CMP.Chen2022VelTypeRain{FT},
         CMP.Parameters0M{FT},
+        CMP.Evaporation0M{FT},
         CMP.ParticlePDFSnow{FT},
         CMP.ParticlePDFIceRain{FT},
         CMP.ParticleMass{FT},

test/microphysics0M_tests.jl

@@ -1,6 +1,6 @@
 import Test as TT
 
-import ClimaParams
+import ClimaParams as CP
 
 import CloudMicrophysics.ThermodynamicsInterface as TDI
 import CloudMicrophysics.Parameters as CMP
@@ -48,7 +48,96 @@ function test_microphysics0M(FT)
     end
 end
 
+function test_microphysics0M_params(FT)
+
+    TT.@testset "0M_parameters" begin
+        toml_dict = CP.create_toml_dict(FT)
+        mp = CMP.Microphysics0MParams(toml_dict)
+
+        TT.@test mp.precip isa CMP.Parameters0M{FT}
+        TT.@test mp.precip_updraft isa CMP.Parameters0M{FT}
+        TT.@test mp.evaporation isa CMP.Evaporation0M{FT}
+        TT.@test mp.precip == CMP.Parameters0M(toml_dict)
+        # default ratio of 1: updrafts use the grid-mean timescale
+        TT.@test mp.precip_updraft.τ_precip == mp.precip.τ_precip
+        TT.@test mp.precip_updraft.qc_0 == mp.precip.qc_0
+        TT.@test mp.precip_updraft.S_0 == mp.precip.S_0
+        TT.@test mp.evaporation == CMP.Evaporation0M(toml_dict)
+        TT.@test CMP.Microphysics0MParams(FT) == mp
+        TT.@test eltype(mp) == FT
+
+        # updraft timescale is the ratio times the grid-mean timescale
+        override = CP.create_toml_dict(
+            FT;
+            override_file = Dict(
+                "precipitation_timescale_updraft_ratio" =>
+                    Dict("value" => 0.2, "type" => "float"),
+                "precipitation_evaporation_coefficient" =>
+                    Dict("value" => 1.36e-4, "type" => "float"),
+            ),
+        )
+        mp_o = CMP.Microphysics0MParams(override)
+        TT.@test mp_o.precip_updraft.τ_precip ≈ FT(0.2) * mp.precip.τ_precip
+        TT.@test mp_o.precip_updraft.qc_0 == mp.precip.qc_0
+        TT.@test mp_o.precip_updraft.S_0 == mp.precip.S_0
+        TT.@test mp_o.precip.τ_precip == mp.precip.τ_precip
+        TT.@test mp_o.evaporation.k_E == FT(1.36e-4)
+
+        # keyword construction without updraft or evaporation parameters
+        mp_kw = CMP.Microphysics0MParams(; precip = mp.precip)
+        TT.@test mp_kw.precip_updraft === mp.precip
+        TT.@test isnothing(mp_kw.evaporation)
+    end
+end
+
+function test_microphysics0M_evaporation(FT)
+
+    TT.@testset "0M_evaporation" begin
+        pe = CMP.Evaporation0M(FT)
+        pe_on = CMP.Evaporation0M(;
+            k_E = FT(5.44e-4),
+            RH_c = pe.RH_c,
+            a_p = pe.a_p,
+            α₂ = pe.α₂,
+            α₃ = pe.α₃,
+        )
+        (; k_E, a_p, α₂, α₃) = pe_on
+
+        deficit = FT(2e-3)
+        p_over_ps = FT(0.8)
+        F = FT(1e-3)
+
+        # off by default
+        TT.@test pe.k_E == FT(0)
+        TT.@test CM0.evaporate_precipitation(pe, deficit, p_over_ps, F) == FT(0)
+
+        # no evaporation without a deficit or a flux
+        TT.@test CM0.evaporate_precipitation(pe_on, FT(0), p_over_ps, F) == FT(0)
+        TT.@test CM0.evaporate_precipitation(pe_on, -deficit, p_over_ps, F) == FT(0)
+        TT.@test CM0.evaporate_precipitation(pe_on, deficit, p_over_ps, FT(0)) == FT(0)
+
+        # rate value
+        E = CM0.evaporate_precipitation(pe_on, deficit, p_over_ps, F)
+        TT.@test E ≈ a_p * k_E * deficit * (sqrt(p_over_ps) * F / (a_p * α₂))^α₃
+        TT.@test E > FT(0)
+        TT.@test E isa FT
+
+        # scaling: linear in k_E and deficit, F^α₃, (p/p_s)^(α₃/2)
+        TT.@test CM0.evaporate_precipitation(pe_on, 2 * deficit, p_over_ps, F) ≈ 2 * E
+        TT.@test CM0.evaporate_precipitation(pe_on, deficit, p_over_ps, 2 * F) ≈
+                 2^α₃ * E
+        TT.@test CM0.evaporate_precipitation(pe_on, deficit, p_over_ps / 4, F) ≈
+                 (FT(1) / 2)^α₃ * E
+
+        # type stability
+        TT.@inferred CM0.evaporate_precipitation(pe_on, deficit, p_over_ps, F)
+        TT.@test @allocated(CM0.evaporate_precipitation(pe_on, deficit, p_over_ps, F)) == 0
+    end
+end
+
 TT.@testset "Microphysics 0M Tests ($FT)" for FT in (Float64, Float32)
     test_microphysics0M(FT)
+    test_microphysics0M_params(FT)
+    test_microphysics0M_evaporation(FT)
 end
 nothing