From 7e7cb43a66f96ab7c17564a2925749fbcb8c665e Mon Sep 17 00:00:00 2001 From: amylu00 Date: Mon, 23 Oct 2023 11:39:54 -0700 Subject: [PATCH] homogeneous freezing to parcel --- Project.toml | 4 +- parcel/Jensen_et_al_2022.jl | 158 ++++++++++++++++++++++++++++++++++++ parcel/Project.toml | 4 +- parcel/parcel.jl | 125 ++++++++++++++++++++-------- src/IceNucleation.jl | 11 ++- 5 files changed, 261 insertions(+), 41 deletions(-) create mode 100644 parcel/Jensen_et_al_2022.jl diff --git a/Project.toml b/Project.toml index 2e86de81da..b062f57f07 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CloudMicrophysics" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" authors = ["Climate Modeling Alliance"] -version = "0.15.1" +version = "0.15.2" [deps] CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" @@ -12,7 +12,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = "0.7.15" +CLIMAParameters = "0.8" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" RootSolvers = "0.3, 0.4" diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl new file mode 100644 index 0000000000..99f5097baa --- /dev/null +++ b/parcel/Jensen_et_al_2022.jl @@ -0,0 +1,158 @@ +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) + +# 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) # Jensen2022 starts with Nₐ = 300 * 1e6 and activates them within parcel +Nᵢ = FT(0) +r₀ = FT(25 * 1e-9) # Value of dry aerosol r from Jensen2022 +#r₀ = FT(7 * 1e-8) # starting when Jensen freezes +p₀ = FT(121 * 1e2) +T₀ = FT(190) +#T₀ = FT(189.7) # starting when Jensen freezes +x_sulph = FT(0) + +# saturation and partial pressure +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +# Sᵢ = FT(1.55) +# Sₗ = S_l(T₀, Sᵢ) +# e = eₛ * Sₗ +# ϵ = R_d / R_v +# starting when Jensen freezes +#q_vap = FT(3.76e-7) +q_vap = FT(5e-6) * FT(0.6) / FT(1.2) +#q_vap = FT(4.1e-7) # this gives more valid ICNC and ice saturation curves +q_liq = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 +q_ice = FT(0) +#q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice) +R_a = TD.gas_constant_air(tps, q) +e = q_vap * p₀ * R_v / R_a +sₗ = e / eₛ # saturation ratio +println("sᵢ is ", s_i(T₀, sₗ)) # saturation ratio + +# 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 + +#qᵥ = mv_v / (md_v + mv_v + ml_v) +#qₗ = ml_v / (md_v + mv_v + ml_v) +#qᵢ = 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) + +#IC = [sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] +IC = [sₗ, p₀, T₀, q_vap, q_liq, q_ice, Nₐ, Nₗ, Nᵢ, x_sulph] # starting when Jensen freezes + +println("initial conditions:") +println(" liquid saturation: ", sₗ) +# println(" vapor mixing ratio: ", qᵥ) +# println(" liquid mixing ratio: ", qₗ) +# println(" ice mixing ratio: ", qᵢ) +println(" vapor mixing ratio: ", q_vap) # starting when Jensen freezes +println(" liquid mixing ratio: ", q_liq) # starting when Jensen freezes + +# 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(140) +#t_max = FT(84) # start when Jensen freezes +aerosol = [] +ice_nucleation_modes = ["HomogeneousFreezing"] # homogeneous freezing only +growth_modes = ["Deposition"] +droplet_size_distribution_list = [["Lognormal"]] + +# 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 = "Ice 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 = "r_liq [m]") + +MK.ylims!(ax2, 188.5, 190) +MK.xlims!(ax2, -5, 150) +MK.xlims!(ax3, -5, 150) +MK.xlims!(ax4, -5, 150) + +MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen 2022", color = :green) +MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green) +MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green) + +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, + ) + # 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 = DSD, color = :blue) + MK.lines!(ax1, sol.t, (sol[1, :]), color = :red) # liq saturation + MK.lines!(ax2, sol.t, sol[3, :]) # 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)# ICNC +end + +#MK.save("Jensen_et_al_2022.svg", fig) +fig \ No newline at end of file diff --git a/parcel/Project.toml b/parcel/Project.toml index 5fbddd66e4..2a1dcf6247 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -4,11 +4,9 @@ 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" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" - -[compat] -CLIMAParameters = "0.7.12" diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 5b91f20a73..5e8b119174 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -3,7 +3,9 @@ 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 +import Random as RD #! format: off """ @@ -21,7 +23,7 @@ function parcel_model(dY, Y, p, t) # TODO - We are solving for both q_ and supersaturation. # We should decide if we want to solve for S (as in the cirrus box model) # or for q_ which is closer to what ClimaAtmos is doing. - S_liq = Y[1] # saturation ratio over liquid water + s_liq = Y[1] # saturation ratio over liquid water p_a = Y[2] # pressure T = Y[3] # temperature q_vap = Y[4] # vapor specific humidity @@ -32,6 +34,8 @@ function parcel_model(dY, Y, p, t) N_ice = Y[9] # number concentration of activated ice crystals x_sulph = Y[10] # percent mass sulphuric acid + + # Constants R_v = TD.Parameters.R_v(tps) grav = TD.Parameters.grav(tps) @@ -47,21 +51,24 @@ function parcel_model(dY, Y, p, t) e_si = TD.saturation_vapor_pressure(tps, T, TD.Ice()) e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) ξ = e_sl / e_si - S_i = ξ * S_liq + s_i = ξ * s_liq + println("s_i = ",s_i) # Constants and variables that depend on the moisture content R_a = TD.gas_constant_air(tps, q) cp_a = TD.cp_m(tps, q) L_subl = TD.latent_heat_sublim(tps, T) + L_fus = TD.latent_heat_fusion(tps, T) L_vap = TD.latent_heat_vapor(tps, T) ρ_air = TD.air_density(tps, 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 - a2 = 1 / q_vap + a2 = R_v * p_a / e_sl / R_a a3 = L_vap^2 / R_v / T^2 / cp_a a4 = L_vap * L_subl / R_v / T^2 / cp_a + 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, ... @@ -73,7 +80,7 @@ function parcel_model(dY, Y, p, t) AF = CMI_het.dust_activated_number_fraction( aerosol, ip.deposition, - S_i, + s_i, T, ) dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt @@ -86,65 +93,117 @@ function parcel_model(dY, Y, p, t) Δa_w = T > FT(185) && T < FT(235) ? 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 + J_immersion = CMI_het.ABIFM_J(aerosol, ip, Δa_w) + 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) - dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + elseif "Gamma" in droplet_size_distribution + λ = cbrt(32 * π * 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 + 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 + + dN_act_dt_hom = FT(0) + dqi_dt_new_hom = FT(0) + if "HomogeneousFreezing" in ice_nucleation_modes + a_w_ice_μ = exp((210368 + 131.438*T - (3.32373e6 /T) - 41729.1*log(T))/(8.31441*T)) + Δa_w = CMO.a_w_eT(tps, e, T) - a_w_ice_μ + println("Δa_w = ", Δa_w) + J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w) + if "Monodisperse" in droplet_size_distribution + #r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air) + V_aero = 4 / 3 * π * r_nuc^3 + dN_act_dt_hom = max(FT(0), J_hom * N_aer * V_aero) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air + elseif "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) - dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air + V_aero = 4 / 3 * π * r_nuc^3 + dN_act_dt_hom = max(FT(0), J_hom * N_aer * V_aero) + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air + elseif "Lognormal" in droplet_size_distribution + #dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(R_mode_aero) + 9 * σ^2)/(2))) + dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(r_nuc) + 36)/(2))) + r_aero = r_nuc * exp(2) # mean radius + dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_aero^3 * ρ_ice / ρ_air end + + println("J_hom = ", J_hom) + println("dqi_dt_new_hom = ", dqi_dt_new_hom) + println("dN_act_dt_hom = ", dN_act_dt_hom) 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 + dN_aer_dt = -dN_act_dt_depo - dN_act_dt_hom + dN_liq_dt = -dN_act_dt_immersion #- dN_act_dt_hom + println("N_ice = ", N_ice) # 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 * π) / ρ_ice * ρ_air) + C_i = r_i + dqi_dt_depo = 4 * π / ρ_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 * π * N_ice / q_ice * ρ_ice / ρ_air) + #A = N_ice* λ^2 + r_i = 2 / λ + C_i = r_i + dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice + elseif "Lognormal" in droplet_size_distribution # TODO - monodisperse for now. need to figure out how to do this lognormally + r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) + dqi_dt_depo = 4 * π / ρ_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 + # if r_l < FT(0) + # r_l = FT(0) + # end + println("r_l = ", r_l) 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) #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 * π / ρ_air * (s_liq - 1) * G_l * r_l * N_liq + println("dql_dt_cond = ", dql_dt_cond) end + dq_liq_dt_lv = dql_dt_cond # from liq-vap transitions + dq_liq_dt_li = -dqi_dt_new_immers #- dqi_dt_new_hom # from liq-ice transitions + dq_ice_dt_iv = dqi_dt_new_depo + dqi_dt_depo # from ice-vap transitions + dq_ice_dt_il = -dq_liq_dt_li + dqi_dt_new_hom # from ice-liq transitions + # Update the tendecies - 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 + dq_ice_dt = dq_ice_dt_iv + dq_ice_dt_il + dq_liq_dt = dq_liq_dt_lv + dq_liq_dt_li + ds_l_dt = a1 * w * s_liq - (a2 + a3 * s_liq) * dq_liq_dt_lv - (a2 + a4 * s_liq) * dq_ice_dt_iv - a5 * s_liq * dq_ice_dt_il 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 - dq_vap_dt = -dq_ice_dt - dq_liq_dt + #dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_lv + L_fus / cp_a * dq_ice_dt_il + L_subl / cp_a * dq_ice_dt_iv + dT_dt = FT(- 1 / 120) + println("------------ dT/dt = ", dT_dt) + println("------------ T = ", T) + println("------------ q_vap = ", q_vap) + dq_vap_dt = -dq_ice_dt_iv - dq_liq_dt_lv # Set tendencies - dY[1] = dS_l_dt # saturation ratio over liquid water + dY[1] = ds_l_dt # saturation ratio over liquid water dY[2] = dp_a_dt # pressure dY[3] = dT_dt # temperature dY[4] = dq_vap_dt # vapor specific humidity @@ -156,7 +215,7 @@ function parcel_model(dY, Y, p, t) dY[10] = FT(0) # sulphuric acid concentration # TODO - add diagnostics output (radius, S, etc) - + println("\n") end #! format: on @@ -167,13 +226,13 @@ Returns the solution of an ODE probelm defined by the parcel model. Inputs: - IC - A vector with the initial conditions for - [S_l, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph] + [s_l, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph] - t_0 - simulation start time - t_end - simulation end time - p - a named tuple with simulation parameters. Initial condition contains (all in base SI units): - - S_l - saturation ratio over liquid water, + - s_l - saturation ratio over liquid water, - p_a - atmospheric pressure - T - temperature - q_vap - water vapor specific humidity diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 2c80312a4b..3e50719bee 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -85,9 +85,14 @@ 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 + if Δa_w < ip.Δa_w_min + Δa_w = ip.Δa_w_min + elseif Δa_w > ip.Δa_w_max + Δa_w = ip.Δa_w_max + end + logJ::FT = ip.c₁ + ip.c₂ * Δa_w - ip.c₃ * Δa_w^2 + ip.c₄ * Δa_w^3 return FT(10)^(logJ) * 1e6