diff --git a/.buildkite/Manifest-v1.11.toml b/.buildkite/Manifest-v1.11.toml index 891d84bf96..e0759210e8 100644 --- a/.buildkite/Manifest-v1.11.toml +++ b/.buildkite/Manifest-v1.11.toml @@ -461,10 +461,10 @@ uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.13" [[deps.CloudMicrophysics]] -deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "LazyArtifacts", "LogExpFunctions", "QuadGK", "RootSolvers", "SpecialFunctions", "Thermodynamics"] -git-tree-sha1 = "cfc9decc96d616ab3b8d911f843eeb2aafd7d2ad" +deps = ["ClimaParams", "DocStringExtensions", "ForwardDiff", "LazyArtifacts", "LogExpFunctions", "QuadGK", "RootSolvers", "SpecialFunctions", "StaticArrays", "Thermodynamics"] +git-tree-sha1 = "8f7137d8bb3098d6c46b2c3f859fb873d5b28abd" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" -version = "0.25.0" +version = "0.26.1" [deps.CloudMicrophysics.extensions] EmulatorModelsExt = ["DataFrames", "MLJ"] diff --git a/Project.toml b/Project.toml index bdefca8f3f..4611b5f871 100644 --- a/Project.toml +++ b/Project.toml @@ -49,7 +49,7 @@ ClimaInterpolations = "0.1.0" ClimaParams = "0.10.35" ClimaTimeSteppers = "0.8.2" ClimaUtilities = "0.1.22" -CloudMicrophysics = "0.25.0" +CloudMicrophysics = "0.26.1" Dates = "1" ForwardDiff = "0.10.38, 1" Insolation = "0.9.2" @@ -64,10 +64,10 @@ NullBroadcasts = "0.1" RRTMGP = "0.21.3" Random = "1" SciMLBase = "2.12" -StaticArrays = "1.7" +StaticArrays = "1.9" Statistics = "1" SurfaceFluxes = "0.11, 0.12" -Thermodynamics = "0.12.13" +Thermodynamics = "0.12.14" UnrolledUtilities = "0.1.9" YAML = "0.4" julia = "1.9" diff --git a/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml b/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml index 2b577ec244..e8f83c1e07 100644 --- a/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml +++ b/config/model_configs/aquaplanet_nonequil_allsky_gw_res_2M.yml @@ -20,7 +20,7 @@ orographic_gravity_wave: "gfdl_restart" surface_setup: "DefaultMoninObukhov" prescribe_ozone: true reproducibility_test: true -prescribed_aerosols: ["SO4", "CB1", "OC1", "DST01", "SSLT01"] +prescribed_aerosols: ["SO4", "CB1", "OC1", "DST01", "SSLT01", "SSLT02", "SSLT03", "SSLT04", "SSLT05"] diagnostics: - short_name: [edt, evu, mmrso4, mmrbcpo, mmrocpo, mmrdust, mmrss, loadss, o3, od550aer, odsc550aer] reduction_time: average diff --git a/docs/src/equations.md b/docs/src/equations.md index 2358177ce0..b9f9c3ffad 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -414,3 +414,14 @@ The source terms functions treat negative specific humidities as zeros, so the simulations should be stable even with small negative numbers. We do not apply hyperdiffusion for precipitation tracers. + +### Aerosol Activation for 2-Moment Microphysics + +Aerosol activation uses functions from the [CloudMicrophysics.jl](https://github.com/CliMA/CloudMicrophysics.jl) library, based on the Abdul-Razzak and Ghan (ARG) parameterization. ARG predicts the number of activated cloud droplets assuming a parcel of clear air rising adiabatically. This formulation is traditionally applied only at cloud base, where the maximum supersaturation typically occurs. + +To enable ARG to be used locally (i.e., without explicitly identifying cloud base), CloudMicrophysics.jl implements a modified equation for the maximum supersaturation that accounts for the presence of pre-existing liquid and ice particles. This allows activation to be applied inside clouds. To ensure that activation occurs only where physically appropriate, we apply additional clipping logic: + +- If the predicted maximum supersaturation is less than the local supersaturation (i.e., supersaturation is decreasing), aerosol activation is not applied. +- If the predicted number of activated droplets is less than the existing local cloud droplet number concentration, activation is also suppressed. + +This ensures that droplet activation occurs only in physically meaningful regions—typically near cloud base—even though the activation routine can be applied throughout the domain. diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl index 6491e2a62e..20b6cba58a 100644 --- a/reproducibility_tests/ref_counter.jl +++ b/reproducibility_tests/ref_counter.jl @@ -1,4 +1,4 @@ -248 +249 # **README** # @@ -20,6 +20,9 @@ #= +249 +- Add aerosol activation for 2M microphysics + 248 - update deps: climacore 0.14.35 diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl index fd0fd6e385..3fa4cfde95 100644 --- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl +++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl @@ -74,13 +74,13 @@ function cloud_condensate_tendency!( (; params, dt) = p FT = eltype(params) thp = CAP.thermodynamics_params(params) - cmc = CAP.microphysics_cloud_params(params) + cmp = CAP.microphysics_2m_params(params) Tₐ = @. lazy(TD.air_temperature(thp, ᶜts)) @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources( - cmc.liquid, + cmp.liquid, thp, specific(Y.c.ρq_tot, Y.c.ρ), specific(Y.c.ρq_liq, Y.c.ρ), @@ -93,7 +93,7 @@ function cloud_condensate_tendency!( ) @. Yₜ.c.ρq_ice += Y.c.ρ * cloud_sources( - cmc.ice, + cmp.ice, thp, specific(Y.c.ρq_tot, Y.c.ρ), specific(Y.c.ρq_liq, Y.c.ρ), @@ -105,19 +105,68 @@ function cloud_condensate_tendency!( dt, ) - @. Yₜ.c.ρn_liq += - Y.c.ρ * aerosol_activation_sources( - cmc.liquid, - thp, - Y.c.ρ, - Tₐ, - specific(Y.c.ρq_tot, Y.c.ρ), - specific(Y.c.ρq_liq, Y.c.ρ), - specific(Y.c.ρq_ice, Y.c.ρ), - specific(Y.c.ρq_rai, Y.c.ρ), - specific(Y.c.ρq_sno, Y.c.ρ), - specific(Y.c.ρn_liq, Y.c.ρ), - cmc.N_cloud_liquid_droplets, - dt, - ) + # Aerosol activation using prescribed aerosol (Sea salt and sulfate) + if !(:prescribed_aerosols_field in propertynames(p.tracers)) + return + end + seasalt_num = p.scratch.ᶜtemp_scalar + seasalt_mean_radius = p.scratch.ᶜtemp_scalar_2 + sulfate_num = p.scratch.ᶜtemp_scalar_3 + @. seasalt_num = 0 + @. seasalt_mean_radius = 0 + @. sulfate_num = 0 + # Get aerosol concentrations if available + seasalt_names = [:SSLT01, :SSLT02, :SSLT03, :SSLT04, :SSLT05] + sulfate_names = [:SO4] + aerosol_field = p.tracers.prescribed_aerosols_field + for aerosol_name in propertynames(aerosol_field) + if aerosol_name in seasalt_names + seasalt_particle_radius = getproperty( + cmp.aerosol, + Symbol(string(aerosol_name) * "_radius"), + ) + seasalt_particle_mass = + FT(4 / 3 * pi) * + seasalt_particle_radius^3 * + cmp.aerosol.seasalt_density + seasalt_mass = getproperty(aerosol_field, aerosol_name) + @. seasalt_num += seasalt_mass / seasalt_particle_mass + @. seasalt_mean_radius += + seasalt_mass / seasalt_particle_mass * + log(seasalt_particle_radius) + elseif aerosol_name in sulfate_names + sulfate_particle_mass = + FT(4 / 3 * pi) * + cmp.aerosol.sulfate_radius^3 * + cmp.aerosol.sulfate_density + sulfate_mass = getproperty(aerosol_field, aerosol_name) + @. sulfate_num += sulfate_mass / sulfate_particle_mass + end + end + # Compute geometric mean radius of the log-normal distribution: + # exp(weighted average of log(radius)) + @. seasalt_mean_radius = + ifelse(seasalt_num == 0, 0, exp(seasalt_mean_radius / seasalt_num)) + + # Compute aerosol activation (ARG 2000) + (; ᶜu) = p.precomputed + ᶜw = p.scratch.ᶜtemp_scalar_4 + Snₗ = p.scratch.ᶜtemp_scalar_5 + @. ᶜw = max(0, w_component.(Geometry.WVector.(ᶜu))) + @. Snₗ = aerosol_activation_sources( + seasalt_num, + seasalt_mean_radius, + sulfate_num, + specific(Y.c.ρq_tot, Y.c.ρ), + specific(Y.c.ρq_liq + Y.c.ρq_rai, Y.c.ρ), + specific(Y.c.ρq_ice + Y.c.ρq_sno, Y.c.ρ), + specific(Y.c.ρn_liq + Y.c.ρn_rai, Y.c.ρ), + Y.c.ρ, + ᶜw, + (cmp,), + thp, + ᶜts, + dt, + ) + @. Yₜ.c.ρn_liq += Y.c.ρ * Snₗ end diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 84a6022f16..f457d1932b 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -1,14 +1,18 @@ # A set of wrappers for using CloudMicrophysics.jl functions inside EDMFX loops import Thermodynamics as TD +import CloudMicrophysics.ThermodynamicsInterface as CMTDI import CloudMicrophysics.Microphysics0M as CM0 import CloudMicrophysics.Microphysics1M as CM1 import CloudMicrophysics.Microphysics2M as CM2 import CloudMicrophysics.MicrophysicsNonEq as CMNe +import CloudMicrophysics.AerosolModel as CMAM +import CloudMicrophysics.AerosolActivation as CMAA import CloudMicrophysics.Parameters as CMP # Define some aliases and functions to make the code more readable const Tₐ = TD.air_temperature +const Pₐ = TD.air_pressure const PP = TD.PhasePartition const qᵥ = TD.vapor_specific_humidity @@ -25,33 +29,6 @@ function limit(q, dt, n::Int) return max(FT(0), q) / float(dt) / n end -""" - ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq) - - - cmc - a struct with cloud and aerosol parameters - - c_dust, c_seasalt, c_SO4 - dust, seasalt and ammonium sulfate mass concentrations [kg/kg] - - q_liq - liquid water specific humidity - -Returns the liquid cloud droplet number concentration diagnosed based on the -aerosol loading and cloud liquid water. -""" -function ml_N_cloud_liquid_droplets(cmc, c_dust, c_seasalt, c_SO4, q_liq) - # We can also add w, T, RH, w' ... - # Also consider lookind only at around cloud base height - (; α_dust, α_seasalt, α_SO4, α_q_liq) = cmc.aml - (; c₀_dust, c₀_seasalt, c₀_SO4, q₀_liq) = cmc.aml - N₀ = cmc.N_cloud_liquid_droplets - - FT = eltype(N₀) - return N₀ * ( - FT(1) + - α_dust * (log(max(c_dust, eps(FT))) - log(c₀_dust)) + - α_seasalt * (log(max(c_seasalt, eps(FT))) - log(c₀_seasalt)) + - α_SO4 * (log(max(c_SO4, eps(FT))) - log(c₀_SO4)) + - α_q_liq * (log(max(q_liq, eps(FT))) - log(q₀_liq)) - ) -end - """ cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt) @@ -386,60 +363,129 @@ end ##### """ - aerosol_activation_sources(cm_params, thp, ρ, Tₐ, qₜ, qₗ, qᵢ, qᵣ, qₛ, n_dp, n_dp_prescribed, dt) - - - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate - - thp - Thermodynamics parameters struct - - ρ - air density - - Tₐ - air temperature - - qₜ - total specific humidity - - qₗ - liquid specific humidity - - qᵢ - ice specific humidity - - qᵣ - rain specific humidity - - qₛ - snow specific humidity - - n_dp - number concentration droplets (liquid or ice) per mass - _ n_dp_prescribed - prescribed number concentration of droplets (liquid or ice) per mass - - dt - model time step + aerosol_activation_sources( + seasalt_num, + seasalt_mean_radius, + sulfate_num, + qₜ, + qₗ, + qᵢ, + nₗ, + ρ, + w, + cmp, + thermo_params, + ts, + dt, + ) -Returns the activation rate. #TODO This function temporarily computes activation rate -based on mass rates and a prescribed droplet mass (no activation parameterization yet). +Computes the source term for cloud droplet number concentration per mass due to aerosol activation, +based on the Abdul-Razzak and Ghan (2000) parameterization. + +This function estimates the number of aerosols activated into cloud droplets per mass of air per second +from a bi-modal aerosol distribution (sea salt and sulfate), given local supersaturation and vertical +velocity. The result is returned as a tendency (per second) of liquid droplet number concentration. + +# Arguments +- `seasalt_num`: Number concentration per mass of sea salt aerosols [kg⁻¹]. +- `seasalt_mean_radius`: Mean dry radius of sea salt aerosol mode [m]. +- `sulfate_num`: Number concentration per mass of sulfate aerosols [kg⁻¹]. +- `qₜ`, `qₗ`, `qᵢ` - total water, liquid (cloud liquid and rain) and ice (cloud ice and snow) specific humidity +- `nₗ` - lqiuid (cloud liquid and rain) number concentration per mass [kg⁻¹] +- `ρ`: Air density [kg/m³]. +- `w`: Vertical velocity [m/s]. +- `cmp`: Microphysics parameters +- `thermo_params`: Thermodynamic parameters for computing saturation, pressure, temperature, etc. +- `ts`: Thermodynamic state (e.g., prognostic variables) used for evaluating the phase partition. +- `dt`: Time step (s) over which the activation tendency is applied. + +# Returns +- Tendency of cloud liquid droplet number concentration per mass of air due to aerosol activation [m⁻³/s]. """ function aerosol_activation_sources( - cm_params::CMP.CloudLiquid{FT}, - thp, - ρ, - Tₐ, + seasalt_num, + seasalt_mean_radius, + sulfate_num, qₜ, qₗ, qᵢ, - qᵣ, - qₛ, - n_dp, - n_dp_prescribed, + nₗ, + ρ, + w, + cmp, + thermo_params, + ts, dt, -) where {FT} - r_dp = FT(2e-6) # 2 μm - m_dp = 4 / 3 * FT(π) * r_dp^3 * cm_params.ρw - Sn = cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt) / m_dp +) + + FT = eltype(nₗ) + air_params = cmp.aps + arg_params = cmp.arg + aerosol_params = cmp.aerosol + T = Tₐ(thermo_params, ts) + p = Pₐ(thermo_params, ts) + S = CMTDI.supersaturation_over_liquid(thermo_params, qₜ, qₗ, qᵢ, ρ, T) + n_aer = seasalt_num + sulfate_num + if S < 0 || n_aer < eps(FT) + return FT(0) + end + + seasalt_mode = CMAM.Mode_κ( + seasalt_mean_radius, + aerosol_params.seasalt_std, + seasalt_num * ρ, + (FT(1),), + (FT(1),), + (FT(0),), + (aerosol_params.seasalt_kappa,), + ) + sulfate_mode = CMAM.Mode_κ( + aerosol_params.sulfate_radius, + aerosol_params.sulfate_std, + sulfate_num * ρ, + (FT(1),), + (FT(1),), + (FT(0),), + (aerosol_params.sulfate_kappa,), + ) + + aerosol_dist = CMAM.AerosolDistribution((seasalt_mode, sulfate_mode)) + + args = ( + arg_params, + aerosol_dist, + air_params, + thermo_params, + T, + p, + w, + qₜ, + qₗ, + qᵢ, + ρ * nₗ, + FT(0), + ) #assuming no ice particles because we don't track n_ice for now + S_max = CMAA.max_supersaturation(args...) + n_act = CMAA.total_N_activated(args...) / ρ return ifelse( - Sn > FT(0), - triangle_inequality_limiter(Sn, limit((n_dp_prescribed - n_dp), dt, 2)), - -triangle_inequality_limiter(abs(Sn), limit(n_dp, dt, 2)), + S_max < S || isnan(n_act) || n_act < nₗ, + FT(0), + (n_act - nₗ) / float(dt), ) end """ compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, nₗ, nᵣ, qₜ, qₗ, qᵢ, qᵣ, qₛ, ts, dt, sb, thp) - - Sᵖ, S₂ᵖ - temporary containters to help compute precipitation source terms - - Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ - cached storage for precipitation source terms - - ρ - air density - - nₗ, nᵣ - cloud liquid and rain number concentration per mass [1 / kg of moist air] - - qₜ, qₗ, qᵢ, qᵣ, qₛ - total water, cloud liquid, cloud ice, rain and snow specific humidity - - ts - thermodynamic state (see td package for details) - - dt - model time step - - thp, mp - structs with thermodynamic and microphysics parameters + - `Sᵖ`, `S₂ᵖ` - temporary containters to help compute precipitation source terms + - `Snₗᵖ`, `Snᵣᵖ`, `Sqₗᵖ`, `Sqᵣᵖ` - cached storage for precipitation source terms + - `ρ` - air density + - `nₗ`, `nᵣ` - cloud liquid and rain number concentration per mass [1 / kg of moist air] + - `qₜ`, `qₗ`, `qᵢ`, `qᵣ`, `qₛ` - total water, cloud liquid, cloud ice, rain and snow specific humidity + - `ts` - thermodynamic state (see td package for details) + - `dt` - model time step + - `thp`, `mp` - structs with thermodynamic and microphysics parameters Computes precipitation number and mass sources due to warm precipitation processes based on the 2-moment [Seifert and Beheng (2006) scheme](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics2M/). diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index 117bef9636..b855da66a5 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -163,6 +163,10 @@ microphys_2m_parameters(toml_dict::CP.AbstractTOMLDict) = (; sb = CM.Parameters.SB2006(toml_dict), aps = CM.Parameters.AirProperties(toml_dict), tv = CM.Parameters.SB2006VelType(toml_dict), + liquid = CM.Parameters.CloudLiquid(toml_dict), + ice = CM.Parameters.CloudIce(toml_dict), + arg = CM.Parameters.AerosolActivationParameters(toml_dict), + aerosol = prescribed_aerosol_parameters(toml_dict), ) function vert_diff_parameters(toml_dict) @@ -195,6 +199,24 @@ function aerosol_ml_parameters(toml_dict) return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") end +function prescribed_aerosol_parameters(toml_dict) + name_map = (; + :MERRA2_seasalt_aerosol_bin01_radius => :SSLT01_radius, + :MERRA2_seasalt_aerosol_bin02_radius => :SSLT02_radius, + :MERRA2_seasalt_aerosol_bin03_radius => :SSLT03_radius, + :MERRA2_seasalt_aerosol_bin04_radius => :SSLT04_radius, + :MERRA2_seasalt_aerosol_bin05_radius => :SSLT05_radius, + :seasalt_aerosol_kappa => :seasalt_kappa, + :seasalt_aerosol_density => :seasalt_density, + :mam3_stdev_coarse => :seasalt_std, + :MERRA2_sulfate_aerosol_radius => :sulfate_radius, + :sulfate_aerosol_kappa => :sulfate_kappa, + :sulfate_aerosol_density => :sulfate_density, + :mam3_stdev_accum => :sulfate_std, + ) + return CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") +end + to_svec(x::AbstractArray) = SA.SVector{length(x)}(x) to_svec(x) = x to_svec(x::NamedTuple) = map(x -> to_svec(x), x)