homogeneous freezing to parcel

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.15.1"
+version = "0.15.2"
 
 [deps]
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
@@ -12,7 +12,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
-CLIMAParameters = "0.7.15"
+CLIMAParameters = "0.8"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"
 RootSolvers = "0.3, 0.4"

parcel/Jensen_et_al_2022.jl

@@ -0,0 +1,158 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+import Distributions as DS
+import CloudMicrophysics.Parameters as CMP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+
+FT = Float64
+
+# Get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+aps = CMP.AirProperties(FT)
+ip = CMP.IceNucleationParameters(FT)
+H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+
+# Constants
+ρₗ = wps.ρw
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+
+# Saturation conversion equations
+ξ(T) =
+TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
+TD.saturation_vapor_pressure(tps, T, TD.Ice())
+s_i(T, s_liq) = @. s_liq * ξ(T) # saturation ratio
+
+# Initial conditions
+Nₐ = FT(300 * 1e6)
+Nₗ = FT(0)       # Jensen2022 starts with Nₐ = 300 * 1e6 and activates them within parcel
+Nᵢ = FT(0)
+r₀ = FT(25 * 1e-9)         # Value of dry aerosol r from Jensen2022
+#r₀ = FT(7 * 1e-8)          # starting when Jensen freezes
+p₀ = FT(121 * 1e2)
+T₀ = FT(190)
+#T₀ = FT(189.7)             # starting when Jensen freezes
+x_sulph = FT(0)
+
+# saturation and partial pressure
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+# Sᵢ = FT(1.55)
+# Sₗ = S_l(T₀, Sᵢ)
+# e = eₛ * Sₗ
+# ϵ = R_d / R_v
+# starting when Jensen freezes
+#q_vap = FT(3.76e-7)
+q_vap = FT(5e-6) * FT(0.6) / FT(1.2)
+#q_vap = FT(4.1e-7) # this gives more valid ICNC and ice saturation curves
+q_liq = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2
+q_ice = FT(0)
+#q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
+R_a = TD.gas_constant_air(tps, q)
+e = q_vap * p₀ * R_v / R_a 
+sₗ = e / eₛ # saturation ratio
+println("sᵢ is ", s_i(T₀, sₗ))  # saturation ratio
+
+# mass per volume for dry air, vapor and liquid
+md_v = (p₀ - e) / R_d / T₀
+mv_v = e / R_v / T₀
+ml_v = @. Nₗ * 4 / 3 * π * ρₗ * r₀^3
+
+#qᵥ = mv_v / (md_v + mv_v + ml_v)
+#qₗ = ml_v / (md_v + mv_v + ml_v)
+#qᵢ = FT(0)
+
+## Moisture dependent initial conditions
+#q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
+#ρₐ = TD.air_density(tps, ts)
+#Rₐ = TD.gas_constant_air(tps, q)
+
+#IC = [sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [sₗ, p₀, T₀, q_vap, q_liq, q_ice, Nₐ, Nₗ, Nᵢ, x_sulph] # starting when Jensen freezes
+
+println("initial conditions:")
+println("   liquid saturation: ", sₗ)
+# println("   vapor mixing ratio: ", qᵥ)
+# println("   liquid mixing ratio: ", qₗ)
+# println("   ice mixing ratio: ", qᵢ)
+println("   vapor mixing ratio: ", q_vap)   # starting when Jensen freezes
+println("   liquid mixing ratio: ", q_liq)  # starting when Jensen freezes
+
+# Simulation parameters passed into ODE solver
+r_nuc = r₀                                 # assumed size of nucleated particles
+w = FT(1)                                  # updraft speed
+α_m = FT(1)                                # accomodation coefficient
+const_dt = FT(0.01)                        # model timestep
+t_max = FT(140)
+#t_max = FT(84)  # start when Jensen freezes
+aerosol = []
+ice_nucleation_modes = ["HomogeneousFreezing"]     # homogeneous freezing only
+growth_modes = ["Deposition"]
+droplet_size_distribution_list = [["Lognormal"]]
+
+# Data from Jensen(2022) Figure 1
+# https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JD036535
+#! format: off
+Jensen_t_sat = [0, 62.71, 70.52, 76.87, 82.4, 84.84, 88.1, 92, 96.07, 100.63, 105.35, 112.51, 119.83]
+Jensen_sat = [1.55, 1.694, 1.7107, 1.7208, 1.725, 1.726, 1.7259, 1.722, 1.715, 1.702, 1.686, 1.653, 1.6126]
+Jensen_t_T = [0, 120]
+Jensen_T = [190, 189]
+Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.658, 94.123, 95.5877, 119.84]
+Jensen_ICNC = [0, 0, 0.282, 0.789, 1.804, 4.1165, 7.218, 12.12, 16.35, 16.8, 16.97, 17.086]
+#! format: on
+
+fig = MK.Figure(resolution = (1000, 1000))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation")
+ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]")
+ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]")
+ax5 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]")
+#ax6 = MK.Axis(fig[3, 2], ylabel = "r_liq [m]")
+
+MK.ylims!(ax2, 188.5, 190)
+MK.xlims!(ax2, -5, 150)
+MK.xlims!(ax3, -5, 150)
+MK.xlims!(ax4, -5, 150)
+
+MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen 2022", color = :green)
+MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green)
+MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green)
+
+for droplet_size_distribution in droplet_size_distribution_list
+    p = (;
+        wps,
+        aps,
+        tps,
+        ip,
+        const_dt,
+        r_nuc,
+        w,
+        α_m,
+        aerosol,
+        ice_nucleation_modes,
+        growth_modes,
+        droplet_size_distribution,
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, p)
+
+    DSD = droplet_size_distribution[1]
+
+    # Plot results
+    MK.lines!(ax1, sol.t, s_i(sol[3, :], (sol[1, :])), label = DSD, color = :blue)
+    MK.lines!(ax1, sol.t, (sol[1, :]), color = :red) # liq saturation
+    MK.lines!(ax2, sol.t, sol[3, :])       # temperature
+    MK.lines!(ax3, sol.t, sol[4, :] * 1e3) # q_vap
+    MK.lines!(ax4, sol.t, sol[5, :] * 1e3) # q_liq
+    MK.lines!(ax5, sol.t, sol[9, :] * 1e-6)# ICNC
+end
+
+#MK.save("Jensen_et_al_2022.svg", fig)
+fig

parcel/Project.toml

@@ -4,11 +4,9 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
-
-[compat]
-CLIMAParameters = "0.7.12"