adding activity based dep to parcel

docs/src/IceNucleation.md

@@ -66,6 +66,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/plots/activity_based_deposition.jl

@@ -0,0 +1,84 @@
+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/
+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]
+
+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,121 @@
+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]
+
+# 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(0.1)                           # model timestep
+t_max = FT(60)
+aerosol_1 = CMP.ArizonaTestDust(FT)             # aersol type for DustDeposition
+aerosol_2 = CMP.Kaolinite(FT)              # aersol type for DepositionNucleation
+ice_nucleation_modes_list = [["DustDeposition"], ["DepositionNucleation"]]
+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 == "DustDeposition"
+        p = (;
+            wps,
+            aps,
+            tps,
+            ip,
+            const_dt,
+            r_nuc,
+            w,
+            α_m,
+            aerosol = aerosol_1,
+            ice_nucleation_modes,
+            growth_modes,
+            droplet_size_distribution,
+        )
+    else
+        p = (;
+            wps,
+            aps,
+            tps,
+            ip,
+            const_dt,
+            r_nuc,
+            w,
+            α_m,
+            aerosol = aerosol_2,
+            ice_nucleation_modes,
+            growth_modes,
+            droplet_size_distribution,
+        )
+    end
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, p)
+
+    # Plot results
+    ξ(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
+
+    MK.lines!(ax1, sol.t, S_i.(sol[3, :], sol[1, :]), label = nuc_mode) # saturation
+    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
+
+MK.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 2,
+    position = :rb,
+)
+
+#MK.save("deposition_nucleation.svg", fig)
+fig

parcel/Project.toml

@@ -11,4 +11,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
-CLIMAParameters = "0.7.12"
+CLIMAParameters = "0.8.4"

parcel/parcel.jl

@@ -81,6 +81,21 @@ function parcel_model(dY, Y, p, t)
         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
+    if "DepositionNucleation" 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
+        if "Gamma" 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
+        println("J_deposition = ", J_deposition)
+    end
 
     dN_act_dt_immersion = FT(0)
     dqi_dt_new_immers = FT(0)
@@ -89,18 +104,18 @@ function parcel_model(dY, Y, p, t)
             CMO.a_w_xT(H2SO4_prs, tps, x_sulph, T) - CMO.a_w_ice(tps, T) :
             CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
         J_immersion = CMI_het.ABIFM_J(aerosol, Δa_w)
-        if "Monodisperse" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
             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)
+            A_l = 4 * π * 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 * π * r_l^3 * ρ_ice / ρ_air
         end
-        if "Gamma" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes
+        if "Gamma" in droplet_size_distribution
             λ = 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)
+            A_l = 4 * π * 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 * π * r_l^3 * ρ_ice / ρ_air
         end
     end