merge with main

docs/src/API.md

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

@@ -334,3 +334,83 @@ 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
+
+The parcel model also includes the parametrization of immersion freezing based on [Frostenberg2023](@cite).
+The concentration of ice nucleating particles (INPs) is randomly selected based on a lognormal relative frequency distribution, depending only on the temperature. If the random INP concentration exceeds the existing concentration of ice crystals, new ice crystals are created, such that their total concentration is equal to the new INP concentration.
+
+Three different implementations of this parametrization are used in the parcel model:
+- **Frostenberg_mean**, in which the concentration of INPs is equal to its mean value defined in [Frostenberg2023](@cite) for a given temperature,
+- **Frostenberg_random**, in which the concentration of INPs is drawn randomly from the distribution defined in [Frostenberg2023](@cite). The time between draws is set by the *drawing_interval* parameter,
+- **Frostenberg_stochastic**, in which the freezing is treated as a stochastic process, with mean defined in [Frostenberg2023](@cite) and Wiener noise. The timescale of the process is set by the `\gamma` parameter. 
+
+
+The stochastic implementation is based on the equation for a stochastic process $x$ with Gaussian random noise:
+
+```math
+\begin{equation}
+  dx = - \gamma \; x \; dt + g dW
+\end{equation}
+```
+
+```math
+\begin{equation}
+  dW = {\bf{N}}(0, dt^2)
+\end{equation}
+```
+where ``\bf{N}`` is the lognormal distribution.
+
+```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.
+```math
+\begin{equation}
+  {\bf{V}}(x) = g^2 \int_0^t e^{2\gamma(s-t)} ds = \frac{g^2}{2\gamma}(1 - e^{-2\gamma t})
+\end{equation}
+```
+where ``\bf{V}`` is the variance.
+We map this notation onto our problem:
+
+```math
+\begin{equation}
+  x = ln(INPC)
+\end{equation}
+```
+
+```math
+\begin{equation}
+  dx = -\gamma (x - \mu(T)) dt + g dW
+\end{equation}
+```
+
+```math
+\begin{equation}
+  x(t) = x e^{-\gamma t} + \mu(T) (1 - e^{-\gamma t}) + g \int_0^t e^{-\gamma (t-s)} dW
+\end{equation}
+```
+
+```math
+\begin{equation}
+  x(t) = {\bf{N}}(x_0 e^{-\gamma t} + (1-e^{-\gamma t}) \mu(T), \frac{g^2}{2\gamma} (1-e^{-2\gamma t}))
+\end{equation}
+```
+If ``\gamma >> 1``, then
+```math
+\begin{equation}
+  x(t) = {\bf{N}}(\mu(T), \frac{g^2}{2\gamma})
+\end{equation}
+```
+Hence ``g = \sigma * \sqrt{2 \gamma}``.
+
+
+The following plot shows resuls of the parcel model with Frostenberg_mean parametrization, Frostenberg_random parametrization with different drawing intervals `t_d` and Frostenberg_stochastic parametrization with different time scales `\gamma`. The timestep of the model is 1 s.
+
+```@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_Frostenberg_Immersion_Freezing.jl

