adding activity based dep to parcel

docs/src/IceNucleation.md

@@ -79,6 +79,18 @@ where ``J`` is in units of ``cm^{-2}s^{-1}``. Note that our source code returns
   ``J`` in SI units. ``m`` and ``c`` are aerosol dependent coefficients. They will
   have different values than those for ABIFM.
 
+### Water activity based deposition nucleation plot
+The following plot shows ``J`` as a function of ``\Delta a_w`` as compared to
+  figure 6 in Alpert et al 2013 and figure S5 in supplementary material of China et al 2017. Intent of this
+  plot is to prove that ``J`` is correctly parameterized as a function
+  of ``\Delta a_w``, with no implications of whether ``\Delta a_w`` is properly
+  parameterized. For more on water activity, please see above section.
+```@example
+include("plots/activity_based_deposition.jl")
+```
+![](water_activity_depo_nuc.svg)
+
+
 ## ABIFM for Sulphuric Acid Containing Droplets
 Water Activity-Based Immersion Freezing Model (ABFIM)
   is a method of parameterizing immersion freezing inspired by the time-dependent

docs/src/IceNucleationParcel0D.md

@@ -198,14 +198,28 @@ where:
 
 ## Deposition Nucleation on dust particles
 There are multiple ways of running deposition nucleation in the parcel.
-  `"MohlerAF_Deposition"` will trigger an activated fraction approach from [Mohler2006](@cite).
-  `"MohlerRate_Deposition"` will trigger a nucleation rate approach from [Mohler2006](@cite).
+  `"MohlerAF_Deposition"` will trigger an activated fraction approach
+  from [Mohler2006](@cite). `"MohlerRate_Deposition"` will trigger a
+  nucleation rate approach from [Mohler2006](@cite). For both approaches,
+  there is no nucleation if saturation over ice exceeds 1.35 as conditions
+  above this value will result in nucleation in a different mode.
+  `"ActivityBasedDeposition"` will trigger a water activity based approach
+  from [Alpert2022](@cite). In this approach, ice production rate ``P_{ice, depo}``
+  is calculated from 
+```math
+\begin{equation}
+  P_{ice} = \left[ \frac{dN_i}{dt} \right]_{dep0} = J_{depo}A_{aero}N_{aero}
+  \label{eq:ActivityBasedDeposition_P_ice}
+\end{equation}
+```
+where ``N_{areo}`` is total number of unactiviated ice nucleating particles and
+  ``A_{aero}`` is surface area of those INP.
 The deposition nucleation methods are parameterized as described in
   [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/).
 
 ## 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}``,
+  nucleation rate coefficient, ``J_{immer}``, can be determined. The ice production rate,``P_{ice, immer}``,
   per second via immersion freezing can then be calculating using
 ```math
 \begin{equation}
@@ -245,7 +259,7 @@ Between each run the water vapor specific humidity is changed,
   of the previous run.
 The prescribed vertical velocity is equal to 3.5 cm/s.
 Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines).
-The pale blue line uses the `"MohlerRate_Deposition"`approach. 
+The pale blue line uses the `"MohlerRate_Deposition"` approach. 
   We only run it for the first GCM timestep because the rate approach requires
   the change in ice saturation over time. With the discontinuous jump in saturation,
   the parameterization is unable to determine a proper nucleation rate. When we force
@@ -258,6 +272,18 @@ include("../../parcel/Tully_et_al_2023.jl")
 ```
 ![](cirrus_box.svg)
 
+The water activity based parameterization for deposition nucleation shows
+  similar outcomes when compared to the `"MohlerRate_Deposition"` approach.
+  Here, we run the parcel for 100 secs for all available aerosol types. The
+  solid lines correspond to the `"MohlerRate_Deposition"` approach while the
+  dashed lines correspond to `"ActivityBasedDeposition"`. Note that there
+  is no common aerosol type between the two parameterizations.
+
+```@example
+include("../../parcel/Deposition_Nucleation.jl")
+```
+![](deposition_nucleation.svg)
+
 In the plots below, the parcel model is ran with only condensation (no ice or freezing)
   assuming either a monodisperse or a gamma distribution of droplets.
 It is compared to [Rogers1975](@cite).

docs/src/plots/activity_based_deposition.jl

