CliMA · bors · Oct 13, 2023 · Sep 28, 2023
diff --git a/docs/src/IceNucleationBox.md b/docs/src/IceNucleationBox.md
@@ -199,6 +199,19 @@ where:
 where:
 - ``N_{aer}`` is the number of available dust aerosol particles.
 
+## Immersion Freezing
+Following the water activity based immersion freezing model (ABIFM), the ABIFM derived
+  nucleation rate coefficient, ``J_{immer}``, can be determined. The ice production rate,``P_{ice}``,
+  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})
+  \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.
+
 ## Example figures
 
 Here we show example simulation results from the adiabatic parcel
@@ -224,3 +237,13 @@ It is compared to [Rogers1975](@cite).
 include("../../parcel/Liquid_only.jl")
 ```
 ![](liquid_only_parcel.svg)
+
+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.
+
+```@example
+include("../../parcel/Immersion_Freezing.jl")
+```
+![](immersion_freezing.svg)
diff --git a/parcel/Immersion_Freezing.jl b/parcel/Immersion_Freezing.jl
@@ -0,0 +1,126 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CloudMicrophysics.CommonTypes as CMT
+import CLIMAParameters as CP
+
+# boilerplate code to get free parameter values
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+# Boiler plate code to have access to model parameters and constants
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+prs = cloud_microphysics_parameters(toml_dict)
+thermo_params = CMP.thermodynamics_params(prs)
+air_props = CMT.AirProperties(FT)
+
+# Constants
+ρₗ = FT(CMP.ρ_cloud_liq(prs))
+R_v = FT(CMP.R_v(prs))
+R_d = FT(CMP.R_d(prs))
+
+# Initial conditions
+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 * π * r₀^3 * ρₗ / 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ᵢ)
+ts = TD.PhaseNonEquil_pTq(thermo_params, p₀, T₀, q)
+ρₐ = TD.air_density(thermo_params, ts)
+Rₐ = TD.gas_constant_air(thermo_params, q)
+eₛ = TD.saturation_vapor_pressure(thermo_params, T₀, TD.Liquid())
+e = qᵥ * p₀ * R_v / Rₐ
+
+Sₗ = FT(e / eₛ)
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+
+# Simulation parameters passed into ODE solver
+r_nuc = FT(0.5 * 1.e-4 * 1e-6)             # assumed size of nucleated particles
+w = FT(0.7)                                # updraft speed
+α_m = FT(0.5)                              # accomodation coefficient
+const_dt = FT(1)                           # model timestep
+t_max = FT(1200)
+aerosol_type = CMT.IlliteType()
+ice_nucleation_modes = ["ImmersionFreezing"]
+growth_modes = ["Condensation", "Deposition"]
+droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]]
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+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 / N_tot",
+    yscale = log10,
+)
+MK.ylims!(ax6, 1e-6, 1)
+
+for droplet_size_distribution in droplet_size_distribution_list
+    p = (;
+        prs,
+        air_props,
+        thermo_params,
+        const_dt,
+        r_nuc,
+        w,
+        α_m,
+        aerosol_type,
+        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
+    ξ(T) =
+        TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) /
+        TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
+    S_i(T, S_liq) = ξ(T) * S_liq - 1
+
+    MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD)
+    MK.lines!(ax2, sol.t / 60, sol[3, :])
+    MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
+    MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)
+
+    sol_Nₗ = sol[8, :]
+    sol_Nᵢ = sol[9, :]
+    sol_qₗ = sol[5, :]
+    if DSD == "Monodisperse"
+        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ)
+    end
+    if DSD == "Gamma"
+        λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
+        rₗ = 2 ./ λ
+    end
+    MK.lines!(ax5, sol.t / 60, sol_Nₗ)
+    MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ))
+end
+
+MK.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 2,
+    position = :rb,
+)
+
+MK.save("immersion_freezing.svg", fig)
diff --git a/parcel/Liquid_only.jl b/parcel/Liquid_only.jl
@@ -58,8 +58,10 @@ w = FT(10)                                 # updraft speed
 α_m = FT(0.5)                              # accomodation coefficient
 const_dt = FT(0.5)                         # model timestep
 t_max = FT(20)
+aerosol_type = []
 ice_nucleation_modes = []                  # no freezing
-growth_modes_list = [["Condensation"], ["Condensation_DSD"]] # switch on condensation
+growth_modes = ["Condensation"]          # switch on condensation
+droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]]
 
 # Data from Rogers(1975) Figure 1
 # https://www.tandfonline.com/doi/abs/10.1080/00046973.1975.9648397
@@ -80,7 +82,7 @@ ax5 = MK.Axis(fig[1, 2], ylabel = "radius [μm]")
 MK.lines!(ax1, Rogers_time_supersat, Rogers_supersat, label = "Rogers_1975")
 MK.lines!(ax5, Rogers_time_radius, Rogers_radius)
 
-for growth_modes in growth_modes_list
+for droplet_size_distribution in droplet_size_distribution_list
     p = (;
         prs,
         air_props,
@@ -89,26 +91,28 @@ for growth_modes in growth_modes_list
         r_nuc,
         w,
         α_m,
+        aerosol_type,
         ice_nucleation_modes,
         growth_modes,
+        droplet_size_distribution,
     )
     # solve ODE
     sol = run_parcel(IC, FT(0), t_max, p)
 
-    gm = growth_modes[1]
+    DSD = droplet_size_distribution[1]
 
     # Plot results
-    MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = gm)
+    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[4, :] * 1e3)
     MK.lines!(ax4, sol.t, sol[5, :] * 1e3)
 
     sol_Nₗ = sol[8, :]
     sol_qₗ = sol[5, :]
-    if gm == "Condensation"
+    if DSD == "Monodisperse"
         rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ)
     end
-    if gm == "Condensation_DSD"
+    if DSD == "Gamma"
         λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
         rₗ = 2 ./ λ
     end

diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl
@@ -78,8 +78,10 @@ function Tully_et_al_2023(FT)
     w = FT(3.5 * 1e-2)                         # updraft speed
     α_m = FT(0.5)                              # accomodation coefficient
     const_dt = 0.1                             # model timestep
+    aerosol_type = CMT.DesertDustType()        # aerosol type
     ice_nucleation_modes = ["DustDeposition"]  # switch on deposition on dust
     growth_modes = ["Deposition"]              # switch on deposition growth
+    droplet_size_distribution = ["Monodisperse"]
     p = (;
         prs,
         air_props,
@@ -88,8 +90,10 @@ function Tully_et_al_2023(FT)
         r_nuc,
         w,
         α_m,
+        aerosol_type,
         ice_nucleation_modes,
         growth_modes,
+        droplet_size_distribution,
     )
 
     # Simulation 1

diff --git a/parcel/parcel.jl b/parcel/parcel.jl
@@ -15,8 +15,8 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 function parcel_model(dY, Y, p, t)
 
     # Get simulation parameters
-    (; prs, air_props, thermo_params, const_dt, r_nuc, w, α_m) = p
-    (; ice_nucleation_modes, growth_modes) = p
+    (; prs, air_props, thermo_params, const_dt, r_nuc, w, α_m, aerosol_type) = p
+    (; ice_nucleation_modes, growth_modes, droplet_size_distribution) = p
     # Numerical precision used in the simulation
     FT = eltype(Y)
 
@@ -40,6 +40,7 @@ function parcel_model(dY, Y, p, t)
     grav = CMP.grav(prs)
     ρ_ice = CMP.ρ_cloud_ice(prs)
     ρ_liq = CMP.ρ_cloud_liq(prs)
+    H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
 
     # Get thermodynamic parameters, phase partition and create thermo state.
     thermo_params = CMP.thermodynamics_params(prs)
@@ -58,6 +59,7 @@ function parcel_model(dY, Y, p, t)
     L_subl = TD.latent_heat_sublim(thermo_params, T)
     L_vap = TD.latent_heat_vapor(thermo_params, T)
     ρ_air = TD.air_density(thermo_params, ts)
+    e = q_vap * p_a * R_v / R_a
 
     # Adiabatic parcel coefficients
     a1 = L_vap * grav / cp_a / T^2 / R_v - grav / R_a / T
@@ -67,23 +69,48 @@ function parcel_model(dY, Y, p, t)
 
     # TODO - we should zero out all tendencies and augemnt them
     # TODO - add immersion, homogeneous, ...
-    #dN_act_dt_immersion = FT(0)
     #dN_act_dt_homogeneous = FT(0)
 
     dN_act_dt_depo = FT(0)
+    dqi_dt_new_depo = FT(0)
     if "DustDeposition" in ice_nucleation_modes
         AF = CMI_het.dust_activated_number_fraction(
             prs,
             S_i,
             T,
-            CMT.DesertDustType(),
+            aerosol_type,
         )
         dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt
+        dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air
+    end
+
+    dN_act_dt_immersion = FT(0)
+    dqi_dt_new_immers = FT(0)
+    if "ImmersionFreezing" in ice_nucleation_modes
+        Δa_w = T > FT(185) && T < FT(235) ?
+            CMO.a_w_xT(H2SO4_prs, thermo_params, x_sulph, T) - CMO.a_w_ice(thermo_params, T) :
+            CMO.a_w_eT(thermo_params, e, T) - CMO.a_w_ice(thermo_params, T)
+        J_immersion = CMI_het.ABIFM_J(prs, aerosol_type, Δa_w)
+        if "Monodisperse" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            A_aer = 4 * π * r_l^2
+            dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer)
+            dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+        end
+        if "Gamma" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes
+            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+            A_aer = 4 * π * r_l^2
+            dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer)
+            dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+        end
     end
-    dN_ice_dt = dN_act_dt_depo
-    dN_aer_dt = -dN_act_dt_depo
-    dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air
 
+    dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion
+    dN_aer_dt = -dN_act_dt_depo
+    dN_liq_dt = -dN_act_dt_immersion
+
     # Growth
     dqi_dt_depo = FT(0)
     if "Deposition" in growth_modes && N_ice > 0
@@ -97,24 +124,24 @@ function parcel_model(dY, Y, p, t)
 
     dql_dt_cond = FT(0)
     if "Condensation" in growth_modes && N_liq > 0
-        # Condensation on existing droplets (assuming all are the same)
-        G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid())
-        r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
-        dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq
-    end
-    if "Condensation_DSD" in growth_modes && N_liq > 0
-        # Condensation on existing droplets
-        # assuming n(r) = A r exp(-λr)
-        G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid())
-        λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
-        #A = N_liq* λ^2
-        r_l = 2 / λ
-        dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq
+        if "Monodisperse" in droplet_size_distribution
+        # Condensation on existing droplets assuming all are the same
+            G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid())
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq
+        elseif "Gamma" in droplet_size_distribution
+        # Condensation on existing droplets assuming n(r) = A r exp(-λr)
+            G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid())
+            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+            dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq
+        end
     end
 
     # Update the tendecies
-    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo
-    dq_liq_dt = dql_dt_cond
+    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers
+    dq_liq_dt = dql_dt_cond - dqi_dt_new_immers
     dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt
     dp_a_dt = -p_a * grav / R_a / T * w
     dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt
@@ -128,7 +155,7 @@ function parcel_model(dY, Y, p, t)
     dY[5] = dq_liq_dt      # liquid water specific humidity
     dY[6] = dq_ice_dt      # ice specific humidity
     dY[7] = dN_aer_dt      # number concentration of interstitial aerosol
-    dY[8] = FT(0)          # mumber concentration of droplets
+    dY[8] = dN_liq_dt      # mumber concentration of droplets
     dY[9] = dN_ice_dt      # number concentration of activated particles
     dY[10] = FT(0)         # sulphuric acid concentration
 
@@ -181,6 +208,9 @@ function run_parcel(IC, t_0, t_end, p)
     if "DustDeposition" in ice_nucleation_modes
         print("Deposition on dust particles ")
     end
+    if "ImmersionFreezing" in ice_nucleation_modes
+        print("Immersion freezing")
+    end
     print("\n")
     print("Growth modes: ")
     if "Condensation" in growth_modes