implement ARG aerosol activation with prescribed aerosol densities for 2M microphysics

Author: sajjadazimi

.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"]

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"

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

src/parameterized_tendencies/microphysics/cloud_condensate.jl

@@ -105,19 +105,67 @@ 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(
+                cmc.aerosol,
+                Symbol(string(aerosol_name) * "_radius"),
+            )
+            seasalt_particle_mass =
+                4 / 3 *
+                FT(pi) *
+                seasalt_particle_radius^3 *
+                cmc.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 =
+                4 / 3 *
+                FT(pi) *
+                cmc.aerosol.sulfate_radius^3 *
+                cmc.aerosol.sulfate_density
+            sulfate_mass = getproperty(aerosol_field, aerosol_name)
+            @. sulfate_num += sulfate_mass / sulfate_particle_mass
+        end
+    end
+    @. seasalt_mean_radius = 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.ρ),
+        specific(Y.c.ρq_ice, Y.c.ρ),
+        specific(Y.c.ρn_liq, Y.c.ρ),
+        Y.c.ρ,
+        ᶜw,
+        (cmc,),
+        thp,
+        ᶜts,
+        dt,
+    )
+    @. Yₜ.c.ρn_liq += Y.c.ρ * Snₗ
 end

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,130 @@ 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,
+        nₗ,
+        ρ,
+        w,
+        air_params,
+        aerosol_params,
+        arg_params,
+        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, cloud liquid and cloud ice specific humidity
+- `nₗ`, - cloud liquid number concentration per mass [kg⁻¹]
+- `ρ`: Air density [kg/m³].
+- `w`: Vertical velocity [m/s].
+- `air_params`: Struct containing physical constants related to air.
+- `aerosol_params`: Struct containing aerosol properties, such as hygroscopicity parameters (κ).
+- `arg_params`: Activation scheme parameters specific to Abdul-Razzak and Ghan (2000).
+- `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 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,
+    cmc,
+    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 = cmc.aps
+    arg_params = cmc.arg
+    aerosol_params = cmc.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
+    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),
+        triangle_inequality_limiter(n_act - nₗ, n_aer - nₗ) / 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/).

src/parameters/create_parameters.jl

@@ -135,6 +135,9 @@ cloud_parameters(toml_dict::CP.AbstractTOMLDict) = (;
         "ClimaAtmos",
     ).prescribed_cloud_droplet_number_concentration,
     aml = aerosol_ml_parameters(toml_dict),
+    aps = CM.Parameters.AirProperties(toml_dict),
+    arg = CM.Parameters.AerosolActivationParameters(toml_dict),
+    aerosol = prescribed_aerosol_parameters(toml_dict),
 )
 
 microphys_1m_parameters(::Type{FT}) where {FT <: AbstractFloat} =
@@ -195,6 +198,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)