@@ -0,0 +1,87 @@
+import Plots as PL
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.HetIceNucleation as IN
+import CloudMicrophysics.Parameters as CMP
+
+FT = Float64
+const tps = TD.Parameters.ThermodynamicsParameters(FT)
+const feldspar = CMP.Feldspar(FT)
+const ferrihydrite = CMP.Ferrihydrite(FT)
+const kaolinite = CMP.Kaolinite(FT)
+
+# Initializing
+Δa_w = range(0, stop = 0.32, length = 50)    # difference in solution and ice water activity
+J_feldspar = @. IN.deposition_J(feldspar, Δa_w)      # J in SI units
+log10J_converted_feldspar = @. log10(J_feldspar * 1e-4)       # converts J into cm^-2 s^-1 and takes log
+J_ferrihydrite = @. IN.deposition_J(ferrihydrite, Δa_w)
+log10J_converted_ferrihydrite = @. log10(J_ferrihydrite * 1e-4)
+J_kaolinite = @. IN.deposition_J(kaolinite, Δa_w)
+log10J_converted_kaolinite = @. log10(J_kaolinite * 1e-4)
+
+# Alpert et al 2022 Figure 6
+# https://doi.org/10.1039/d1ea00077b
+# China et al 2017 Supplementary Figure S5
+# https://doi.org/10.1002/2016JD025817
+
+# data read from Fig 6 in Alpert et al 2022 and
+# China et al 2017 figure S5
+# using https://automeris.io/WebPlotDigitizer/
+
+#! format: off
+Alpert2022_Feldspar_Delta_a = [0.019459, 0.256216]
+Alpert2022_Feldspar_log10J = [1.039735, 4.165563]
+
+Alpert2022_Ferrihydrite_Delta_a = [0.0989189, 0.336486]
+Alpert2022_Ferrihydrite_log10J = [1.2781457, 4.21854]
+
+China2017_Delta_a = [0.01918, 0.02398, 0.02518, 0.03537, 0.07314, 0.07794, 0.10252, 0.10492, 0.1217, 0.1307, 0.15588, 0.16787, 0.20264, 0.23981, 0.256, 0.27638]
+China2017_log10J = [-4.36923, -5.07692, -1.38462, -0.64615, 1.2, 1.35385, -0.58462, 1.90769, 2.06154, 2.24615, 2.52308, 0, 2.46154, 3.32308, 4, 4.43077, 5.26154]
+#! format: on
+
+PL.plot(
+    Δa_w,
+    log10J_converted_feldspar,
+    label = "CM.jl Feldspar",
+    xlabel = "Delta a_w [unitless]",
+    ylabel = "log10(J) [cm^-2 s^-1]",
+    linestyle = :dash,
+    linecolor = :cyan,
+)
+PL.plot!(
+    Δa_w,
+    log10J_converted_ferrihydrite,
+    label = "CM.jl Ferrihydrite",
+    linestyle = :dash,
+    linecolor = :pink,
+)
+PL.plot!(
+    Δa_w,
+    log10J_converted_kaolinite,
+    label = "CM.jl Kaolinite",
+    linestyle = :dash,
+    linecolor = :orange,
+)
+PL.plot!(
+    Alpert2022_Feldspar_Delta_a,
+    Alpert2022_Feldspar_log10J,
+    linecolor = :blue,
+    label = "Alpert2022 Feldspar",
+)
+PL.plot!(
+    Alpert2022_Ferrihydrite_Delta_a,
+    Alpert2022_Ferrihydrite_log10J,
+    linecolor = :red,
+    label = "Alpert2022 Ferrihydrite",
+)
+PL.scatter!(
+    China2017_Delta_a,
+    China2017_log10J,
+    markercolor = :yellow,
+    markersize = 2,
+    label = "China2017 Kaolinite",
+)
+
+PL.savefig("water_activity_depo_nuc.svg")

parcel/Deposition_Nucleation.jl

