From a31a5cb9c857578d127e584242853947bc01589d Mon Sep 17 00:00:00 2001 From: amylu00 Date: Thu, 1 Feb 2024 16:30:47 -0800 Subject: [PATCH] replacing Mohler AF with nucleation rate --- docs/src/IceNucleation.md | 23 ++-- docs/src/IceNucleationParcel0D.md | 27 ++--- parcel/Tully_et_al_2023.jl | 188 +++++++++++++++++------------- parcel/parcel.jl | 20 +++- 4 files changed, 144 insertions(+), 114 deletions(-) diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index f1987dacd1..0d1059e1f9 100644 --- a/docs/src/IceNucleation.md +++ b/docs/src/IceNucleation.md @@ -22,28 +22,25 @@ The parameterization for deposition on dust particles is an implementation of freezing modes. ## Activated fraction for deposition freezing on dust -The parameterization models the activated fraction +The parameterization models the nucleation rate of ice as an empirical function of ice saturation ratio, - see eq. (3) in [Mohler2006](@cite). + see eq. (5) in [Mohler2006](@cite). ```math \begin{equation} -f_i(S_i) = exp[a(S_i - S_0)] - 1 +\frac{dn_{ice}}{dt} = N_{aer} a \frac{dS_i}{dt} \end{equation} ``` where: - - ``f_i`` is the activated fraction - (the ratio of aerosol particles acting as ice nuclei to the total number of aerosol particles), + - ``N_{aer}`` is the number of available aerosol/ice nuclei, + - ``a`` is a scaling parameter dependent on aerosol properties and temperature, - ``S_i`` is the ice saturation ratio - (the ratio of water vapor partial pressure and the water vapor partial pressure at saturation over ice), - - ``S_0`` is the threshold ice saturation ratio, - - ``a`` is a scaling parameter dependent on aerosol properties and temperature. + (the ratio of water vapor partial pressure and the water vapor partial pressure at saturation over ice). -Limited experimental values for the free parameters ``S_0`` and ``a`` can be found in [Mohler2006](@cite). -Both parameters are dependent on aerosol properties and temperature. +Limited experimental values for the free parameter ``a`` can be found in [Mohler2006](@cite). Parameter ``a`` is strongly dependent on aerosol properties and temperature. !!! note - For a ``f_i`` values above 0.08 or ``S_i`` between 1.35 and 1.5, + For ``S_i`` between 1.35 and 1.5, freezing occurs in a different ice nucleation mode (either a second deposition or other condensation type mode). @@ -97,10 +94,10 @@ Once ``J_{ABIFM}`` is calculated, it can be used to determine the ice production per second via immersion freezing. ```math \begin{equation} - P_{ice} = J_{ABIFM}A(N_{tot} - N_{ice}) + P_{ice} = J_{ABIFM}A(N_{aer} - N_{ice}) \end{equation} ``` -where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number +where ``A`` is surface area of an individual ice nuclei, ``N_{aer}`` is total number of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. ### ABIFM Example Figures diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 1117889b7a..ed726b6606 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -121,7 +121,7 @@ where ``\xi = \frac{e_{sl}}{e_{si}}`` is the ratio of saturation vapor pressure The crux of the problem is modeling the ``\frac{dq_l}{dt}`` and ``\frac{dq_i}{dt}`` for different homogeneous and heterogeneous ice nucleation paths. -## Condensation +## Condensation growth The diffusional growth of individual cloud droplet is described by ([see discussion](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics1M/#Snow-autoconversion)), @@ -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 nucleation on dust particles +## Deposition growth Similarly, for a case of a spherical ice particle growing through water vapor deposition ```math @@ -196,15 +196,12 @@ It follows that 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/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i`` -```math -\begin{equation} - N_{act} = N_{aer} f_i -\end{equation} -``` -where: -- ``N_{aer}`` is the number of available dust aerosol particles. +## Deposition Nucleation on dust particles +There are multiple ways of running deposition nucleation in the parcel. + ``"MohlerAF_Deposition"`` will trigger an activated fraction approach from [Mohler2006](@cite). + ``"MohlerRate_Deposition"`` will trigger a nucleation rate approach from [Mohler2006](@cite). +The deposition nucleation methods are parameterized as described in + [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/). ## Immersion Freezing Following the water activity based immersion freezing model (ABIFM), the ABIFM derived @@ -241,13 +238,17 @@ Here we show various example simulation results from the adiabatic parcel 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). +The model is run three times using the ``"MohlerAF_Deposition"`` approach + for 30 minutes simulation time, (shown by three different colors on the plot). Between each run the water vapor specific humidity is changed, while keeping all other state variables the same as at the last time step of the previous run. The prescribed vertical velocity is equal to 3.5 cm/s. Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines). +The pale blue line uses the ``"MohlerRate_Deposition"``approach. When we force + the initial ice crystal number concentration for this simulation to match + that in the ``"MohlerAF_Deposition"`` approach, we obtain the same results as + in the `"MohlerAF_Deposition"` approach. ```@example include("../../parcel/Tully_et_al_2023.jl") diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl index 95fe557f8c..418a306e4c 100644 --- a/parcel/Tully_et_al_2023.jl +++ b/parcel/Tully_et_al_2023.jl @@ -74,76 +74,11 @@ function Tully_et_al_2023(FT) α_m = FT(0.5) # accomodation coefficient const_dt = 0.1 # model timestep aerosol = CMP.DesertDust(FT) # aerosol type - ice_nucleation_modes = ["DustDeposition"] # switch on deposition on dust + ice_nucleation_modes_list = [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]] # ice nucleation modes growth_modes = ["Deposition"] # switch on deposition growth droplet_size_distribution = ["Monodisperse"] - p = (; - wps, - aps, - tps, - ip, - const_dt, - r_nuc, - w, - α_m, - aerosol, - ice_nucleation_modes, - growth_modes, - droplet_size_distribution, - ) - - # Simulation 1 - IC1 = get_initial_condition( - tps, - p_0, - T_0, - q_vap_0, - q_liq_0, - q_ice_0, - N_aerosol, - N_droplets, - N_0, - x_sulph, - ) - sol1 = run_parcel(IC1, 0, t_max, p) - - # Simulation 2 - # (alternatively set T and take q_vap from the previous simulation) - #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end]) - IC2 = get_initial_condition( - tps, - sol1[2, end], - #sol1[3, end], - T2, - q_vap2, - q_liq_0, - sol1[6, end], - sol1[7, end], - sol1[8, end], - sol1[9, end], - x_sulph, - ) - sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p) - - # Simulation 3 - # (alternatively set T and take q_vap from the previous simulation) - #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end]) - IC3 = get_initial_condition( - tps, - sol2[2, end], - #sol2[3, end], - T3, - q_vap3, - q_liq_0, - sol2[6, end], - sol2[7, end], - sol2[8, end], - sol2[9, end], - x_sulph, - ) - sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p) - - # Plot results + + # Plots fig = MK.Figure(resolution = (800, 600)) ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]") ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") @@ -160,23 +95,108 @@ function Tully_et_al_2023(FT) TD.saturation_vapor_pressure(tps, T, TD.Ice()) S_i(T, S_liq) = ξ(T) * S_liq - 1 - sol = [sol1, sol2, sol3] - clr = ["blue", "orange", "green"] - for it in [1, 2, 3] - MK.lines!(ax1, sol[it].t * w, sol[it][1, :] .- 1, color = clr[it]) - MK.lines!(ax2, sol[it].t * w, sol[it][3, :], color = clr[it]) - MK.lines!(ax3, sol[it].t * w, sol[it][9, :] * 1e-3, color = clr[it]) - MK.lines!(ax4, sol[it].t * w, sol[it][7, :] * 1e-3, color = clr[it]) - MK.lines!(ax5, sol[it].t * w, sol[it][4, :] * 1e3, color = clr[it]) - MK.lines!(ax6, sol[it].t * w, sol[it][6, :] * 1e3, color = clr[it]) - - MK.lines!( - ax1, - sol[it].t * w, - S_i.(sol[it][3, :], sol[it][1, :]), - linestyle = :dash, - color = clr[it], + for ice_nucleation_modes in ice_nucleation_modes_list + p = (; + wps, + aps, + tps, + ip, + const_dt, + r_nuc, + w, + α_m, + aerosol, + ice_nucleation_modes, + growth_modes, + droplet_size_distribution, + ) + + # Simulation 1 + IC1 = get_initial_condition( + tps, + p_0, + T_0, + q_vap_0, + q_liq_0, + q_ice_0, + N_aerosol, + N_droplets, + N_0, + x_sulph, ) + sol1 = run_parcel(IC1, 0, t_max, p) + if ice_nucleation_modes == ["MohlerAF_Deposition"] + # Simulation 2 + # (alternatively set T and take q_vap from the previous simulation) + #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end]) + IC2 = get_initial_condition( + tps, + sol1[2, end], + #sol1[3, end], + T2, + q_vap2, + q_liq_0, + sol1[6, end], + sol1[7, end], + sol1[8, end], + sol1[9, end], + x_sulph, + ) + sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p) + + # Simulation 3 + # (alternatively set T and take q_vap from the previous simulation) + #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end]) + IC3 = get_initial_condition( + tps, + sol2[2, end], + #sol2[3, end], + T3, + q_vap3, + q_liq_0, + sol2[6, end], + sol2[7, end], + sol2[8, end], + sol2[9, end], + x_sulph, + ) + sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p) + + # Plot results + sol = [sol1, sol2, sol3] + clr = ["blue", "orange", "green"] + for it in [1, 2, 3] + MK.lines!(ax1, sol[it].t * w, sol[it][1, :] .- 1, color = clr[it]) + MK.lines!(ax2, sol[it].t * w, sol[it][3, :], color = clr[it]) + MK.lines!(ax3, sol[it].t * w, sol[it][9, :] * 1e-3, color = clr[it]) + MK.lines!(ax4, sol[it].t * w, sol[it][7, :] * 1e-3, color = clr[it]) + MK.lines!(ax5, sol[it].t * w, sol[it][4, :] * 1e3, color = clr[it]) + MK.lines!(ax6, sol[it].t * w, sol[it][6, :] * 1e3, color = clr[it]) + MK.lines!( + ax1, + sol[it].t * w, + S_i.(sol[it][3, :], sol[it][1, :]), + linestyle = :dash, + color = clr[it], + ) + end + + elseif ice_nucleation_modes == ["MohlerRate_Deposition"] + MK.lines!(ax1, sol1.t * w, sol1[1, :] .- 1, color = :lightblue) + MK.lines!(ax2, sol1.t * w, sol1[3, :], color = :lightblue) + MK.lines!(ax3, sol1.t * w, sol1[9, :] * 1e-3, color = :lightblue) + MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3, color = :lightblue) + MK.lines!(ax5, sol1.t * w, sol1[4, :] * 1e3, color = :lightblue) + MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3, color = :lightblue) + + MK.lines!( + ax1, + sol1.t * w, + S_i.(sol1[3, :], sol1[1, :]), + linestyle = :dash, + color = :lightblue, + ) + end end MK.save("cirrus_box.svg", fig) end diff --git a/parcel/parcel.jl b/parcel/parcel.jl index 8ed506fa24..1f30d5e20a 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -75,7 +75,7 @@ function parcel_model(dY, Y, p, t) dN_act_dt_depo = FT(0) dqi_dt_new_depo = FT(0) - if "DustDeposition" in ice_nucleation_modes + if "MohlerAF_Deposition" in ice_nucleation_modes AF = CMI_het.dust_activated_number_fraction( aerosol, ip.deposition, @@ -85,6 +85,11 @@ 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 "MohlerRate_Deposition" in ice_nucleation_modes + a = T > ip.deposition.T_thr ? aerosol.a_warm : aerosol.a_cold + dN_act_dt_depo = max(FT(0), N_aer * a * (dY[1] * ξ)) + dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air + end dN_act_dt_immersion = FT(0) dqi_dt_new_immers = FT(0) @@ -268,14 +273,21 @@ function run_parcel(IC, t_0, t_end, p) println(" ") println("Running parcel model with: ") print("Ice nucleation modes: ") - if "DustDeposition" in ice_nucleation_modes - print("Deposition on dust particles ") + if "MohlerAF_Deposition" in ice_nucleation_modes + print("Deposition on dust particles using AF ") + timestepper = ODE.Euler() + end + if "MohlerRate_Deposition" in ice_nucleation_modes + print("Deposition nucleation based on rate eqn in Mohler 2006 ") + timestepper = ODE.Tsit5() end if "ImmersionFreezing" in ice_nucleation_modes print("Immersion freezing ") + timestepper = ODE.Tsit5() end if "HomogeneousFreezing" in ice_nucleation_modes print("Homogeneous freezing ") + timestepper = ODE.Tsit5() end print("\n") print("Growth modes: ") @@ -306,7 +318,7 @@ function run_parcel(IC, t_0, t_end, p) problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p) sol = ODE.solve( problem, - ODE.Euler(), + timestepper, dt = const_dt, reltol = 10 * eps(FT), abstol = 10 * eps(FT),