diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index ea3c24ae95..177aaa43f6 100644 --- a/docs/src/IceNucleation.md +++ b/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 diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 6be24cb69a..47c5b7368e 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/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). diff --git a/docs/src/plots/activity_based_deposition.jl b/docs/src/plots/activity_based_deposition.jl new file mode 100644 index 0000000000..276ed8d41f --- /dev/null +++ b/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") diff --git a/parcel/Deposition_Nucleation.jl b/parcel/Deposition_Nucleation.jl new file mode 100644 index 0000000000..aad3527538 --- /dev/null +++ b/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) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index c11b2b97d0..afb64ac4a8 100644 --- a/parcel/parcel.jl +++ b/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,6 +301,9 @@ 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 @@ -298,6 +311,7 @@ function run_parcel(IC, t_0, t_end, p) 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")