@@ -0,0 +1,194 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+import Random as RAND
+random_seeds = [0, 1234, 5678, 3443]
+
+# 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 should be ρₐ
+qᵢ = FT(0)
+x_sulph = FT(0.01)
+
+# 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]
+
+# Simulation parameters passed into ODE solver
+w = FT(0.7)                                # updraft speed
+const_dt = FT(1)                           # model timestep
+t_max = FT(1200)
+aerosol = CMP.Illite(FT)
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+DSD = "Monodisperse"
+
+# Plotting
+fig = MK.Figure(resolution = (900, 700))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [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_liq")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice")
+
+colors = [:blue, :green, :darkorange]
+
+function plot_results(
+    sol_t,
+    sol,
+    variable,
+    label = "",
+    unit = "",
+    linestyle = :solid,
+    color = :black,
+)
+
+    MK.lines!(
+        ax1,
+        sol_t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+    MK.lines!(
+        ax2,
+        sol_t / 60,
+        sol[3, :],
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+    MK.lines!(
+        ax3,
+        sol_t / 60,
+        sol[6, :] * 1e3,
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+    MK.lines!(
+        ax4,
+        sol_t / 60,
+        sol[5, :] * 1e3,
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+    MK.lines!(
+        ax5,
+        sol_t / 60,
+        sol[8, :],
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+    MK.lines!(
+        ax6,
+        sol_t / 60,
+        sol[9, :],
+        label = label * string(variable) * unit,
+        linestyle = linestyle,
+        color = color,
+    )
+
+    MK.axislegend(ax2, position = :lt)
+end
+
+#Frostenberg_mean
+params = parcel_params{FT}(
+    const_dt = const_dt,
+    w = w,
+    aerosol = aerosol,
+    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)
+plot_results(sol.t, sol, "mean")
+
+
+# Frostenberg_random with different drawing frequencies
+drawing_interval_range = range(FT(1), stop = FT(5), length = 3)
+
+for (drawing_interval, color) in zip(drawing_interval_range, colors)
+
+    # creating an ensamble of solutions
+    solutions = []
+    for random_seed in random_seeds
+
+        RAND.seed!(trunc(Int, random_seed)) #set the random seed
+
+        params = parcel_params{FT}(
+            const_dt = const_dt,
+            w = w,
+            aerosol = aerosol,
+            heterogeneous = "Frostenberg_random",
+            condensation_growth = condensation_growth,
+            deposition_growth = deposition_growth,
+            size_distribution = DSD,
+            drawing_interval = drawing_interval,
+        )
+        # solve ODE
+        sol = run_parcel(IC, FT(0), t_max, params)
+        push!(solutions, sol)
+    end
+    mean_sol = sum(solutions) / length(solutions)
+    plot_results(sol.t, mean_sol, drawing_interval, "t_d=", " s", :dot, color)
+end
+
+
+# Frostenberg_stochastic with different timescales γ
+γ_range = range(FT(1), stop = FT(5), length = 3)
+
+for (γ, color) in zip(γ_range, colors)
+
+    # creating an ensamble of solutions
+    solutions = []
+    for random_seed in random_seeds
+
+        RAND.seed!(trunc(Int, random_seed)) #set the random seed
+
+        params = parcel_params{FT}(
+            const_dt = const_dt,
+            w = w,
+            aerosol = aerosol,
+            heterogeneous = "Frostenberg_stochastic",
+            condensation_growth = condensation_growth,
+            deposition_growth = deposition_growth,
+            size_distribution = DSD,
+            γ = γ,
+        )
+        # solve ODE
+        sol = run_parcel(IC, FT(0), t_max, params)
+        push!(solutions, sol)
+    end
+    mean_sol = sum(solutions) / length(solutions)
+    plot_results(sol.t, mean_sol, γ, "γ=", " s", :solid, color)
+end
+
+MK.save("frostenberg_immersion_freezing.svg", fig)

parcel/ParcelModel.jl

@@ -23,6 +23,9 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     w = FT(1)
     r_nuc = FT(0.5 * 1.e-4 * 1e-6)
     A_aer = FT(1e-9)
+    drawing_interval = FT(1)
+    γ = FT(1)
+    ip = CMP.Frostenberg2023(FT)
 end
 
 """
@@ -52,6 +55,7 @@ function parcel_model(dY, Y, p, t)
         Nₗ = clip!(Y[8]),
         Nᵢ = clip!(Y[9]),
         xS = Y[10],
+        t = t,
     )
 
     # Constants
@@ -61,7 +65,7 @@ function parcel_model(dY, Y, p, t)
     ρₗ = wps.ρw
 
     # Get the state values
-    (; Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₗ, Nᵢ) = state
+    (; Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₗ, Nᵢ, t) = state
     # Get thermodynamic parameters, phase partition and create thermo state.
     q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
     ts = TD.PhaseNonEquil_pTq(tps, p_air, T, q)
@@ -163,6 +167,7 @@ Parcel state vector (all variables are in base SI units):
  - Nₗ    - number concentration of existing water droplets
  - Nᵢ    - concentration of activated ice crystals
  - xS    - percent mass sulphuric acid
+ - t     - time in seconds since the beginning of the simulation
 
 The parcel parameters struct comes with default values that can be overwritten:
  - deposition - string with deposition ice nucleation parameterization choice ["None" (default), "MohlerAF", "MohlerRate", "ActivityBased", "P3_dep"]
@@ -177,6 +182,9 @@ The parcel parameters struct comes with default values that can be overwritten:
  - w - parcel vertical velocity [m/s]. Default value is 1 m/s
  - r_nuc - assumed size of nucleating ice crystals. Default value is 5e-11 [m]
  - A_aer - assumed surface area of ice nucleating aerosol. Default value assumes radius of 5e-11 [m]
+ - ip - parameters of INPC frequency distribution for Frostenberg parametrization of immerison freezing
+ - drawing_interval - time in seconds between random draws in Frostenberg_random parametrization of immerison freezing
+ - γ - timescale for the stochastic process in Frostenberg_stochastic parametrization of immerison freezing
 """
 function run_parcel(IC, t_0, t_end, pp)
 
@@ -216,6 +224,12 @@ function run_parcel(IC, t_0, t_end, pp)
         imm_params = ABIFM{FT}(pp.H₂SO₄ps, pp.tps, pp.aerosol, pp.A_aer)
     elseif pp.heterogeneous == "P3_het"
         imm_params = P3_het{FT}(pp.ips, pp.const_dt)
+    elseif pp.heterogeneous == "Frostenberg_random"
+        imm_params = Frostenberg_random{FT}(pp.ip, pp.drawing_interval)
+    elseif pp.heterogeneous == "Frostenberg_mean"
+        imm_params = Frostenberg_mean{FT}(pp.ip)
+    elseif pp.heterogeneous == "Frostenberg_stochastic"
+        imm_params = Frostenberg_stochastic{FT}(pp.ip, pp.γ)
     else
         throw("Unrecognized heterogeneous mode")
     end

parcel/ParcelParameters.jl

@@ -39,6 +39,20 @@ struct P3_het{FT} <: CMP.ParametersType{FT}
     const_dt::FT
 end
 
+struct Frostenberg_random{FT} <: CMP.ParametersType{FT}
+    ip::CMP.ParametersType{FT}
+    drawing_interval::FT
+end
+
+struct Frostenberg_stochastic{FT} <: CMP.ParametersType{FT}
+    ip::CMP.ParametersType{FT}
+    γ::FT
+end
+
+struct Frostenberg_mean{FT} <: CMP.ParametersType{FT}
+    ip::CMP.ParametersType{FT}
+end
+
 struct ABHOM{FT} <: CMP.ParametersType{FT}
     tps::TDP.ThermodynamicsParameters{FT}
     ips::CMP.ParametersType{FT}