Parcel model with "rainbow" ice nucleation from Frostenberg et al 2023

docs/src/API.md

@@ -102,6 +102,7 @@ HetIceNucleation.ABIFM_J
 HetIceNucleation.P3_deposition_N_i
 HetIceNucleation.P3_het_N_i
 HetIceNucleation.INP_concentration_frequency
+HetIceNucleation.INP_concentration_mean
 ```
 
 # Homogeneous ice nucleation

docs/src/IceNucleationParcel0D.md

@@ -296,7 +296,7 @@ include("../../parcel/Example_Liquid_only.jl")
 The plots below are the results of the adiabatic parcel model
   with immersion freezing, condensation growth, and deposition growth for
   both a monodisperse and gamma size distribution. Note that this has not
-  yet been validated against literature. Results are very sensitive to 
+  yet been validated against literature. Results are very sensitive to
   initial conditions.
 
 ```@example
@@ -334,3 +334,65 @@ Shown below are three separate parcel simulations for deposition nucleation,
 include("../../parcel/Example_P3_ice_nuc.jl")
 ```
 ![](P3_ice_nuc.svg)
+
+
+## Immersion Freezing Parametrization based on Frostenberg et al. 2023
+
+Here we show the parcel model results when using the parametrization of immersion freezing
+  based on [Frostenberg2023](@cite).
+The concentration of ice nucleating particles (INPC) depends only on air temperature,
+  and is based on a lognormal relative frequency distribution.
+New ice crystals are created if the INPC exceeds the existing concentration of ice crystals,
+  provided there are sufficient numbers of cloud liquid droplets to freeze.
+
+Three different implementations of this parametrization are used in the parcel model:
+- `mean` - in which INPC is equal to its mean value defined in [Frostenberg2023](@cite).
+- `random` - in which INPC is sampled randomly from the distribution defined in [Frostenberg2023](@cite).
+  The number of model time steps between sampling is set by `sampling_interval`.
+- `stochastic` - in which INPC is solved for as a stochastic process,
+  with the mean and standard deviation defined in [Frostenberg2023](@cite).
+  The inverse timescale of the process is set by ``\gamma``.
+
+The stochastic implementation is based on the equation for a generic
+  stochastic process ``x`` with Gaussian random noise (i.e. Wiener process):
+```math
+\begin{equation}
+  dx = - \gamma \; x \; dt + g dW; \;\;\;\;\;\;\;   dW = {\bf{N}}(0, dt^2),
+\end{equation}
+```
+where ``\bf{N}`` is the normal distribution.
+For constant ``\gamma`` and ``g``, ``x`` has solution
+```math
+\begin{equation}
+  x(t) = x_0 \; e^{-\gamma t} + g \; \int_0^t e^{\gamma(s-t)} dW
+\end{equation}
+```
+where ``1/\gamma`` is the assumed timescale of the process.
+We can calculate the variance ``\bf{V}`` as,
+```math
+\begin{equation}
+  {\bf{V}}(t) = g^2 \int_0^t e^{2\gamma(s-t)} ds = \frac{g^2}{2\gamma}(1 - e^{-2\gamma t}).
+\end{equation}
+```
+
+We use this process to model ``x=\log(\text{INPC})``,
+  which tends toward the mean ``\mu(T)``.
+If we denote ``\tau = 1 / \gamma`` as the process timescale and
+  ``\sigma = g / \sqrt{2 / \gamma}`` as the process uncertainty,
+  then
+```math
+\begin{equation}
+  \frac{d\log(\text{INPC})}{dt} = - \frac{\log(\text{INPC}) - μ}{\tau} + \sigma \sqrt{\frac{2}{\tau \; dt}} \; {\bf{N}}(0, 1))
+\end{equation}
+```
+
+The following plot shows resuls of the parcel model with the `mean` (black line),
+  `random` (dotted lines) and `stochastic` (solid lines) parameterization options.
+We show results for two sampling intervals ``\Delta t`` (random),
+  two process time scales ``\tau`` (stochastic),
+  and two model time steps `dt`.
+
+```@example
+include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl")
+```
+![](frostenberg_immersion_freezing.svg)

docs/src/plots/Frostenberg_fig1.jl

@@ -8,20 +8,19 @@ import CloudMicrophysics.Parameters as CMP
 FT = Float32
 ip = CMP.Frostenberg2023(FT)
 
-T_range_celsius = range(-40, stop = -2, length = 500)  # air temperature [°C]
-T_range_kelvin = T_range_celsius .+ 273.15  # air temperature [K]
+T_range = range(233, stop = 271, length = 500) # air temperature [K]
 INPC_range = 10.0 .^ (range(-5, stop = 7, length = 500)) #ice nucleating particle concentration
 
 frequency = [
     IN.INP_concentration_frequency(ip, INPC, T) > 0.015 ?
     IN.INP_concentration_frequency(ip, INPC, T) : missing for
-    INPC in INPC_range, T in T_range_kelvin
+    INPC in INPC_range, T in T_range
 ]
-mu = [exp(log(-(ip.b * T)^9 * 10^(-9))) for T in T_range_celsius]
+mu = [exp(IN.INP_concentration_mean(T)) for T in T_range]
 
 
 PL.contourf(
-    T_range_kelvin,
+    T_range,
     INPC_range,
     frequency,
     xlabel = "T [K]",
@@ -35,7 +34,7 @@ PL.contourf(
 )
 
 PL.plot!(
-    T_range_kelvin,
+    T_range,
     mu,
     label = "median INPC",
     legend = :topright,

parcel/Example_Deposition_Nucleation.jl

@@ -21,6 +21,7 @@ qᵥ = FT(3.3e-4)
 qₗ = FT(0)
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -32,7 +33,7 @@ eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ)
 
 Sₗ = FT(e / eₛ)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles

parcel/Example_Frostenberg_Immersion_Freezing.jl

@@ -0,0 +1,199 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CloudMicrophysics.HetIceNucleation as CMI
+import ClimaParams as CP
+
+import Random as RAND
+RAND.seed!(44)
+N_ensemble = 32
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+FT = Float32
+# get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+
+# Initial conditions
+ρₗ = wps.ρw
+Nₐ = FT(0)
+Nₗ = FT(500 * 1e3)
+Nᵢ = FT(0)
+r₀ = FT(1e-6)
+p₀ = FT(800 * 1e2)
+T₀ = FT(251)
+qᵥ = FT(8.1e-4)
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # ρₐ ~ 1.2
+qᵢ = FT(0)
+x_sulph = FT(0.01)
+ln_INPC₀ = FT(CMI_het.INP_concentration_mean(T₀))
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_v = TD.Parameters.R_v(tps)
+Rₐ = TD.gas_constant_air(tps, q)
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+e = eᵥ(qᵥ, p₀, Rₐ, R_v)
+Sₗ = FT(e / eₛ)
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC₀]
+
+# Simulation parameters passed into ODE solver
+w = FT(0.7)
+t_max = FT(1200)
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+DSD = "Monodisperse"
+
+# Model time step
+const_dt_range = [FT(1), FT(0.25)]
+# Every how many time steps will I sample for the "random" option
+n_dt_range_12 = [[1, 60], [4, 240]]
+# Assumed process timescale for the stochastic option
+τ_range = [1, 100]
+
+results_mean = []
+labels_mean = []
+time_mean = []
+for const_dt in const_dt_range
+    # Frostenberg_mean
+    local params = parcel_params{FT}(
+        const_dt = const_dt,
+        w = w,
+        heterogeneous = "Frostenberg_mean",
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        size_distribution = DSD,
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, params)
+    push!(time_mean, sol.t)
+    push!(results_mean, sol)
+    push!(labels_mean, "dt = " * string(const_dt))
+end
+
+results_random = []
+labels_random = []
+time_random = []
+for (const_dt, n_dt_range) in zip(const_dt_range, n_dt_range_12)
+    # Frostenberg_random with different sampling frequencies
+    for n_dt in n_dt_range
+        # creating an ensamble of solutions
+        sampling_interval = FT(n_dt * const_dt)
+        solutions_rd = []
+        for ensemble_member in range(1, length = N_ensemble)
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                w = w,
+                heterogeneous = "Frostenberg_random",
+                condensation_growth = condensation_growth,
+                deposition_growth = deposition_growth,
+                size_distribution = DSD,
+                sampling_interval = sampling_interval,
+            )
+            # solve ODE
+            sol_rd = run_parcel(IC, FT(0), t_max, params)
+            push!(solutions_rd, sol_rd)
+        end
+        mean_sol = sum(solutions_rd) / length(solutions_rd)
+        push!(time_random, solutions_rd[1].t)
+        push!(results_random, mean_sol)
+        push!(
+            labels_random,
+            "dt = " * string(const_dt) * " Δt = " * string(sampling_interval),
+        )
+    end
+end
+
+results_stoch = []
+labels_stoch = []
+time_stoch = []
+for const_dt in const_dt_range
+    # Frostenberg_stochastic with different timescales γ
+    for τ in τ_range
+        # creating an ensamble of solutions
+        γ = FT(1) / τ
+        solutions_st = []
+        for ensemble_member in range(1, length = N_ensemble)
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                w = w,
+                heterogeneous = "Frostenberg_stochastic",
+                condensation_growth = condensation_growth,
+                deposition_growth = deposition_growth,
+                size_distribution = DSD,
+                γ = γ,
+            )
+            # solve ODE
+            sol_st = run_parcel(IC, FT(0), t_max, params)
+            push!(solutions_st, sol_st)
+        end
+        mean_sol = sum(solutions_st) / length(solutions_st)
+        push!(time_stoch, solutions_st[1].t)
+        push!(results_stoch, mean_sol)
+        push!(labels_stoch, "dt = " * string(const_dt) * " τ = " * string(τ))
+    end
+end
+
+# Plotting
+fig = MK.Figure(size = (1200, 900))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "T [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]")
+ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_ice [1/cm3]")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_liq [1/cm3]")
+ax7 = MK.Axis(
+    fig[4, 1:2],
+    xlabel = "T [K]",
+    ylabel = "INPC [1/m3]",
+    yscale = log10,
+)
+
+colors = [:red, :green, :orange, :limegreen]
+colors_mean = [:black, :gray]
+
+function plot_results(sol, t, ll, cl, ls)
+    MK.lines!(
+        ax1,
+        t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        linestyle = ls,
+        color = cl,
+    )
+    MK.lines!(ax2, t / 60, sol[3, :], linestyle = ls, color = cl)
+    MK.lines!(ax3, t / 60, sol[6, :] * 1e3, linestyle = ls, color = cl)
+    MK.lines!(ax4, t / 60, sol[5, :] * 1e3, linestyle = ls, color = cl)
+    MK.lines!(ax5, t / 60, sol[9, :] * 1e-6, linestyle = ls, color = cl)
+    MK.lines!(
+        ax6,
+        t / 60,
+        sol[8, :] * 1e-6,
+        linestyle = ls,
+        color = cl,
+        label = ll,
+    )
+end
+for (sol, t, ll, cl) in zip(results_stoch, time_stoch, labels_stoch, colors)
+    ls = :solid
+    plot_results(sol, t, ll, cl, ls)
+    MK.lines!(ax7, sol[3, :], exp.(sol[11, :]), linestyle = ls, color = cl)
+end
+for (sol, t, ll, cl) in zip(results_random, time_random, labels_random, colors)
+    ls = :dot
+    plot_results(sol, t, ll, cl, ls)
+end
+for (sol, t, ll, cl) in zip(results_mean, time_mean, labels_mean, colors_mean)
+    ls = :solid
+    plot_results(sol, t, ll, cl, ls)
+    MK.lines!(
+        ax7,
+        sol[3, :],
+        exp.(CMI_het.INP_concentration_mean.(sol[3, :])),
+        linestyle = ls,
+        color = cl,
+    )
+end
+fig[1:2, 3] = MK.Legend(fig, ax6, framevisible = false)
+MK.save("frostenberg_immersion_freezing.svg", fig)

parcel/Example_Immersion_Freezing.jl

@@ -24,6 +24,7 @@ qᵥ = FT(8.1e-4)
 qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -32,7 +33,7 @@ Rₐ = TD.gas_constant_air(tps, q)
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 e = eᵥ(qᵥ, p₀, Rₐ, R_v)
 Sₗ = FT(e / eₛ)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(0.4)                                # updraft speed

parcel/Example_Jensen_et_al_2022.jl

@@ -22,6 +22,7 @@ Nᵢ = FT(0)
 T₀ = FT(190)
 cᵥ₀ = FT(5 * 1e-6)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Constants
 ρₗ = wps.ρw
@@ -41,7 +42,7 @@ Sᵢ = FT(1.55)
 Sₗ = Sᵢ / ξ(tps, T₀)
 e = Sₗ * eₛ
 p₀ = e / cᵥ₀
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(1)           # updraft speed

parcel/Example_Liquid_only.jl

@@ -26,6 +26,7 @@ r₀ = FT(8e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(273.15 + 7.0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 e = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 Sₗ = FT(1)
 md_v = (p₀ - e) / R_d / T₀
@@ -34,7 +35,7 @@ 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)
 qᵢ = FT(0)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(10)                                 # updraft speed