diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index f1987dacd1..bbc147cea8 100644 --- a/docs/src/IceNucleation.md +++ b/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 diff --git a/docs/src/plots/activity_based_deposition.jl b/docs/src/plots/activity_based_deposition.jl new file mode 100644 index 0000000000..1e50dcf353 --- /dev/null +++ b/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") diff --git a/parcel/Deposition_Nucleation.jl b/parcel/Deposition_Nucleation.jl new file mode 100644 index 0000000000..ad59848ca4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/parcel/Project.toml b/parcel/Project.toml index 5fbddd66e4..32d948c706 100644 --- a/parcel/Project.toml +++ b/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" diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 8fc387460b..c74c5e66db 100644 --- a/parcel/parcel.jl +++ b/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