diff --git a/docs/bibliography.bib b/docs/bibliography.bib index bf5a6b65d..22b52ccee 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -182,6 +182,20 @@ @article{Heymsfield2003 year = {2003} } +@article{Jensen2022, +author = {Jensen, E. J. and Diskin, G. S. and DiGangi, J. and Woods, S. and Lawson, R. P. and Bui, T. V.}, +title = {Homogeneous Freezing Events Sampled in the Tropical Tropopause Layer}, +journal = {Journal of Geophysical Research: Atmospheres}, +volume = {127}, +number = {17}, +pages = {e2022JD036535}, +keywords = {cirrus, nucleation, freezing}, +doi = {https://doi.org/10.1029/2022JD036535}, +url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022JD036535}, +note = {e2022JD036535 2022JD036535}, +year = {2022} +} + @article{Karcher2006, author = {Kärcher, B. and Hendricks, J. and Lohmann, U.}, title = {Physically based parameterization of cirrus cloud formation for use in global atmospheric models}, diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 8aaf7c9ed..f5f0abe6f 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -165,7 +165,7 @@ For a gamma distribution of droplets ``n(r) = A \; r \; exp(-\lambda r)``, ``\bar{r} = \frac{2}{\lambda}`` where ``\lambda = \left(\frac{32 \pi N_{tot} \rho_l}{q_l \rho_a}\right)^{1/3}``. -## Deposition on dust particles +## Deposition nucleation on dust particles Similarly, for a case of a spherical ice particle growing through water vapor deposition ```math @@ -197,7 +197,7 @@ where: - ``N_{act}`` is the number of activated ice particles. ``N_{act}`` can be computed for example from - [activated fraction](https://clima.github.io/CloudMicrophysics.jl/previews/PR103/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i`` + [activated fraction](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i`` ```math \begin{equation} N_{act} = N_{aer} f_i @@ -212,17 +212,35 @@ Following the water activity based immersion freezing model (ABIFM), the ABIFM d 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}) + P_{ice} = \left[ \frac{dN_i}{dt} \right]_{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. +## Homogeneous Freezing +Homogeneous freezing follows the water-activity based model described in the + [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/) section which gives a nucleation rate coefficient of + units ``cm^{-3}s^{-1}``. +The ice production rate from homogeneous freezing can then be determined: +```math +\begin{equation} + P_{ice} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq}) + \label{eq:hom_P_ice} +\end{equation} +``` +where ``N_{liq}`` is total number of ice nuclei containing droplets and + ``V`` is the volume of those droplets. + ## Example figures -Here we show example simulation results from the adiabatic parcel - model with deposition freezing on dust. +Here we show various example simulation results from the adiabatic parcel + model. This includes examples with deposition nucleation on dust, + liquid processes only, immersion freezing with condensation and deposition growth, + and homogeneous freezing with deposition growth. + +We start with deposition freezing on dust. The model is run three times for 30 minutes simulation time, (shown by three different colors on the plot). Between each run the water vapor specific humidity is changed, @@ -254,3 +272,16 @@ The plots below are the results of the adiabatic parcel model include("../../parcel/Immersion_Freezing.jl") ``` ![](immersion_freezing.svg) + +The following plots show the parcel model running with homogeneous freezing and + depositional growth assuming a lognormal distribution of aerosols. + It is compared against [Jensen2022](@cite). Note that running with the initial + conditions described in [Jensen2022](@cite) results in a ``\Delta a_w`` smaller + than the minimum valid value for the ``J_{hom}`` parameterization. We have forced + the ``\Delta a_w`` to be equal to the minimum valid value in this example only + for demonstrative purposes. + +```@example +include("../../parcel/Jensen_et_al_2022.jl") +``` +![](Jensen_et_al_2022.svg) diff --git a/parcel/Immersion_Freezing.jl b/parcel/Immersion_Freezing.jl index 12a413f9b..3c2630a5e 100644 --- a/parcel/Immersion_Freezing.jl +++ b/parcel/Immersion_Freezing.jl @@ -26,7 +26,7 @@ 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ₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ qᵢ = FT(0) x_sulph = FT(0.01) @@ -102,10 +102,10 @@ for droplet_size_distribution in droplet_size_distribution_list sol_Nᵢ = sol[9, :] sol_qₗ = sol[5, :] if DSD == "Monodisperse" - rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ) + rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) end if DSD == "Gamma" - λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) + λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) rₗ = 2 ./ λ end MK.lines!(ax5, sol.t / 60, sol_Nₗ) diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl new file mode 100644 index 000000000..e6c728dca --- /dev/null +++ b/parcel/Jensen_et_al_2022.jl @@ -0,0 +1,159 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP +import Distributions as DS +import CloudMicrophysics.Parameters as CMP + +# 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) +wps = CMP.WaterProperties(FT) +aps = CMP.AirProperties(FT) +ip = CMP.IceNucleationParameters(FT) +H2SO4_prs = CMP.H2SO4SolutionParameters(FT) + +# Constants +ρₗ = wps.ρw +R_v = TD.Parameters.R_v(tps) +R_d = TD.Parameters.R_d(tps) +ϵₘ = R_d / R_v + +# Saturation conversion equations +ξ(T) = + TD.saturation_vapor_pressure(tps, T, TD.Liquid()) / + TD.saturation_vapor_pressure(tps, T, TD.Ice()) +S_i(T, S_liq) = @. S_liq * ξ(T) # saturation ratio + +# Initial conditions +Nₐ = FT(300 * 1e6) +Nₗ = FT(0) +Nᵢ = FT(0) +r₀ = FT(25 * 1e-9) +T₀ = FT(190) +cᵥ₀ = FT(5 * 1e-6) +x_sulph = FT(0) + +# saturation and partial pressure +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀) +qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 +qᵢ = FT(0) +q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +R_a = TD.gas_constant_air(tps, q) +Sᵢ = FT(1.55) +Sₗ = Sᵢ / ξ(T₀) +e = Sₗ * eₛ +p₀ = e / cᵥ₀ + +## Moisture dependent initial conditions +ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) + +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + +# Simulation parameters passed into ODE solver +r_nuc = r₀ # assumed size of nucleated particles +w = FT(1) # updraft speed +α_m = FT(1) # accomodation coefficient +const_dt = FT(0.01) # model timestep +t_max = FT(120) # total time +aerosol = [] # aerosol type. fill in if immersion or deposition freezing +R_mode_liq = r_nuc # lognormal distribution mode radius +σ = FT(2) # lognormal distribution std deviation +ice_nucleation_modes = ["HomogeneousFreezing"] # homogeneous freezing only +growth_modes = ["Deposition"] +droplet_size_distribution_list = [["Jensen_Example"]] + +# Data from Jensen(2022) Figure 1 +# https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JD036535 +#! format: off +Jensen_t_sat = [0, 62.71, 70.52, 76.87, 82.4, 84.84, 88.1, 92, 96.07, 100.63, 105.35, 112.51, 119.83] +Jensen_sat = [1.55, 1.694, 1.7107, 1.7208, 1.725, 1.726, 1.7259, 1.722, 1.715, 1.702, 1.686, 1.653, 1.6126] +Jensen_t_T = [0, 120] +Jensen_T = [190, 189] +Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.658, 94.123, 95.5877, 119.84] +Jensen_ICNC = [0, 0, 0.282, 0.789, 1.804, 4.1165, 7.218, 12.12, 16.35, 16.8, 16.97, 17.086] +#! format: on + +fig = MK.Figure(resolution = (1000, 1000)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Saturation") +ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]") +ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]") +ax5 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]") +ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]") + +MK.ylims!(ax2, 188.5, 190) +MK.xlims!(ax2, -5, 150) +MK.xlims!(ax3, -5, 150) +MK.xlims!(ax4, -5, 150) + +#! format: off +MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen et al 2022", color = :green) +MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green, label = "Jensen et al 2022") +MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green, label = "Jensen et al 2022") + +for droplet_size_distribution in droplet_size_distribution_list + p = (; + wps, + aps, + tps, + ip, + const_dt, + r_nuc, + w, + α_m, + aerosol, + ice_nucleation_modes, + growth_modes, + droplet_size_distribution, + R_mode_liq, + σ, + ) + # solve ODE + sol = run_parcel(IC, FT(0), t_max, p) + + DSD = droplet_size_distribution[1] + + # Plot results + MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = "ice", color = :blue) + MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation + MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl") # temperature + MK.lines!(ax3, sol.t, sol[4, :] * 1e3) # q_vap + MK.lines!(ax4, sol.t, sol[5, :] * 1e3) # q_liq + MK.lines!(ax5, sol.t, sol[9, :] * 1e-6, label = "CM.jl") # ICNC + MK.lines!(ax6, sol.t, sol[6, :] * 1e3) # q_ice + + MK.axislegend( + ax1, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, + ) + MK.axislegend( + ax5, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, + ) + MK.axislegend( + ax2, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, + ) +#! format: on +end + +MK.save("Jensen_et_al_2022.svg", fig) diff --git a/parcel/Liquid_only.jl b/parcel/Liquid_only.jl index ebd45ad76..7f6480d42 100644 --- a/parcel/Liquid_only.jl +++ b/parcel/Liquid_only.jl @@ -36,7 +36,7 @@ e = eₛ # mass per volume for dry air, vapor and liquid md_v = (p₀ - e) / R_d / T₀ mv_v = e / R_v / T₀ -ml_v = Nₗ * 4 / 3 * π * ρₗ * r₀^3 +ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀^3 qᵥ = mv_v / (md_v + mv_v + ml_v) qₗ = ml_v / (md_v + mv_v + ml_v) @@ -110,10 +110,10 @@ for droplet_size_distribution in droplet_size_distribution_list sol_Nₗ = sol[8, :] sol_qₗ = sol[5, :] if DSD == "Monodisperse" - rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ) + rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) end if DSD == "Gamma" - λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) + λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ) rₗ = 2 ./ λ end MK.lines!(ax5, sol.t, rₗ * 1e6) diff --git a/parcel/Project.toml b/parcel/Project.toml index f05e751c6..69532b37a 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -4,6 +4,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 8fc387460..8ed506fa2 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -3,6 +3,7 @@ import CloudMicrophysics as CM import CloudMicrophysics.Common as CMO import CloudMicrophysics.HetIceNucleation as CMI_het +import CloudMicrophysics.HomIceNucleation as CMI_hom import CloudMicrophysics.Parameters as CMP #! format: off @@ -14,6 +15,11 @@ function parcel_model(dY, Y, p, t) # Get simulation parameters (; aps, wps, tps, ip, const_dt, r_nuc, w, α_m, aerosol) = p (; ice_nucleation_modes, growth_modes, droplet_size_distribution) = p + if "Jensen_Example" in droplet_size_distribution + (; σ) = p + elseif "Lognormal" in droplet_size_distribution + (; R_mode_liq, σ) = p + end # Numerical precision used in the simulation FT = eltype(Y) @@ -30,7 +36,7 @@ function parcel_model(dY, Y, p, t) N_aer = Y[7] # number concentration of interstitial aerosol N_liq = Y[8] # number concentration of existing water droplets N_ice = Y[9] # number concentration of activated ice crystals - x_sulph = Y[10] # percent mass sulphuric acid + x_sulph = Y[10] # percent mass sulphuric acid # Constants R_v = TD.Parameters.R_v(tps) @@ -66,8 +72,6 @@ function parcel_model(dY, Y, p, t) a5 = L_vap * L_fus / R_v / cp_a / (T^2) # TODO - we should zero out all tendencies and augemnt them - # TODO - add immersion, homogeneous, ... - #dN_act_dt_homogeneous = FT(0) dN_act_dt_depo = FT(0) dqi_dt_new_depo = FT(0) @@ -79,7 +83,7 @@ function parcel_model(dY, Y, p, t) 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 + dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air end dN_act_dt_immersion = FT(0) @@ -89,65 +93,118 @@ 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 - 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 + if "Monodisperse" in droplet_size_distribution + r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) + elseif "Gamma" in droplet_size_distribution + λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air) + #A = N_liq* λ^2 + r_l = 2 / λ end - if "Gamma" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes - λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) + A_aer = 4 * FT(π) * 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 * FT(π) * r_l^3 * ρ_ice / ρ_air + end + + dN_act_dt_hom = FT(0) + dqi_dt_new_hom = FT(0) + if "HomogeneousFreezing" in ice_nucleation_modes + a_w_ice = CMO.a_w_ice(tps, T) + Δa_w = CMO.a_w_eT(tps, e, T) - a_w_ice + if "Monodisperse" in droplet_size_distribution + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) + V_l = 4 /3 * FT(π) * r_l^3 + dN_act_dt_hom = max(FT(0), J_hom * N_liq * V_l) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air + elseif "Gamma" in droplet_size_distribution + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + λ = cbrt(32 * FT(π) * 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 + V_l = 4 / 3 * FT(π) * r_l^3 + dN_act_dt_hom = max(FT(0), J_hom * N_aer * V_l) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air + elseif "Lognormal" in droplet_size_distribution + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + dN_act_dt_hom = max(FT(0), J_hom * N_liq * 4 / 3 * FT(π) * exp((6*log(R_mode_liq) + 9 * σ^2)/(2))) + r_l = R_mode_liq * exp(σ) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air + elseif "Jensen_Example" in droplet_size_distribution + Δa_w = ifelse(Δa_w < ip.homogeneous.Δa_w_min, ip.homogeneous.Δa_w_min, ifelse(Δa_w > ip.homogeneous.Δa_w_max, ip.homogeneous.Δa_w_max, Δa_w)) + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * FT(π) * exp((6*log(r_nuc) + 9 * σ^2)/(2))) + r_aero = r_nuc * exp(σ) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_aero^3 * ρ_ice / ρ_air end end - 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 + dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion + dN_act_dt_hom + if "Jensen_Example" in droplet_size_distribution + dN_aer_dt = -dN_act_dt_depo - dN_act_dt_hom + dN_liq_dt = -dN_act_dt_immersion + else + dN_aer_dt = -dN_act_dt_depo + dN_liq_dt = -dN_act_dt_immersion - dN_act_dt_hom + end # Growth dqi_dt_depo = FT(0) if "Deposition" in growth_modes && N_ice > 0 - # Deposition on existing crystals (assuming all are the same...) G_i = CMO.G_func(aps, tps, T, TD.Ice()) - r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) - C_i = r_i - dqi_dt_depo = 4 * π / ρ_air * (S_i - 1) * G_i * r_i * N_ice + if "Monodisperse" in droplet_size_distribution + # Deposition on existing crystals (assuming all are the same...) + r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) + C_i = r_i + dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice + elseif "Gamma" in droplet_size_distribution + # Condensation on existing droplets assuming n(r) = A r exp(-λr) + λ = cbrt(32 * FT(π) * N_ice / q_ice * ρ_ice / ρ_air) + #A = N_ice* λ^2 + r_i = 2 / λ + C_i = r_i + dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice + elseif "Lognormal" in droplet_size_distribution || "Jensen_Example" in droplet_size_distribution + r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) * exp(σ) + dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice + end # TODO - check if r_i is non-zero if initial conditions say q_ice = 0 end dql_dt_cond = FT(0) if "Condensation" in growth_modes && N_liq > 0 + G_l = CMO.G_func(aps, tps, T, TD.Liquid()) if "Monodisperse" in droplet_size_distribution # Condensation on existing droplets assuming all are the same - G_l = CMO.G_func(aps, tps, 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 + r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air) elseif "Gamma" in droplet_size_distribution # Condensation on existing droplets assuming n(r) = A r exp(-λr) - G_l = CMO.G_func(aps, tps, T, TD.Liquid()) - λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air) + λ = cbrt(32 * FT(π) * 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 + dql_dt_cond = 4 * FT(π) / ρ_air * (S_liq - 1) * G_l * r_l * N_liq end - dq_liq_dt_vap_to_liq = dql_dt_cond # from liq-vap transitions - dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo # from ice-vap transitions - dq_ice_dt_liq_to_ice = dqi_dt_new_immers # from ice-liq transitions + dq_liq_dt_vap_to_liq = dql_dt_cond # change in q_liq from vap-liq transitions + if "Jensen_Example" in droplet_size_distribution + # change in q_ice from liq-ice transitions + dq_ice_dt_liq_to_ice = dqi_dt_new_immers + # change in q_ice from vap-ice transitions + dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_hom + else + # change in q_ice from liq-ice transitions + dq_ice_dt_liq_to_ice = dqi_dt_new_immers + dqi_dt_new_hom + # change in q_ice from vap-ice transitions + dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo + end # Update the tendecies dq_ice_dt = dq_ice_dt_vap_to_ice + dq_ice_dt_liq_to_ice dq_liq_dt = dq_liq_dt_vap_to_liq - dq_ice_dt_liq_to_ice + dq_vap_dt = -dq_ice_dt - dq_liq_dt dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * S_liq * dq_ice_dt_vap_to_ice - a5 * S_liq * dq_ice_dt_liq_to_ice dp_a_dt = -p_a * grav / R_a / T * w dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_vap_to_liq + L_fus / cp_a * dq_ice_dt_liq_to_ice + L_subl / cp_a * dq_ice_dt_vap_to_ice - dq_vap_dt = -dq_ice_dt - dq_liq_dt # Set tendencies dY[1] = dS_l_dt # saturation ratio over liquid water @@ -162,7 +219,6 @@ function parcel_model(dY, Y, p, t) dY[10] = FT(0) # sulphuric acid concentration # TODO - add diagnostics output (radius, S, etc) - end #! format: on @@ -205,9 +261,9 @@ The named tuple p should contain: - droplet_size_distribution - a vector with assumed droplet size distribution. Possible options: ("Monodisperse", "Gamma") """ function run_parcel(IC, t_0, t_end, p) - + #! format: off FT = eltype(IC) - (; const_dt, ice_nucleation_modes, growth_modes) = p + (; const_dt, ice_nucleation_modes, growth_modes, droplet_size_distribution) = p println(" ") println("Running parcel model with: ") @@ -216,24 +272,37 @@ function run_parcel(IC, t_0, t_end, p) print("Deposition on dust particles ") end if "ImmersionFreezing" in ice_nucleation_modes - print("Immersion freezing") + print("Immersion freezing ") + end + if "HomogeneousFreezing" in ice_nucleation_modes + print("Homogeneous freezing ") end print("\n") print("Growth modes: ") if "Condensation" in growth_modes - print( - "Condensation on liquid droplets with monodisperse size distribution", - ) + print("Condensation on liquid droplets",) end if "Condensation_DSD" in growth_modes - print("Condensation on liquid droplets with gamma distribution") + print("Condensation on liquid droplets") end if "Deposition" in growth_modes - print("Deposition on ice crystals with monodisperse size distribution") + print("Deposition on ice crystals") end print("\n") + print("Size distribution modes: ") + if "Monodisperse" in droplet_size_distribution + print("Monodisperse size distribution") + end + if "Gamma" in droplet_size_distribution + print("Gamma size distribution") + end + if "Lognormal" in droplet_size_distribution + print("Lognormal size distribution") + end println(" ") + #! format: on + problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) sol = ODE.solve( problem, diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index f8ea26268..2c6270d47 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -107,8 +107,8 @@ Parameterization based on Koop 2000, DOI: 10.1038/35020537. """ function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT} - @assert Δa_w > ip.Δa_w_min - @assert Δa_w < ip.Δa_w_max + @assert Δa_w >= ip.Δa_w_min + @assert Δa_w <= ip.Δa_w_max logJ::FT = ip.c₁ + ip.c₂ * Δa_w - ip.c₃ * Δa_w^2 + ip.c₄ * Δa_w^3 diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl index 773cf5043..0ef8b14e5 100644 --- a/test/homogeneous_ice_nucleation_tests.jl +++ b/test/homogeneous_ice_nucleation_tests.jl @@ -36,11 +36,11 @@ function test_homogeneous_J(FT) ) # If Δa_w out of range - TT.@test_throws AssertionError("Δa_w > ip.Δa_w_min") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J( ip.homogeneous, Δa_w_too_small, ) - TT.@test_throws AssertionError("Δa_w < ip.Δa_w_max") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J( ip.homogeneous, Δa_w_too_large, )