Implement local 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

@@ -48,8 +48,8 @@ ClimaDiagnostics = "0.2.12"
 ClimaInterpolations = "0.1.0"
 ClimaParams = "0.10.35"
 ClimaTimeSteppers = "0.8.2"
-ClimaUtilities = "0.1.22"
-CloudMicrophysics = "0.25.0"
+ClimaUtilities = "0.1.23"
+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

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.

perf/flame.jl

@@ -47,7 +47,7 @@ allocs_limit["flame_aquaplanet_progedmf"] = 774_712
 allocs_limit["flame_aquaplanet_progedmf_dense_autodiff"] = 774_712
 allocs_limit["flame_diffusion"] = 138_432
 allocs_limit["flame_threaded"] = 2047_736
-allocs_limit["flame_callbacks"] = 391_864
+allocs_limit["flame_callbacks"] = 391_942
 allocs_limit["flame_gravity_wave"] = 581_381_976
 # Ideally, we would like to track all the allocations, but this becomes too
 # expensive there is too many of them. Here, we set the default sample rate to

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-249
+250
 
 # **README**
 #
@@ -20,6 +20,11 @@
 
 
 #=
+250
+- Add ARG aerosol activation for 2M microphysics; increase the allocation limit for `flame_callbacks`
+  from 391864 to 391942 to account for additional allocations when constructing `ClimaAtmosParameters` 
+  with the new aerosol parameters in `microphysics_2m_parameters`.
+
 249
 - Remove viscosity, diffusivity, and mixing length from precomputed quantities
 

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