ARG aerosol activation applied at cloud base

Author: sajjadazimi

Project.toml

@@ -44,7 +44,7 @@ AtmosphericProfilesLibrary = "0.1.7"
 ClimaComms = "0.6.8"
 ClimaCore = "0.14.34"
 ClimaDiagnostics = "0.2.12"
-ClimaParams = "0.10.32"
+ClimaParams = "0.10.33"
 ClimaTimeSteppers = "0.8.2"
 ClimaUtilities = "0.1.22"
 CloudMicrophysics = "0.24.1"

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/cache/temporary_quantities.jl

@@ -22,6 +22,10 @@ function temporary_quantities(Y, atmos)
         temp_field_level = Fields.level(Fields.Field(FT, center_space), 1),
         temp_field_level_2 = Fields.level(Fields.Field(FT, center_space), 1),
         temp_field_level_3 = Fields.level(Fields.Field(FT, center_space), 1),
+        temp_tuple_field_level = similar(
+            Fields.level(Fields.Field(FT, center_space), 1),
+            Tuple{FT, FT},
+        ),
         temp_data = Fields.field_values(Fields.Field(FT, center_space)),
         temp_data_face_level = Fields.field_values(
             Fields.level(Fields.Field(FT, face_space), half),

src/parameterized_tendencies/microphysics/cloud_condensate.jl

@@ -105,19 +105,70 @@ 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 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.ρn_liq, Y.c.ρ),
+        Y.c.ρ,
+        ᶜw,
+        (cmc,),
+        thp,
+        ᶜts,
+        dt,
+    )
+    # Use the Snₗ tendnecy at cloud base
+    cloud_base_Snₗ_and_z = p.scratch.temp_tuple_field_level
+    z = p.scratch.ᶜtemp_scalar
+    z = Fields.coordinate_field(Snₗ).z
+    find_cloud_base!(Snₗ, z, cloud_base_Snₗ_and_z)
+    @. Snₗ = ifelse(z == last(cloud_base_Snₗ_and_z), Snₗ, FT(0))
+    @. Yₜ.c.ρn_liq += Y.c.ρ * Snₗ
 end

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -1,10 +1,13 @@
 # A set of wrappers for using CloudMicrophysics.jl functions inside EDMFX loops
 
 import Thermodynamics as TD
+import ClimaCore.Operators as Operators
 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
@@ -25,33 +28,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)
 
@@ -373,46 +349,129 @@ end
 #####
 
 """
-    aerosol_activation_sources(cm_params, thp, ρ, Tₐ, qₜ, qₗ, qᵢ, qᵣ, qₛ, n_dp, n_dp_prescribed, dt)
+    find_cloud_base!(Snₗ, z, cloud_base_Snₗ_and_z)
 
- - 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
+Finds the cloud base in each vertical column based on the first occurrence of aerosol activation.
+
+This function performs a vertical column-wise reduction to identify the first level (from bottom to
+top) where `Snₗ > 0`, indicating the onset of cloud droplet activation (i.e., cloud base). The result
+is stored in `cloud_base_Snₗ_and_z` as a tuple `(Snₗ, z)` for each column.
+
+# Arguments
+- `Snₗ`: Field of cloud droplet number concentration tendencies [m⁻³/s], from an activation parameterization.
+- `z`: Field of geometric height [m] corresponding to each grid point.
+- `cloud_base_Snₗ_and_z`: Output field to store the `(Snₗ, z)` values at cloud base for each column.
+"""
+function find_cloud_base!(Snₗ, z, cloud_base_Snₗ_and_z)
+    Operators.column_reduce!(
+        ((target_Snₗ_value, target_z_value), (Snₗ_value, z_value)) -> ifelse(
+            target_Snₗ_value == 0 && Snₗ_value > 0,
+            (Snₗ_value, z_value),
+            (target_Snₗ_value, target_z_value),
+        ), # reduction function
+        cloud_base_Snₗ_and_z, # destination for output (a field of tuples)
+        Base.broadcasted(tuple, Snₗ, z),
+    )
+end
 
-Returns the activation rate. #TODO This function temporarily computes activation rate
-based on mass rates and a prescribed droplet mass (no activation parameterization yet).
+"""
+    aerosol_activation_sources(
+        seasalt_num,
+        seasalt_mean_radius,
+        sulfate_num,
+        nₗ,
+        ρ,
+        w,
+        air_params,
+        aerosol_params,
+        arg_params,
+        thermo_params,
+        ts,
+        dt,
+    )
+
+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⁻¹].
+- `nₗ`: Cloud droplet 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,
+    seasalt_num,
+    seasalt_mean_radius,
+    sulfate_num,
+    nₗ,
     ρ,
-    Tₐ,
-    qₜ,
-    qₗ,
-    qᵢ,
-    qᵣ,
-    qₛ,
-    n_dp,
-    n_dp_prescribed,
+    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
+    q = TD.PhasePartition(thermo_params, ts)
+    T = TD.air_temperature(thermo_params, ts)
+    p = TD.air_pressure(thermo_params, ts)
+    S::FT = TD.supersaturation(thermo_params, q, ρ, T, TD.Liquid())
+
+    seasalt_mode = CMAM.Mode_κ(
+        seasalt_mean_radius,
+        aerosol_params.seasalt_std,
+        seasalt_num * ρ,
+        (FT(1.0),),
+        (FT(1.0),),
+        (FT(0.0),),
+        (aerosol_params.seasalt_kappa,),
+    )
+    sulfate_mode = CMAM.Mode_κ(
+        aerosol_params.sulfate_radius,
+        aerosol_params.sulfate_std,
+        sulfate_num * ρ,
+        (FT(1.0),),
+        (FT(1.0),),
+        (FT(0.0),),
+        (aerosol_params.sulfate_kappa,),
+    )
+
+    aerosol_dist = CMAM.AerosolDistribution((seasalt_mode, sulfate_mode))
+    n_act::FT =
+        CMAA.total_N_activated(
+            arg_params,
+            aerosol_dist,
+            air_params,
+            thermo_params,
+            T,
+            p,
+            w,
+            q,
+        ) / ρ
 
     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 < 0 || isnan(n_act),
+        FT(0),
+        triangle_inequality_limiter(n_act, n_act - nₗ) / dt,
     )
 end
 

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)