diff --git a/docs/src/API.md b/docs/src/API.md index 9e9109fda..0a7fcf584 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -87,6 +87,7 @@ AerosolActivation.total_M_activated ```@docs HetIceNucleation HetIceNucleation.dust_activated_number_fraction +HetIceNucleation.MohlerDepositionRate HetIceNucleation.deposition_J HetIceNucleation.ABIFM_J ``` diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index f1987dacd..ea3c24ae9 100644 --- a/docs/src/IceNucleation.md +++ b/docs/src/IceNucleation.md @@ -22,28 +22,41 @@ 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 - as an empirical function of ice saturation ratio, - see eq. (3) in [Mohler2006](@cite). +There are 2 parameterizations from [Mohler2006](@cite) available: one + which calculates the activated fraction and one which calculates nucleation + rate. +The activated fraction parameterization follows eq. (3) in the paper. ```math \begin{equation} f_i(S_i) = exp[a(S_i - S_0)] - 1 \end{equation} ``` -where: +where - ``f_i`` is the activated fraction - (the ratio of aerosol particles acting as ice nuclei to the total number of aerosol particles), + (the ratio of aerosol particles acting as ice nuclei to the total number of aerosol particles), + - ``a`` is a scaling parameter dependent on aerosol properties and temperature, + - ``S_i`` is the ice saturation ratio, + - ``S_0`` is an empirically derived threshold ice saturation ratio. +The other parameterization models the nucleation rate of ice + as an empirical function of ice saturation ratio, + see eq. (5) in [Mohler2006](@cite). +```math +\begin{equation} +\frac{dn_{ice}}{dt} = N_{aer} a \frac{dS_i}{dt} +\end{equation} +``` +where: + - ``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 parameters ``a`` and ``S_0`` can be found in [Mohler2006](@cite). These + free parameters are 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 ``f_i`` values above 0.08 or ``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 +110,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 1117889b7..6be24cb69 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,20 @@ 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. + We only run it for the first GCM timestep because the rate approach requires + the change in ice saturation over time. With the discontinuous jump in saturation, + the parameterization is unable to determine a proper nucleation rate. 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 for the first GCM timestep. ```@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 95fe557f8..eba2d7049 100644 --- a/parcel/Tully_et_al_2023.jl +++ b/parcel/Tully_et_al_2023.jl @@ -74,76 +74,12 @@ 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 = # ice nucleation modes + [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]] 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 +96,110 @@ 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"] + #! format: off + 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 + #! format: on + + 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 8ed506fa2..c11b2b97d 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -75,16 +75,30 @@ 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 - AF = CMI_het.dust_activated_number_fraction( - aerosol, - ip.deposition, - S_i, - T, - ) + if "MohlerAF_Deposition" in ice_nucleation_modes + if S_i > ip.deposition.Sᵢ_max + println("Supersaturation exceeds the allowed value. No dust will be activated.") + AF = FT(0) + else + AF = CMI_het.dust_activated_number_fraction( + aerosol, + ip.deposition, + S_i, + T, + ) + end 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 + if S_i > ip.deposition.Sᵢ_max + println("Supersaturation exceeds the allowed value. No dust will be activated.") + dN_act_dt_depo = FT(0) + else + dN_act_dt_depo = CMI_het.MohlerDepositionRate(aerosol, ip.deposition, S_i, T, (dY[1] * ξ), N_aer) + end + 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) @@ -264,12 +278,18 @@ function run_parcel(IC, t_0, t_end, p) #! format: off FT = eltype(IC) (; const_dt, ice_nucleation_modes, growth_modes, droplet_size_distribution) = p + timestepper = ODE.Tsit5() 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.Euler() end if "ImmersionFreezing" in ice_nucleation_modes print("Immersion freezing ") @@ -306,7 +326,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), diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 2c6270d47..40044ba01 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -7,6 +7,7 @@ import ..Parameters as CMP import Thermodynamics as TD export dust_activated_number_fraction +export MohlerDepositionRate export deposition_J export ABIFM_J @@ -30,15 +31,39 @@ function dust_activated_number_fraction( T::FT, ) where {FT} - if Si > ip.Sᵢ_max - @warn "Supersaturation exceeds the allowed value." - @warn "No dust particles will be activated" - return FT(0) - else - S₀::FT = T > ip.T_thr ? dust.S₀_warm : dust.S₀_cold - a::FT = T > ip.T_thr ? dust.a_warm : dust.a_cold - return max(0, exp(a * (Si - S₀)) - 1) - end + @assert Si < ip.Sᵢ_max + + S₀::FT = T > ip.T_thr ? dust.S₀_warm : dust.S₀_cold + a::FT = T > ip.T_thr ? dust.a_warm : dust.a_cold + return max(0, exp(a * (Si - S₀)) - 1) +end + +""" + MohlerDepositionRate(dust, ip, S_i, T, dSi_dt, N_aer) + + - `dust` - a struct with dust parameters + - `ip` - a struct with ice nucleation parameters + - `Si` - ice saturation + - `T` - ambient temperature + - `dSi_dt` - change in ice saturation over time + - `N_aer` - number of unactivated aerosols + +Returns the ice nucleation rate from deposition. +From Mohler et al 2006 equation 5. +""" +function MohlerDepositionRate( + dust::Union{CMP.DesertDust, CMP.ArizonaTestDust}, + ip::CMP.Mohler2006, + Si::FT, + T::FT, + dSi_dt::FT, + N_aer::FT, +) where {FT} + + @assert Si < ip.Sᵢ_max + + a::FT = T > ip.T_thr ? dust.a_warm : dust.a_cold + return max(0, N_aer * a * dSi_dt) end """ diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 3dc7f7165..bdcb66e36 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -480,6 +480,62 @@ end end end +@kernel function IceNucleation_dust_activated_number_fraction_kernel!( + output::AbstractArray{FT}, + desert_dust, + arizona_test_dust, + ip, + S_i, + T, +) where {FT} + + i = @index(Group, Linear) + + @inbounds begin + output[1] = + CMI_het.dust_activated_number_fraction(desert_dust, ip, S_i, T) + output[2] = CMI_het.dust_activated_number_fraction( + arizona_test_dust, + ip, + S_i, + T, + ) + end +end + +@kernel function IceNucleation_MohlerDepositionRate_kernel!( + output::AbstractArray{FT}, + desert_dust, + arizona_test_dust, + ip, + S_i, + T, + dSi_dt, + N_aero, +) where {FT} + + i = @index(Group, Linear) + + @inbounds begin + output[1] = CMI_het.MohlerDepositionRate( + desert_dust, + ip, + S_i, + T, + dSi_dt, + N_aero, + ) + output[2] = CMI_het.MohlerDepositionRate( + arizona_test_dust, + ip, + S_i, + T, + dSi_dt, + N_aero, + ) + end +end + @kernel function IceNucleation_deposition_J_kernel!( output::AbstractArray{FT}, kaolinite, @@ -591,6 +647,8 @@ function test_gpu(FT) # ice nucleation H2SO4_prs = CMP.H2SO4SolutionParameters(FT) + desert_dust = CMP.DesertDust(FT) + arizona_test_dust = CMP.ArizonaTestDust(FT) illite = CMP.Illite(FT) kaolinite = CMP.Kaolinite(FT) feldspar = CMP.Feldspar(FT) @@ -840,6 +898,54 @@ function test_gpu(FT) end @testset "Ice Nucleation kernels" begin + dims = (1, 2) + (; output, ndrange) = setup_output(dims, FT) + + T = FT(240) + S_i = FT(1.2) + + kernel! = IceNucleation_dust_activated_number_fraction_kernel!( + backend, + work_groups, + ) + kernel!( + output, + desert_dust, + arizona_test_dust, + ip.deposition, + S_i, + T; + ndrange, + ) + + # test if dust_activated_number_fraction is callable and returns reasonable values + @test Array(output)[1] ≈ FT(0.0129835639) + @test Array(output)[2] ≈ FT(1.2233164999) + + dims = (1, 2) + (; output, ndrange) = setup_output(dims, FT) + + dSi_dt = FT(0.03) + N_aero = FT(3000) + + kernel! = + IceNucleation_MohlerDepositionRate_kernel!(backend, work_groups) + kernel!( + output, + desert_dust, + arizona_test_dust, + ip.deposition, + S_i, + T, + dSi_dt, + N_aero; + ndrange, + ) + + # test if MohlerDepositionRate is callable and returns reasonable values + @test Array(output)[1] ≈ FT(38.6999999999) + @test Array(output)[2] ≈ FT(423) + dims = (1, 3) (; output, ndrange) = setup_output(dims, FT) diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl index f6f63c74e..bb881c4d9 100644 --- a/test/heterogeneous_ice_nucleation_tests.jl +++ b/test/heterogeneous_ice_nucleation_tests.jl @@ -32,18 +32,9 @@ function test_heterogeneous_ice_nucleation(FT) Si_med = FT(1.2) Si_hgh = FT(1.34) Si_too_hgh = FT(1.5) - - # No activation below critical supersaturation - for dust in [ATD, desert_dust] - for T in [T_warm, T_cold] - TT.@test CMI_het.dust_activated_number_fraction( - dust, - ip.deposition, - Si_low, - T, - ) == FT(0) - end - end + dSi_dt = FT(0.05) + dSi_dt_negative = FT(-0.3) + N_aer = FT(3000) # Activate more in cold temperatures and higher supersaturations for dust in [ATD, desert_dust] @@ -69,15 +60,53 @@ function test_heterogeneous_ice_nucleation(FT) Si_med, T_warm, ) + TT.@test CMI_het.MohlerDepositionRate( + dust, + ip.deposition, + Si_med, + T_cold, + dSi_dt, + N_aer, + ) > CMI_het.MohlerDepositionRate( + dust, + ip.deposition, + Si_med, + T_warm, + dSi_dt, + N_aer, + ) end + # no activation if saturation exceeds allowed value for dust in [ATD, desert_dust] for T in [T_warm, T_cold] - TT.@test CMI_het.dust_activated_number_fraction( + TT.@test_throws AssertionError("Si < ip.Sᵢ_max") CMI_het.dust_activated_number_fraction( + dust, + ip.deposition, + Si_too_hgh, + T, + ) + TT.@test_throws AssertionError("Si < ip.Sᵢ_max") CMI_het.MohlerDepositionRate( dust, ip.deposition, Si_too_hgh, T, + dSi_dt, + N_aer, + ) + end + end + + # no activation if dSi_dt is negative + for dust in [ATD, desert_dust] + for T in [T_warm, T_cold] + TT.@test CMI_het.MohlerDepositionRate( + dust, + ip.deposition, + Si_low, + T, + dSi_dt_negative, + N_aer, ) == FT(0) end end diff --git a/test/performance_tests.jl b/test/performance_tests.jl index dfa04750d..842eda4a9 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -88,6 +88,7 @@ function benchmark_test(FT) T_air_2 = FT(250) T_air_cold = FT(230) S_ice = FT(1.2) + dSi_dt = FT(0.05) w_air = FT(0.5) p_air = FT(1e5) @@ -137,6 +138,11 @@ function benchmark_test(FT) (desert_dust, ip.deposition, S_ice, T_air_2), 50, ) + bench_press( + CMI_het.MohlerDepositionRate, + (desert_dust, ip.deposition, S_ice, T_air_2, dSi_dt, N_aer), + 80, + ) bench_press(CMI_het.deposition_J, (kaolinite, Delta_a_w), 230) bench_press(CMI_het.ABIFM_J, (desert_dust, Delta_a_w), 230) bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230)