diff --git a/Project.toml b/Project.toml index 4c314fe55..2fadf4068 100644 --- a/Project.toml +++ b/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" diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 36fdc290c..0bcb77f01 100644 --- a/docs/bibliography.bib +++ b/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}, diff --git a/docs/src/API.md b/docs/src/API.md index dde6defa4..832eeda57 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -26,6 +26,7 @@ MicrophysicsNonEq.gamma_helper Microphysics0M Microphysics0M.remove_precipitation Microphysics0M.∂remove_precipitation_∂q_tot +Microphysics0M.evaporate_precipitation ``` ## Parameters diff --git a/docs/src/Microphysics0M.md b/docs/src/Microphysics0M.md index 0de710e0c..e6d52b6b8 100644 --- a/docs/src/Microphysics0M.md +++ b/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``. diff --git a/src/Microphysics0M.jl b/src/Microphysics0M.jl index e6b8687c6..f41a6c6a5 100644 --- a/src/Microphysics0M.jl +++ b/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 diff --git a/src/parameters/Microphysics0M.jl b/src/parameters/Microphysics0M.jl index 90c9882da..d9fc566d2 100644 --- a/src/parameters/Microphysics0M.jl +++ b/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 diff --git a/src/parameters/Microphysics0MParams.jl b/src/parameters/Microphysics0MParams.jl index cff826315..80664e1f3 100644 --- a/src/parameters/Microphysics0MParams.jl +++ b/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)) diff --git a/test/common_types_tests.jl b/test/common_types_tests.jl index c7bf6eb1d..c84298bee 100644 --- a/test/common_types_tests.jl +++ b/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}, diff --git a/test/microphysics0M_tests.jl b/test/microphysics0M_tests.jl index 020ec888d..8a6c028e2 100644 --- a/test/microphysics0M_tests.jl +++ b/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