From 456a3e452bd022d4dd4ace8bc97d9114c8d25eb5 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Thu, 28 Sep 2023 14:54:58 -0700 Subject: [PATCH] ABIFM to parcel --- docs/src/IceNucleationBox.md | 23 +++++++ parcel/Immersion_Freezing.jl | 126 +++++++++++++++++++++++++++++++++++ parcel/Liquid_only.jl | 16 +++-- parcel/Tully_et_al_2023.jl | 4 ++ parcel/parcel.jl | 76 ++++++++++++++------- 5 files changed, 216 insertions(+), 29 deletions(-) create mode 100644 parcel/Immersion_Freezing.jl diff --git a/docs/src/IceNucleationBox.md b/docs/src/IceNucleationBox.md index 808c127712..9405b9991d 100644 --- a/docs/src/IceNucleationBox.md +++ b/docs/src/IceNucleationBox.md @@ -199,6 +199,19 @@ where: where: - ``N_{aer}`` is the number of available dust aerosol particles. +## 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}``, + per second via immersion freezing can then be calculating using +```math +\begin{equation} + P_{ice} = [\frac{dN_i}{dt}]_{immer} = J_{immer}A(N_{liq}) + \label{eq:ABIFM_P_ice} +\end{equation} +``` +where ``N_{liq}`` is total number of ice nuclei containing droplets and + ``A`` is surface area of those droplets. + ## Example figures Here we show example simulation results from the adiabatic parcel @@ -224,3 +237,13 @@ It is compared to [Rogers1975](@cite). include("../../parcel/Liquid_only.jl") ``` ![](liquid_only_parcel.svg) + +The plots below are the results of the adiabatic parcel model + with immersion freezing, condensation growth, and deposition growth for + both a monodisperse and gamma size distribution. Note that this has not + yet been validated against literature. + +```@example +include("../../parcel/Immersion_Freezing.jl") +``` +![](immersion_freezing.svg) diff --git a/parcel/Immersion_Freezing.jl b/parcel/Immersion_Freezing.jl new file mode 100644 index 0000000000..80c3850926 --- /dev/null +++ b/parcel/Immersion_Freezing.jl @@ -0,0 +1,126 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CloudMicrophysics.CommonTypes as CMT +import CLIMAParameters as CP + +# boilerplate code to get free parameter values +include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) +# definition of the ODE problem for parcel model +include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +# Boiler plate code to have access to model parameters and constants +FT = Float64 +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +prs = cloud_microphysics_parameters(toml_dict) +thermo_params = CMP.thermodynamics_params(prs) +air_props = CMT.AirProperties(FT) + +# Constants +ρₗ = FT(CMP.ρ_cloud_liq(prs)) +R_v = FT(CMP.R_v(prs)) +R_d = FT(CMP.R_d(prs)) + +# Initial conditions +Nₐ = FT(0) +Nₗ = FT(500 * 1e3) +Nᵢ = FT(0) +r₀ = FT(1e-6) +p₀ = FT(800 * 1e2) +T₀ = FT(251) +qᵥ = FT(8.1e-4) +qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ +qᵢ = FT(0) +x_sulph = FT(0.01) + +# Moisture dependent initial conditions +q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +ts = TD.PhaseNonEquil_pTq(thermo_params, p₀, T₀, q) +ρₐ = TD.air_density(thermo_params, ts) +Rₐ = TD.gas_constant_air(thermo_params, q) +eₛ = TD.saturation_vapor_pressure(thermo_params, 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(1) # model timestep +t_max = FT(1200) +aerosol_type = CMT.IlliteType() +ice_nucleation_modes = ["ImmersionFreezing"] +growth_modes = ["Condensation", "Deposition"] +droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]] + +# Plotting +fig = MK.Figure(resolution = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]") +ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]") +ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]") +ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq") +ax6 = MK.Axis( + fig[3, 2], + xlabel = "Time [min]", + ylabel = "N_ice / N_tot", + yscale = log10, +) +MK.ylims!(ax6, 1e-6, 1) + +for droplet_size_distribution in droplet_size_distribution_list + p = (; + prs, + air_props, + thermo_params, + const_dt, + r_nuc, + w, + α_m, + aerosol_type, + ice_nucleation_modes, + growth_modes, + droplet_size_distribution, + ) + # solve ODE + sol = run_parcel(IC, FT(0), t_max, p) + + DSD = droplet_size_distribution[1] + + # Plot results + ξ(T) = + TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) / + TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) + S_i(T, S_liq) = ξ(T) * S_liq - 1 + + MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD) + MK.lines!(ax2, sol.t / 60, sol[3, :]) + MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3) + MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3) + + sol_Nₗ = sol[8, :] + sol_Nᵢ = sol[9, :] + sol_qₗ = sol[5, :] + if DSD == "Monodisperse" + rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ) + end + if DSD == "Gamma" + λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) + rₗ = 2 ./ λ + end + MK.lines!(ax5, sol.t / 60, sol_Nₗ) + MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ)) +end + +MK.axislegend( + ax1, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :rb, +) + +MK.save("immersion_freezing.svg", fig) diff --git a/parcel/Liquid_only.jl b/parcel/Liquid_only.jl index c46a985e2b..669e746908 100644 --- a/parcel/Liquid_only.jl +++ b/parcel/Liquid_only.jl @@ -58,8 +58,10 @@ w = FT(10) # updraft speed α_m = FT(0.5) # accomodation coefficient const_dt = FT(0.5) # model timestep t_max = FT(20) +aerosol_type = [] ice_nucleation_modes = [] # no freezing -growth_modes_list = [["Condensation"], ["Condensation_DSD"]] # switch on condensation +growth_modes = ["Condensation"] # switch on condensation +droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]] # Data from Rogers(1975) Figure 1 # https://www.tandfonline.com/doi/abs/10.1080/00046973.1975.9648397 @@ -80,7 +82,7 @@ ax5 = MK.Axis(fig[1, 2], ylabel = "radius [μm]") MK.lines!(ax1, Rogers_time_supersat, Rogers_supersat, label = "Rogers_1975") MK.lines!(ax5, Rogers_time_radius, Rogers_radius) -for growth_modes in growth_modes_list +for droplet_size_distribution in droplet_size_distribution_list p = (; prs, air_props, @@ -89,26 +91,28 @@ for growth_modes in growth_modes_list r_nuc, w, α_m, + aerosol_type, ice_nucleation_modes, growth_modes, + droplet_size_distribution, ) # solve ODE sol = run_parcel(IC, FT(0), t_max, p) - gm = growth_modes[1] + DSD = droplet_size_distribution[1] # Plot results - MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = gm) + MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = DSD) MK.lines!(ax2, sol.t, sol[3, :]) MK.lines!(ax3, sol.t, sol[4, :] * 1e3) MK.lines!(ax4, sol.t, sol[5, :] * 1e3) sol_Nₗ = sol[8, :] sol_qₗ = sol[5, :] - if gm == "Condensation" + if DSD == "Monodisperse" rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ) end - if gm == "Condensation_DSD" + if DSD == "Gamma" λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) rₗ = 2 ./ λ end diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl index bec763be5f..23a4f5976b 100644 --- a/parcel/Tully_et_al_2023.jl +++ b/parcel/Tully_et_al_2023.jl @@ -78,8 +78,10 @@ function Tully_et_al_2023(FT) w = FT(3.5 * 1e-2) # updraft speed α_m = FT(0.5) # accomodation coefficient const_dt = 0.1 # model timestep + aerosol_type = CMT.DesertDustType() # aerosol type ice_nucleation_modes = ["DustDeposition"] # switch on deposition on dust growth_modes = ["Deposition"] # switch on deposition growth + droplet_size_distribution = ["Monodisperse"] p = (; prs, air_props, @@ -88,8 +90,10 @@ function Tully_et_al_2023(FT) r_nuc, w, α_m, + aerosol_type, ice_nucleation_modes, growth_modes, + droplet_size_distribution, ) # Simulation 1 diff --git a/parcel/parcel.jl b/parcel/parcel.jl index f7e78c22a3..c3f0c67882 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -15,8 +15,8 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) function parcel_model(dY, Y, p, t) # Get simulation parameters - (; prs, air_props, thermo_params, const_dt, r_nuc, w, α_m) = p - (; ice_nucleation_modes, growth_modes) = p + (; prs, air_props, thermo_params, const_dt, r_nuc, w, α_m, aerosol_type) = p + (; ice_nucleation_modes, growth_modes, droplet_size_distribution) = p # Numerical precision used in the simulation FT = eltype(Y) @@ -40,6 +40,7 @@ function parcel_model(dY, Y, p, t) grav = CMP.grav(prs) ρ_ice = CMP.ρ_cloud_ice(prs) ρ_liq = CMP.ρ_cloud_liq(prs) + H2SO4_prs = CMP.H2SO4SolutionParameters(FT) # Get thermodynamic parameters, phase partition and create thermo state. thermo_params = CMP.thermodynamics_params(prs) @@ -58,6 +59,7 @@ function parcel_model(dY, Y, p, t) L_subl = TD.latent_heat_sublim(thermo_params, T) L_vap = TD.latent_heat_vapor(thermo_params, T) ρ_air = TD.air_density(thermo_params, ts) + e = q_vap * p_a * R_v / R_a # Adiabatic parcel coefficients a1 = L_vap * grav / cp_a / T^2 / R_v - grav / R_a / T @@ -67,23 +69,48 @@ function parcel_model(dY, Y, p, t) # TODO - we should zero out all tendencies and augemnt them # TODO - add immersion, homogeneous, ... - #dN_act_dt_immersion = FT(0) #dN_act_dt_homogeneous = FT(0) dN_act_dt_depo = FT(0) + dqi_dt_new_depo = FT(0) if "DustDeposition" in ice_nucleation_modes AF = CMI_het.dust_activated_number_fraction( prs, S_i, T, - CMT.DesertDustType(), + aerosol_type, ) 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 + + dN_act_dt_immersion = FT(0) + dqi_dt_new_immers = FT(0) + if "ImmersionFreezing" in ice_nucleation_modes + Δa_w = T > FT(185) && T < FT(235) ? + CMO.a_w_xT(H2SO4_prs, thermo_params, x_sulph, T) - CMO.a_w_ice(thermo_params, T) : + CMO.a_w_eT(thermo_params, e, T) - CMO.a_w_ice(thermo_params, T) + J_immersion = CMI_het.ABIFM_J(prs, aerosol_type, Δa_w) + if "Monodisperse" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes + 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) + 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 + λ = 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) + dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + end end - dN_ice_dt = dN_act_dt_depo - dN_aer_dt = -dN_act_dt_depo - dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air + dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion + dN_aer_dt = -dN_act_dt_depo + dN_liq_dt = -dN_act_dt_immersion + # Growth dqi_dt_depo = FT(0) if "Deposition" in growth_modes && N_ice > 0 @@ -97,24 +124,24 @@ function parcel_model(dY, Y, p, t) dql_dt_cond = FT(0) if "Condensation" in growth_modes && N_liq > 0 - # Condensation on existing droplets (assuming all are the same) - G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid()) - r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) - dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq - end - if "Condensation_DSD" in growth_modes && N_liq > 0 - # Condensation on existing droplets - # assuming n(r) = A r exp(-λr) - G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid()) - λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) - #A = N_liq* λ^2 - r_l = 2 / λ - dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq + if "Monodisperse" in droplet_size_distribution + # Condensation on existing droplets assuming all are the same + G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid()) + r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) + dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq + elseif "Gamma" in droplet_size_distribution + # Condensation on existing droplets assuming n(r) = A r exp(-λr) + G_l = CMO.G_func(air_props, thermo_params, T, TD.Liquid()) + λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) + #A = N_liq* λ^2 + r_l = 2 / λ + dql_dt_cond = 4 * π / ρ_air * (S_liq - 1) * G_l * r_l * N_liq + end end # Update the tendecies - dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo - dq_liq_dt = dql_dt_cond + dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers + dq_liq_dt = dql_dt_cond - dqi_dt_new_immers dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt dp_a_dt = -p_a * grav / R_a / T * w dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt @@ -128,7 +155,7 @@ function parcel_model(dY, Y, p, t) dY[5] = dq_liq_dt # liquid water specific humidity dY[6] = dq_ice_dt # ice specific humidity dY[7] = dN_aer_dt # number concentration of interstitial aerosol - dY[8] = FT(0) # mumber concentration of droplets + dY[8] = dN_liq_dt # mumber concentration of droplets dY[9] = dN_ice_dt # number concentration of activated particles dY[10] = FT(0) # sulphuric acid concentration @@ -181,6 +208,9 @@ function run_parcel(IC, t_0, t_end, p) if "DustDeposition" in ice_nucleation_modes print("Deposition on dust particles ") end + if "ImmersionFreezing" in ice_nucleation_modes + print("Immersion freezing") + end print("\n") print("Growth modes: ") if "Condensation" in growth_modes