homogeneous freezing to parcel

docs/bibliography.bib

@@ -182,6 +182,20 @@ @article{Heymsfield2003
   year = {2003}
 }
 
+@article{Jensen2022,
+author = {Jensen, E. J. and Diskin, G. S. and DiGangi, J. and Woods, S. and Lawson, R. P. and Bui, T. V.},
+title = {Homogeneous Freezing Events Sampled in the Tropical Tropopause Layer},
+journal = {Journal of Geophysical Research: Atmospheres},
+volume = {127},
+number = {17},
+pages = {e2022JD036535},
+keywords = {cirrus, nucleation, freezing},
+doi = {https://doi.org/10.1029/2022JD036535},
+url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022JD036535},
+note = {e2022JD036535 2022JD036535},
+year = {2022}
+}
+
 @article{Karcher2006,
   author = {Kärcher, B. and Hendricks, J. and Lohmann, U.},
   title = {Physically based parameterization of cirrus cloud formation for use in global atmospheric models},

docs/src/IceNucleationParcel0D.md

@@ -165,7 +165,7 @@ For a gamma distribution of droplets ``n(r) = A \; r \; exp(-\lambda r)``,
 ``\bar{r} = \frac{2}{\lambda}``
 where ``\lambda = \left(\frac{32 \pi N_{tot} \rho_l}{q_l \rho_a}\right)^{1/3}``.
 
-## Deposition on dust particles
+## Deposition nucleation on dust particles
 
 Similarly, for a case of a spherical ice particle growing through water vapor deposition
 ```math
@@ -197,7 +197,7 @@ where:
 - ``N_{act}`` is the number of activated ice particles.
 
 ``N_{act}`` can be computed for example from
-  [activated fraction](https://clima.github.io/CloudMicrophysics.jl/previews/PR103/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i``
+  [activated fraction](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i``
 ```math
 \begin{equation}
   N_{act} = N_{aer} f_i
@@ -212,17 +212,35 @@ Following the water activity based immersion freezing model (ABIFM), the ABIFM d
   per second via immersion freezing can then be calculating using
 ```math
 \begin{equation}
-  P_{ice} = [\frac{dN_i}{dt}]_{immer} = J_{immer}A(N_{liq})
+  P_{ice} = \left[ \frac{dN_i}{dt} \right]_{immer} = J_{immer}A(N_{liq})
   \label{eq:ABIFM_P_ice}
 \end{equation}
 ```
 where ``N_{liq}`` is total number of ice nuclei containing droplets and 
   ``A`` is surface area of those droplets.
 
+## Homogeneous Freezing
+Homogeneous freezing follows the water-activity based model described in the
+  [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/) section which gives a nucleation rate coefficient of
+  units ``cm^{-3}s^{-1}``.
+The ice production rate from homogeneous freezing can then be determined:
+```math
+\begin{equation}
+  P_{ice} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq})
+  \label{eq:hom_P_ice}
+\end{equation}
+```
+where ``N_{liq}`` is total number of ice nuclei containing droplets and 
+  ``V`` is the volume of those droplets.
+
 ## Example figures
 
-Here we show example simulation results from the adiabatic parcel
-  model with deposition freezing on dust.
+Here we show various example simulation results from the adiabatic parcel
+  model. This includes examples with deposition nucleation on dust,
+  liquid processes only, immersion freezing with condensation and deposition growth,
+  and homogeneous freezing with deposition growth.
+
+We start with deposition freezing on dust.
 The model is run three times for 30 minutes simulation time,
   (shown by three different colors on the plot).
 Between each run the water vapor specific humidity is changed,
@@ -254,3 +272,16 @@ The plots below are the results of the adiabatic parcel model
 include("../../parcel/Immersion_Freezing.jl")
 ```
 ![](immersion_freezing.svg)
+
+The following plots show the parcel model running with homogeneous freezing and 
+  depositional growth assuming a lognormal distribution of aerosols. 
+  It is compared against [Jensen2022](@cite). Note that running with the initial
+  conditions described in [Jensen2022](@cite) results in a ``\Delta a_w`` smaller
+  than the minimum valid value for the ``J_{hom}`` parameterization. We have forced
+  the ``\Delta a_w`` to be equal to the minimum valid value in this example only
+  for demonstrative purposes.
+
+```@example
+include("../../parcel/Jensen_et_al_2022.jl")
+```
+![](Jensen_et_al_2022.svg)

parcel/Immersion_Freezing.jl

@@ -26,7 +26,7 @@ r₀ = FT(1e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
-qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ
 qᵢ = FT(0)
 x_sulph = FT(0.01)
 
@@ -102,10 +102,10 @@ for droplet_size_distribution in droplet_size_distribution_list
     sol_Nᵢ = sol[9, :]
     sol_qₗ = sol[5, :]
     if DSD == "Monodisperse"
-        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ)
+        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
     end
     if DSD == "Gamma"
-        λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
+        λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
         rₗ = 2 ./ λ
     end
     MK.lines!(ax5, sol.t / 60, sol_Nₗ)

parcel/Jensen_et_al_2022.jl

@@ -0,0 +1,159 @@
+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)
+ϵₘ = R_d / R_v
+
+# 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)
+Nᵢ = FT(0)
+r₀ = FT(25 * 1e-9)
+T₀ = FT(190)
+cᵥ₀ = FT(5 * 1e-6)
+x_sulph = FT(0)
+
+# saturation and partial pressure
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀)
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2
+qᵢ = FT(0)
+q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_a = TD.gas_constant_air(tps, q)
+Sᵢ = FT(1.55)
+Sₗ = Sᵢ / ξ(T₀)
+e = Sₗ * eₛ
+p₀ = e / cᵥ₀
+
+## Moisture dependent initial conditions
+ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
+
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+# 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(120)     # total time
+aerosol = []        # aerosol type. fill in if immersion or deposition freezing
+R_mode_liq = r_nuc  # lognormal distribution mode radius
+σ = FT(2)           # lognormal distribution std deviation
+ice_nucleation_modes = ["HomogeneousFreezing"]     # homogeneous freezing only
+growth_modes = ["Deposition"]
+droplet_size_distribution_list = [["Jensen_Example"]]
+
+# 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 = "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 = "q_ice [g/kg]")
+
+MK.ylims!(ax2, 188.5, 190)
+MK.xlims!(ax2, -5, 150)
+MK.xlims!(ax3, -5, 150)
+MK.xlims!(ax4, -5, 150)
+
+#! format: off
+MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen et al 2022", color = :green)
+MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green, label = "Jensen et al 2022")
+MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green, label = "Jensen et al 2022")
+
+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,
+        R_mode_liq,
+        σ,
+    )
+    # 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 = "ice", color = :blue)
+    MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation
+    MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl")                  # 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, label = "CM.jl")           # ICNC
+    MK.lines!(ax6, sol.t, sol[6, :] * 1e3) # q_ice
+
+    MK.axislegend(
+        ax1,
+        framevisible = false,
+        labelsize = 12,
+        orientation = :horizontal,
+        nbanks = 2,
+        position = :lc,
+    )
+    MK.axislegend(
+        ax5,
+        framevisible = false,
+        labelsize = 12,
+        orientation = :horizontal,
+        nbanks = 2,
+        position = :lc,
+    )
+    MK.axislegend(
+        ax2,
+        framevisible = false,
+        labelsize = 12,
+        orientation = :horizontal,
+        nbanks = 2,
+        position = :lc,
+    )
+#! format: on
+end
+
+MK.save("Jensen_et_al_2022.svg", fig)

parcel/Liquid_only.jl

@@ -36,7 +36,7 @@ e = eₛ
 # 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
+ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀^3
 
 qᵥ = mv_v / (md_v + mv_v + ml_v)
 qₗ = ml_v / (md_v + mv_v + ml_v)
@@ -110,10 +110,10 @@ for droplet_size_distribution in droplet_size_distribution_list
     sol_Nₗ = sol[8, :]
     sol_qₗ = sol[5, :]
     if DSD == "Monodisperse"
-        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ)
+        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
     end
     if DSD == "Gamma"
-        λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
+        λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
         rₗ = 2 ./ λ
     end
     MK.lines!(ax5, sol.t, rₗ * 1e6)

parcel/Project.toml

@@ -4,6 +4,7 @@ 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"