Nonequilibrium in the parcel model

parcel/Example_NonEq.jl

@@ -0,0 +1,131 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+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)
+liquid = CMP.CloudLiquid(τ, ρₗ)
+
+override_file = Dict(
+    "cloud_ice_apparent_density" => Dict("value" => ρᵢ, "type" => "float"),
+    "sublimation_deposition_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
+
+toml_dict = CP.create_toml_dict(FT; override_file)
+
+ice = CMP.CloudIce(toml_dict)
+@info("relaxations:", liquid.τ_relax, ice.τ_relax)
+
+# 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)
+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", "Gamma"]
+
+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, label = DSD)
+    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,
+        label = DSD,
+    )
+    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.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 2,
+    position = :rb,
+)
+
+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_Anna"
+        ce_params = NonEqCondParams_Anna{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_Anna"
+        ds_params = NonEqDepParams_Anna{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_Anna{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_Anna{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_Anna, 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, 2)),
+        -min(abs(cond_rate), limit(qₗ, dt, 2)),
+    )
+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_Anna, 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, 2)),
+        -min(abs(dep_rate), limit(qᵢ, dt, 2)),
+    )
+end