From aa4c9eab4b3c1f117c3a1e67bb4781afc85f86d7 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Tue, 30 Jan 2024 15:39:01 -0800 Subject: [PATCH] adding activity based dep to parcel --- docs/src/IceNucleation.md | 12 ++ docs/src/plots/activity_based_deposition.jl | 84 ++++++++++++++ parcel/Deposition_Nucleation.jl | 121 ++++++++++++++++++++ parcel/Project.toml | 2 +- parcel/parcel.jl | 31 ++++- 5 files changed, 247 insertions(+), 3 deletions(-) create mode 100644 docs/src/plots/activity_based_deposition.jl create mode 100644 parcel/Deposition_Nucleation.jl 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..62c675862f --- /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(5e-11) # 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.ArizonaTestDust(FT) # aersol type for DustDeposition +aerosol_2 = CMP.Feldspar(FT) # aersol type for DepositionNucleation +ice_nucleation_modes_list = [["DustDeposition"], ["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 == "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 69532b37ad..14e8c2d57f 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -12,4 +12,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = "0.8" +CLIMAParameters = "0.8.4" diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 8ed506fa24..46e7e551cf 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -85,6 +85,33 @@ 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 * FT(π) * 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 + 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 + println("N_aer = ", N_aer) + println("A_aer = ", A_aer) + println("J_deposition = ", J_deposition) + end dN_act_dt_immersion = FT(0) dqi_dt_new_immers = FT(0) @@ -100,8 +127,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