Sa/add aerosol activation

.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

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
 

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,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(
+                cmp.aerosol,
+                Symbol(string(aerosol_name) * "_radius"),
+            )
+            seasalt_particle_mass =
+                4 / 3 *
+                FT(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 =
+                4 / 3 *
+                FT(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
+    @. 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,
+        (cmp,),
+        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,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, cloud liquid and cloud ice specific humidity
+- `nₗ`, - cloud liquid 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 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
+    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ₗ) / 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/).