Nonequilibrium in the parcel model

docs/src/IceNucleationParcel0D.md

@@ -179,7 +179,9 @@ Aerosol activation is described by ([see discussion](https://clima.github.io/Clo
 
 ### Condensation growth
 
-The diffusional growth of individual cloud droplet is described by
+There are both equilibrium and non-equilibrium options for condensation. The non-equilibrium equations are described in the non-equilibrium documentation ([here](https://clima.github.io/CloudMicrophysics.jl/dev/MicrophysicsNonEq)).
+
+For the equilibrium formulation, the diffusional growth of individual cloud droplet is described by
 ([see discussion](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics1M/#Snow-autoconversion)),
 ```math
 \begin{equation}
@@ -217,7 +219,10 @@ where:
 
 ### Deposition growth
 
-Similarly, for a case of a spherical ice particle growing through water vapor deposition
+Again, there are both equilibrium and non-equilibrium options for deposition. The non-equilibrium equations are described in the non-equilibrium documentation ([here](https://clima.github.io/CloudMicrophysics.jl/dev/MicrophysicsNonEq)).
+
+For the equilibrium formulation, we use a similar process as with condensation.
+For a case of a spherical ice particle growing through water vapor deposition
 ```math
 \begin{equation}
   \frac{dm_i}{dt} = 4 \pi \, r_i \, \alpha_m \, (S_i - 1) \, G_i(T)
@@ -457,3 +462,9 @@ and `stochastic` (solid lines) parameterization options. We show results for two
   include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl")
 ```
 ![](frostenberg_immersion_freezing.svg)
+
+Below, we show how the non-equilibrium formulation is able to represent the Wegener-Bergeron-Findeisen (WBF) regime in the parcel model. These plots show the liquid supersaturation percent, ice supersaturation percent, temperature, specific humidity of water vapor `q_{vap}`, specific humidity of liquid `q_{liq}`, and specific humidity of ice `q_{ice}` over time. When the supersaturation is negative evaporation/sublimation occurs and when it is positive condensation/deposition occurs. Because the water vapor pressure in the parcel is greater than the saturation water vapor pressure of liquid but larger than that of ice, liquid evaporates while ice grows.
+
+```@example
+  include("../../parcel/Example_NonEq.jl")
+```

parcel/Example_NonEq.jl

@@ -0,0 +1,125 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CloudMicrophysics.CloudDiagnostics as CD
+import CloudMicrophysics.Parameters as CMP
+import ClimaParams as CP
+
+FT = Float64
+
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+
+
+# choosing a value for the ice relaxation timescale
+τ = FT(10)
+@info("using τ value ", τ)
+
+# Get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+# Constants
+ρₗ = wps.ρw
+ρᵢ = wps.ρi
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+
+# Initial conditions
+Nₐ = FT(0)
+Nₗ = FT(200 * 1e6)
+Nᵢ = FT(1e6)
+r₀ = FT(8e-6)
+p₀ = FT(800 * 1e2)
+T₀ = FT(243)
+ln_INPC = FT(0)
+e_sat = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+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)
+
+override_file = Dict(
+    "condensation_evaporation_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
+
+liquid_toml_dict = CP.create_toml_dict(FT; override_file)
+liquid = CMP.CloudLiquid(liquid_toml_dict)
+
+override_file = Dict(
+    "sublimation_deposition_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
+
+ice_toml_dict = CP.create_toml_dict(FT; override_file)
+
+ice = CMP.CloudIce(ice_toml_dict)
+@info("relaxations:", liquid.τ_relax, ice.τ_relax)
+
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, ln_INPC]
+simple = false
+
+# Simulation parameters passed into ODE solver
+w = FT(1)                                 # updraft speed
+const_dt = FT(0.001)                         # model timestep
+t_max = FT(20)#FT(const_dt*1)
+size_distribution_list = ["Monodisperse"]
+
+condensation_growth = "NonEq_Condensation"
+deposition_growth = "NonEq_Deposition"
+
+# 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_liq [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_ice [g/kg]")
+
+for DSD in size_distribution_list
+    local params = parcel_params{FT}(
+        liq_size_distribution = DSD,
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        const_dt = const_dt,
+        w = w,
+        liquid = liquid,
+        ice = ice,
+    )
+
+    # solve ODE
+    local 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)
+
+
+    sol_Nₗ = sol[8, :]
+    sol_Nᵢ = sol[9, :]
+    sol_T = sol[3, :]
+    sol_p = sol[2, :]
+    sol_qᵥ = sol[4, :]
+    sol_qₗ = sol[5, :]
+    sol_qᵢ = sol[6, :]
+
+    q = TD.PhasePartition.(sol_qᵥ + sol_qₗ + sol_qᵢ, sol_qₗ, sol_qᵢ)
+
+    local q = TD.PhasePartition.(sol_qᵥ + sol_qₗ + sol_qᵢ, sol_qₗ, sol_qᵢ)
+    local ts = TD.PhaseNonEquil_pTq.(tps, sol_p, sol_T, q)
+    local ρₐ = TD.air_density.(tps, ts)
+
+end
+
+MK.save("noneq_parcel.svg", fig)
+nothing

parcel/ParcelModel.jl

@@ -28,6 +28,8 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     tps = TD.Parameters.ThermodynamicsParameters(FT)
     aap = CMP.AerosolActivationParameters(FT)
     ips = CMP.IceNucleationParameters(FT)
+    liquid = CMP.CloudLiquid(FT)
+    ice = CMP.CloudIce(FT)
     H₂SO₄ps = CMP.H2SO4SolutionParameters(FT)
     const_dt = 1
     w = FT(1)
@@ -289,6 +291,10 @@ function run_parcel(IC, t_0, t_end, pp)
         ce_params = Empty{FT}()
     elseif pp.condensation_growth == "Condensation"
         ce_params = CondParams{FT}(pp.aps, pp.tps, pp.const_dt)
+    elseif pp.condensation_growth == "NonEq_Condensation_simple"
+        ce_params = NonEqCondParams_simple{FT}(pp.tps, pp.liquid)
+    elseif pp.condensation_growth == "NonEq_Condensation"
+        ce_params = NonEqCondParams{FT}(pp.tps, pp.liquid, pp.const_dt)
     else
         throw("Unrecognized condensation growth mode")
     end
@@ -298,6 +304,10 @@ function run_parcel(IC, t_0, t_end, pp)
         ds_params = Empty{FT}()
     elseif pp.deposition_growth == "Deposition"
         ds_params = DepParams{FT}(pp.aps, pp.tps)
+    elseif pp.deposition_growth == "NonEq_Deposition_simple"
+        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)
     else
         throw("Unrecognized deposition growth mode")
     end

parcel/ParcelParameters.jl

@@ -82,7 +82,29 @@ struct CondParams{FT} <: CMP.ParametersType{FT}
     const_dt::FT
 end
 
+struct NonEqCondParams_simple{FT} <: CMP.ParametersType{FT}
+    tps::TDP.ThermodynamicsParameters{FT}
+    liquid::CMP.CloudLiquid{FT}
+end
+
+struct NonEqCondParams{FT} <: CMP.ParametersType{FT}
+    tps::TDP.ThermodynamicsParameters{FT}
+    liquid::CMP.CloudLiquid{FT}
+    dt::FT
+end
+
 struct DepParams{FT} <: CMP.ParametersType{FT}
     aps::CMP.ParametersType{FT}
     tps::TDP.ThermodynamicsParameters{FT}
 end
+
+struct NonEqDepParams_simple{FT} <: CMP.ParametersType{FT}
+    tps::TDP.ThermodynamicsParameters{FT}
+    ice::CMP.CloudIce{FT}
+end
+
+struct NonEqDepParams{FT} <: CMP.ParametersType{FT}
+    tps::TDP.ThermodynamicsParameters{FT}
+    ice::CMP.CloudIce{FT}
+    dt::FT
+end

parcel/ParcelTendencies.jl

@@ -2,10 +2,16 @@ import Thermodynamics as TD
 import CloudMicrophysics.Common as CMO
 import CloudMicrophysics.HetIceNucleation as CMI_het
 import CloudMicrophysics.HomIceNucleation as CMI_hom
+import CloudMicrophysics.MicrophysicsNonEq as MNE
 import CloudMicrophysics.Parameters as CMP
 import Distributions as DS
 import SpecialFunctions as SF
 
+# helper function to limit the tendency for noneq
+function limit(q, dt, n::Int)
+    return q / dt / n
+end
+
 function aerosol_activation(::Empty, state)
     FT = eltype(state)
     return FT(0)
@@ -223,6 +229,40 @@ function condensation(params::CondParams, PSD_liq, state, ρ_air)
     return qₗ + dqₗ_dt * const_dt > 0 ? dqₗ_dt : -qₗ / const_dt
 end
 
+function condensation(params::NonEqCondParams_simple, PSD, state, ρ_air)
+
+    FT = eltype(state)
+    (; Sₗ, T, qₗ, qᵥ, qᵢ) = state
+    (; tps, liquid) = params
+
+    q_sat_liq = max(Sₗ * qᵥ - qᵥ, 0)
+    q_sat = TD.PhasePartition(FT(0), q_sat_liq, FT(0))
+
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+
+    new_q = MNE.conv_q_vap_to_q_liq_ice(liquid, q_sat, q)
+
+    return new_q
+end
+
+function condensation(params::NonEqCondParams, PSD, state, ρ_air)
+    FT = eltype(state)
+    (; T, qₗ, qᵥ, qᵢ) = state
+
+    (; tps, liquid, dt) = params
+
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+
+    cond_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, q, ρ_air, T)
+
+    # using same limiter as ClimaAtmos for now
+    return ifelse(
+        cond_rate > FT(0),
+        min(cond_rate, limit(qᵥ, dt, 1)),
+        -min(abs(cond_rate), limit(qₗ, dt, 1)),
+    )
+end
+
 function deposition(::Empty, PSD_ice, state, ρ_air)
     FT = eltype(state)
     return FT(0)
@@ -236,3 +276,39 @@ function deposition(params::DepParams, PSD_ice, state, ρ_air)
     Gᵢ = CMO.G_func(aps, tps, T, TD.Ice())
     return 4 * FT(π) / ρ_air * (Sᵢ - 1) * Gᵢ * PSD_ice.r * Nᵢ
 end
+
+function deposition(params::NonEqDepParams_simple, PSD, state, ρ_air)
+    FT = eltype(state)
+    (; T, Sₗ, qₗ, qᵥ, qᵢ) = state
+
+    (; tps, ice) = params
+
+    Sᵢ = S_i(tps, T, Sₗ)
+    q_sat_ice = max(Sᵢ * qᵥ - qᵥ, 0)
+
+    q_sat = TD.PhasePartition(FT(0), FT(0), q_sat_ice)
+
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+
+    new_q = MNE.conv_q_vap_to_q_liq_ice(ice, q_sat, q)
+
+    return new_q
+end
+
+function deposition(params::NonEqDepParams, PSD, state, ρ_air)
+    FT = eltype(state)
+    (; T, qₗ, qᵥ, qᵢ) = state
+
+    (; tps, ice, dt) = params
+
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+
+    dep_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, q, ρ_air, T)
+
+    # using same limiter as ClimaAtmos for now
+    return ifelse(
+        dep_rate > FT(0),
+        min(dep_rate, limit(qᵥ, dt, 1)),
+        -min(abs(dep_rate), limit(qᵢ, dt, 1)),
+    )
+end