Add options for NonEq relaxation timescales as functions of N and T

docs/src/API.md

@@ -10,6 +10,8 @@ CurrentModule = CloudMicrophysics
 MicrophysicsNonEq
 MicrophysicsNonEq.τ_relax
 MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl
+MicrophysicsNonEq.conv_q_vap_to_q_lcl_MM2015
+MicrophysicsNonEq.conv_q_vap_to_q_icl_MM2015
 MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl_MM2015
 MicrophysicsNonEq.INP_limiter
 MicrophysicsNonEq.terminal_velocity

docs/src/MicrophysicsNonEq.md

@@ -99,7 +99,7 @@ To see this, it is necessary to use the definitions of ``\tau``, ``q_{sl}``, and
 
 ```math
 \begin{equation}
-  \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\;
+  \tau = \left( 4 \pi N_{tot} D_v \bar{r} \right)^{-1}, \;\;\;\;\;\;\;\;
   q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\;
   D_v = \frac{K}{\rho c_p}.
 \end{equation}
@@ -127,6 +127,41 @@ include("plots/NonEqCondEvapRate.jl")
 
 
 
+### Liquid and ice relaxation timescales as functions of number concentration
+
+We also provide functionality to set the relaxation timescales of liquid and ice, ``\tau_l`` and ``\tau_i``, as functions of cloud droplet number concentrations ``N_l`` and ``N_i`` following [MorrisonGrabowski2008_supersat](@cite) and [MorrisonMilbrandt2015](@cite). First, we approximate average radius of droplets from mass and number:
+
+```math
+\begin{equation}
+\bar{r} = \left( \frac{3q_c}{4 π N(T) \rho_c} \right)^{1/3} 
+\end{equation}
+```
+
+where ``q_c`` is the mass of condensate and ``ρ_c`` is the density of the condensate (either liquid or water). Then, we calculate the relaxation timescale:
+
+```math
+\begin{equation}
+\tau = \left( 4 \pi N_c D_v \bar{r} \right)^{-1}
+\end{equation}
+```
+
+where ``D_v`` is the thermal diffusivity and ``N_c`` is the droplet number concentration of the condensate, either ``N_l`` or ``N_i``.
+
+### Ice relaxation timescale as function of temperature through Frostenberg 2023
+
+We can also calculate ``\tau_i`` as a function of temperature using the [Frostenberg2023](@cite) parameterization. Namely, we approximate the concentration of ice crystals as equal to the mean amount of ice nucleating particles for a given temperature using [Frostenberg2023](@cite): 
+```math
+N_i(T) = ln(-(b \cdot T)^9 \times 10^{-9})
+```
+
+ ``\sigma^2`` is the variance, ``a`` and ``b`` are coefficients. The parameters defined in [Frostenberg2023](@cite) for marine data sets are ``\sigma=1.37``, ``a=1`` m``^3``, ``b=1``/C.
+
+Then, we use ``N_i`` to solve for ``\tau_i`` as demonstrated above. Here we show an example of how ``\tau_i`` changes with temperature for a set ``q_{icl}=1e-6``:
+```@example
+include("plots/NonEqFrostenberg.jl")
+```
+![](taui_Frostenberg.svg)
+
 
 ## Cloud condensate sedimentation
 

docs/src/plots/NonEqCondEvapRate.jl

@@ -50,10 +50,8 @@ function generate_cond_evap_rate(::HydrostaticBalance_qₗ_z)
     ρ = 1; title *= "ρ=$(ρ)kg/m³, "
     dt = 1; title *= "dt=$(dt)s, "
     τ_relax = 1; title *= "τ_relax=$(τ_relax)s, "
-    cm_params = CM.Parameters.CloudLiquid(FT)
-    cm_params = CM.Parameters.CloudLiquid(FT(τ_relax), cm_params.ρw, cm_params.r_eff, cm_params.N_0) # overwrite τ_relax
 
-    data = @. CMNe.conv_q_vap_to_q_lcl_icl_MM2015(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T)
+    data = @. CMNe.conv_q_vap_to_q_lcl_MM2015(thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, τ_relax)
     colorrange = extrema(data)
     @. data[iszero(data)] = NaN  # set zero values to NaN, then use `nan_color=:gray` to show them as gray. These are clipped values.
     colormap = spliced_cmap(:blues, :reds, colorrange...; mid = 0, categorical = true, symmetrize_color_ranges = true)
@@ -85,7 +83,7 @@ function generate_cond_evap_rate(::Range_qₗ_T)
     cm_params = CM.Parameters.CloudLiquid(FT)
     cm_params = CM.Parameters.CloudLiquid(FT(τ_relax), cm_params.ρw, cm_params.r_eff, cm_params.N_0) # overwrite τ_relax
 
-    data = @. CMNe.conv_q_vap_to_q_lcl_icl_MM2015(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T')
+    data = @. CMNe.conv_q_vap_to_q_lcl_MM2015(thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T', τ_relax)
     colorrange = extrema(data)
     @. data[iszero(data)] = NaN  # set zero values to NaN, then use `nan_color=:gray` to show them as gray. These are clipped values.
     colormap = spliced_cmap(:blues, :reds, colorrange...; mid = 0, categorical = true, symmetrize_color_ranges = true)

docs/src/plots/NonEqFrostenberg.jl

@@ -0,0 +1,35 @@
+import Plots as PL
+import ClimaParams as CP
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.MicrophysicsNonEq as CMNe
+
+FT = Float32
+
+# set constants
+ip = CMP.Frostenberg2023(FT)
+aps = CMP.AirProperties(FT)
+q_icl = FT(1e-6)
+
+A = FT(1e-5)
+override_file = Dict(
+    "sublimation_deposition_timescale" =>
+        Dict("value" => A, "type" => "float"),
+)
+override_toml_dict = CP.create_toml_dict(FT; override_file)
+ice = CMP.CloudIce(override_toml_dict)
+
+T_range = range(233, stop = 271, length = 500)
+
+τᵢ_values = [CMNe.τ_Frostenberg(ice, aps, ip, q_icl, T) for T in T_range]
+
+PL.plot(
+    T_range,
+    τᵢ_values,
+    xlabel = "T [K]",
+    ylabel = "τᵢ [s]",
+    label = "τᵢ(T)",
+    legend = :topright,
+    yaxis = :log,
+)
+
+PL.savefig("τᵢ_Frostenberg.svg")

parcel/Example_NonEq_Frostenberg.jl

@@ -0,0 +1,96 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+
+import CloudMicrophysics as CM
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.ThermodynamicsInterface as TDI
+import ClimaParams as CP
+
+FT = Float64
+
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+
+τₗ = FT(10)
+
+@info("using τ value ", τ)
+override_file = Dict(
+    "condensation_evaporation_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+    "sublimation_deposition_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
+override_toml_dict = CP.create_toml_dict(FT; override_file)
+liquid = CMP.CloudLiquid(override_toml_dict)
+ice = CMP.CloudIce(override_toml_dict)
+@info("relaxations:", liquid.τ_relax, ice.τ_relax)
+
+# Get free parameters
+tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+aps = CMP.AirProperties(FT)
+ip = CMP.Frostenberg2023(FT)
+
+# Constants
+ρₗ = wps.ρw
+ρᵢ = wps.ρi
+R_v = TDI.Rᵥ(tps)
+R_d = TDI.Rd(tps)
+
+# Initial conditions
+Nₐ = FT(0)
+Nₗ = FT(200 * 1e6)
+Nᵢ = FT(1e6)
+r₀ₗ = FT(1e-6)
+r₀ᵢ = FT(8e-6)
+p₀ = FT(800 * 1e2)
+T₀ = FT(243)
+ln_INPC = FT(0)
+e_sat = TDI.saturation_vapor_pressure_over_liquid(tps, T₀)
+Sₗ = FT(1)
+e = Sₗ * e_sat
+md_v = (p₀ - e) / R_d / T₀
+mv_v = e / R_v / T₀
+ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀ₗ^3
+mi_v = Nᵢ * 4 / 3 * FT(π) * ρᵢ * r₀ᵢ^3
+qᵥ = mv_v / (md_v + mv_v + ml_v + mi_v)
+qₗ = ml_v / (md_v + mv_v + ml_v + mi_v)
+qᵢ = mi_v / (md_v + mv_v + ml_v + mi_v)
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, ln_INPC]
+
+# Simulation parameters passed into ODE solver
+w = FT(1)                     # updraft speed
+const_dt = FT(0.001)          # model timestep
+t_max = FT(20)
+condensation_growth = "NonEq_Condensation"
+deposition_growth = "NonEq_Deposition_Frostenberg"
+
+# Setup the plots
+fig = MK.Figure(size = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Liquid Supersaturation [%]")
+ax2 = MK.Axis(fig[2, 1], ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "q_lcl [g/kg]")
+ax4 = MK.Axis(fig[1, 2], ylabel = "Ice Supersaturation [%]")
+ax5 = MK.Axis(fig[2, 2], ylabel = "q_vap [g/kg]")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [s]", ylabel = "q_icl [g/kg]")
+
+params = parcel_params{FT}(
+    condensation_growth = condensation_growth,
+    deposition_growth = deposition_growth,
+    const_dt = const_dt,
+    w = w,
+    liquid = liquid,
+    ice = ice,
+)
+
+# solve ODE
+sol = run_parcel(IC, FT(0), t_max, params)
+
+# Plot results
+MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0)
+MK.lines!(ax2, sol.t, sol[3, :])
+MK.lines!(ax3, sol.t, sol[5, :] * 1e3)
+MK.lines!(ax4, sol.t, (S_i.(tps, sol[3, :], sol[1, :]) .- 1) * 100.0)
+MK.lines!(ax5, sol.t, sol[4, :] * 1e3)
+MK.lines!(ax6, sol.t, sol[6, :] * 1e3)
+
+MK.save("noneq_frostenberg_parcel.svg", fig)

parcel/ParcelModel.jl

@@ -328,6 +328,8 @@ function run_parcel(IC, t_0, t_end, pp)
         ds_params = NonEqDepParams_simple{FT}(pp.tps, pp.ice)
     elseif pp.deposition_growth == "NonEq_Deposition"
         ds_params = NonEqDepParams{FT}(pp.tps, pp.ice, pp.const_dt)
+    elseif pp.deposition_growth == "NonEq_Deposition_Frostenberg"
+        ds_params = NonEqDepFrostenbergParams{FT}(pp.aps, pp.tps, pp.ip, pp.ice, pp.const_dt)
     else
         throw("Unrecognized deposition growth mode")
     end

parcel/ParcelParameters.jl

@@ -110,3 +110,11 @@ struct NonEqDepParams{FT} <: CMP.ParametersType
     ice::CMP.CloudIce{FT}
     dt::FT
 end
+
+struct NonEqDepFrostenbergParams{FT} <: CMP.ParametersType
+    aps::CMP.AirProperties{FT}
+    tps::TDP.ThermodynamicsParameters{FT}
+    ip::CMP.Frostenberg2023{FT}
+    ice::CMP.CloudIce{FT}
+    dt::FT
+end

parcel/ParcelTendencies.jl

@@ -253,7 +253,7 @@ function condensation(params::NonEqCondParams, PSD, state, ρ_air)
 
         qₜ = qᵥ + qₗ + qᵢ
 
-        cond_rate = MNE.conv_q_vap_to_q_lcl_icl_MM2015(liquid, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
+        cond_rate = MNE.conv_q_vap_to_q_lcl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, liquid.τ_relax)
 
         # Using same limiter as ClimaAtmos for now
         # Not sure why, but without intermediate storing of the tendencies for the
@@ -310,7 +310,31 @@ function deposition(params::NonEqDepParams, PSD, state, ρ_air)
     if qᵥ + qᵢ > FT(0)
         qₜ = qᵥ + qₗ + qᵢ
 
-        dep_rate = MNE.conv_q_vap_to_q_lcl_icl_MM2015(ice, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
+        dep_rate = MNE.conv_q_vap_to_q_icl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, ice.τ_relax)
+
+        # Using same limiter as ClimaAtmos for now
+        # Not sure why, but without intermediate storing of the tendencies for the
+        # if/else branch this code segfaults on julia v1.11 (works fine on v1.10)
+        dep_limit = min(dep_rate, limit(qᵥ, dt, 1))
+        sub_limit = min(abs(dep_rate), limit(qᵢ, dt, 1))
+        ret = ifelse(dep_rate > FT(0), dep_limit, -1 * sub_limit)
+        return ret
+    else
+        return FT(0)
+    end
+end
+
+function deposition(params::NonEqDepFrostenbergParams, PSD, state, ρ_air)
+    FT = eltype(state)
+    (; T, qₗ, qᵥ, qᵢ) = state
+
+    (; aps, tps, ip, ice, dt) = params
+
+    if qᵥ + qᵢ > FT(0)
+        qₜ = qᵥ + qₗ + qᵢ
+
+        τᵢ = MNE.τ_Frostenberg(ice, aps, ip, qᵢ, T)
+        dep_rate = MNE.conv_q_vap_to_q_icl_MM2015(tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T, τᵢ)
 
         # Using same limiter as ClimaAtmos for now
         # Not sure why, but without intermediate storing of the tendencies for the