@@ -0,0 +1,160 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# 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)
+aps = CMP.AirProperties(FT)
+wps = CMP.WaterProperties(FT)
+ip = CMP.IceNucleationParameters(FT)
+
+# Constants
+ρₗ = wps.ρw
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+
+# Initial conditions
+Nₐ = FT(2000 * 1e3)
+Nₗ = FT(0)
+Nᵢ = FT(0)
+p₀ = FT(20000)
+T₀ = FT(230)
+qᵥ = FT(3.3e-4)
+qₗ = FT(0)
+qᵢ = FT(0)
+x_sulph = 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)
+eₛ = TD.saturation_vapor_pressure(tps, 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]
+
+ξ(T) =
+    TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
+    TD.saturation_vapor_pressure(tps, T, TD.Ice())
+S_i(T, S_liq) = ξ(T) * S_liq
+
+# Simulation parameters passed into ODE solver
+r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles
+w = FT(3.5 * 1e-2)                      # updraft speed
+α_m = FT(0.5)                           # accomodation coefficient
+const_dt = FT(0.1)                      # model timestep
+t_max = FT(100)
+aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)]                     # aersol type for DustDeposition
+aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)]   # aersol type for DepositionNucleation
+ice_nucleation_modes_list =
+    [["MohlerRate_Deposition"], ["ActivityBasedDeposition"]]
+growth_modes = ["Deposition"]
+droplet_size_distribution_list = [["Monodisperse"]]
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time")
+ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time")
+
+for ice_nucleation_modes in ice_nucleation_modes_list
+    nuc_mode = ice_nucleation_modes[1]
+    droplet_size_distribution = droplet_size_distribution_list[1]
+
+    if nuc_mode == "MohlerRate_Deposition"
+        for aerosol in aerosol_1
+            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)
+
+            # Plot results
+            if aerosol == CMP.DesertDust(FT)
+                aero_label = "DesertDust"
+            elseif aerosol == CMP.ArizonaTestDust(FT)
+                aero_label = "ArizonaTestDust"
+            end
+            MK.lines!(                               # saturation
+                ax1,
+                sol.t,
+                S_i.(sol[3, :], sol[1, :]),
+                label = aero_label
+            )
+            MK.lines!(ax2, sol.t, sol[3, :])         # temperature
+            MK.lines!(ax3, sol.t, sol[6, :] * 1e3)   # q_ice
+            MK.lines!(ax4, sol.t, sol[9, :])         # N_ice
+        end
+
+    elseif nuc_mode == "ActivityBasedDeposition"
+        for aerosol in aerosol_2
+            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)
+
+            # Plot results
+            if aerosol == CMP.Feldspar(FT)
+                aero_label = "Feldspar"
+            elseif aerosol == CMP.Ferrihydrite(FT)
+                aero_label = "Ferrihydrite"
+            elseif aerosol == CMP.Kaolinite(FT)
+                aero_label = "Kaolinite"
+            end
+            MK.lines!(                                                      # saturation
+                ax1,
+                sol.t,
+                S_i.(sol[3, :], sol[1, :]),
+                label = aero_label,
+                linestyle = :dash
+            )  
+            MK.lines!(ax2, sol.t, sol[3, :], linestyle = :dash)             # temperature
+            MK.lines!(ax3, sol.t, sol[6, :] * 1e3, linestyle = :dash)       # q_ice
+            MK.lines!(ax4, sol.t, sol[9, :], linestyle = :dash)             # N_ice
+        end
+    end
+end
+
+MK.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 3,
+    position = :lt,
+)
+
+MK.save("deposition_nucleation.svg", fig)

parcel/parcel.jl

@@ -99,6 +99,15 @@ function parcel_model(dY, Y, p, t)
         end
         dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
     end
+    if "ActivityBasedDeposition" in ice_nucleation_modes
+        Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
+        J_deposition = CMI_het.deposition_J(aerosol, Δa_w)
+        if "Monodisperse" in droplet_size_distribution
+            A_aer = 4 * FT(π) * r_nuc^2
+            dN_act_dt_depo = max(FT(0), J_deposition * N_aer * A_aer)
+            dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
+        end
+    end
 
     dN_act_dt_immersion = FT(0)
     dqi_dt_new_immers = FT(0)
@@ -114,8 +123,8 @@ function parcel_model(dY, Y, p, t)
             #A = N_liq* λ^2
             r_l = 2 / λ
         end
-        A_aer = 4 * FT(π) * r_l^2
-        dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer)
+        A_l = 4 * FT(π) * r_l^2
+        dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_l)
         dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
     end
 
@@ -282,6 +291,7 @@ function run_parcel(IC, t_0, t_end, p)
 
     println(" ")
     println("Running parcel model with: ")
+
     print("Ice nucleation modes: ")
     if "MohlerAF_Deposition" in ice_nucleation_modes
         print("Deposition on dust particles using AF ")
@@ -291,13 +301,17 @@ function run_parcel(IC, t_0, t_end, p)
         print("Deposition nucleation based on rate eqn in Mohler 2006 ")
         timestepper = ODE.Euler()
     end
+    if "ActivityBasedDeposition" in ice_nucleation_modes
+        print("Water activity based deposition nucleation ")
+    end
     if "ImmersionFreezing" in ice_nucleation_modes
         print("Immersion freezing ")
     end
     if "HomogeneousFreezing" in ice_nucleation_modes
         print("Homogeneous freezing ")
     end
     print("\n")
+
     print("Growth modes: ")
     if "Condensation" in growth_modes
         print("Condensation on liquid droplets",)
@@ -309,6 +323,7 @@ function run_parcel(IC, t_0, t_end, p)
         print("Deposition on ice crystals")
     end
     print("\n")
+
     print("Size distribution modes: ")
     if "Monodisperse" in droplet_size_distribution
         print("Monodisperse size distribution")