Modify supersaturation equation to include preexisting droplets; update alpha and gamma per Korolev & Mazin; update parcel example and docs (#615)

Author: sajjadazimi

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.26.0"
+version = "0.26.1"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

docs/src/AerosolActivation.md

@@ -140,11 +140,9 @@ where:
 
 ## Maximum Supersaturation
 
-Maximum supersaturation reached by the system ``S_{max}`` governs
-  what aerosol sizes are activated and what aren't.
-We estimate ``S_{max}`` by considering a parcel of air raising adiabatically
-  with a constant velocity, see for example [Rogers1975](@cite).
-The time rate of change of ``S`` is given by eq (10) in [Rogers1975](@cite)
+Maximum supersaturation reached by the system ``S_{max}`` governs what aerosol sizes are activated and what aren't.
+We estimate ``S_{max}`` by considering a parcel of air raising adiabatically with a constant velocity, see for example [KorolevMazin2003](@cite).
+The time rate of change of ``S`` is given by eq (A11) in [KorolevMazin2003](@cite)
 
 ```math
 \begin{equation}
@@ -157,25 +155,24 @@ where:
   - ``d\chi / dt`` is the water condensation rate during aerosol activation and growth,
   - ``\alpha`` and ``\gamma`` are coefficients that do not depend on aerosol properties.
 
-The parameters ``\alpha`` and ``\gamma`` are defined by eq. (11) and (12)
-  in [Abdul-Razzaketal1998](@cite) .
-Here, they are implemented following the [Rogers1975](@cite) eq. (10):
+The parameters ``\alpha`` and ``\gamma`` are defined similarly to eq. (11) and (12) in [Abdul-Razzaketal1998](@cite), but we follow the full formulation from [KorolevMazin2003](@cite) eq. (A11), without assuming ``S + 1 ≈ 1``. While the difference is often small, we retain the full expressions for completeness and consistency with our implementation. The parameters are:
 ```math
 \begin{equation}
-\alpha = \frac{L_v \, g \, \epsilon}{R_m \, c_{pm} \, T^2} - \frac{g}{R_m \, T}
+\alpha = \frac{p_{vap}}{p_{vap}^{sat}}\left[\frac{L_v \, g}{R_v \, c_{pm} \, T^2} - \frac{g}{R_m \, T}\right]
 \end{equation}
 ```
 ```math
 \begin{equation}
-\gamma = \frac{R_m T}{\epsilon \, p_{vap}^{sat}} + \frac{\epsilon \, L_v^2}{c_{pm} \, T \, p}
+\gamma = \frac{R_v T}{p_{vap}^{sat}} + \frac{p_{vap}}{p_{vap}^{sat}}\frac{R_m \, L_v^2}{R_v \, c_{pm} \, T \, p}
 \end{equation}
 ```
 where:
   - ``L_v`` is the latent heat of vaporization,
   - ``g`` is gravitational acceleration,
-  - ``\epsilon`` is the ratio of water molar mass to dry air molar mass,
+  - ``R_v`` is the gas constant of vapor,
   - ``R_m`` is the gas constant of air (dry air and moisture)
   - ``c_{pm}`` is the specific heat of air (dry air and moisture),
+  - ``p_{vap}`` is the vapor pressure,
   - ``p_{vap}^{sat}`` is the saturation vapor pressure,
   - ``p`` is the total air pressure.
 
@@ -285,3 +282,108 @@ The CloudMicrophysics package offers an advanced feature for predicting the aero
 In addition, we provide functionality and a demonstration for re-training the free parameters of the traditional ARG activation scheme using Ensemble Kalman Processes. This functionality is showcased in `test/aerosol_activation_calibration.jl`.
 
 Using ML emulators or calibration with Ensemble Kalman Processes for predicting aerosol activation is provided through an extension to the main package. This extension will be loaded with `CloudMicrophysics.jl` if both `MLJ.jl`, `DataFrames.jl`, and `EnsembleKalmanProcesses.jl` are loaded by the user as well.
+
+## Modification for Local Supersaturation with Preexisting Hydrometeors
+
+The original ARG activation scheme assumes adiabatic ascent with no preexisting cloud or ice particles, which is a valid approximation near cloud base. As a result, the ARG formulation is traditionally applied only at cloud base, where no vapor has yet condensed and no hydrometeors exist to affect the supersaturation evolution.
+
+However, in many applications—such as large-eddy simulations or cloud-resolving models—it is desirable to apply the activation scheme locally, throughout the vertical extent of the cloud. In such cases, preexisting cloud droplets and ice crystals act as persistent vapor sinks, suppressing supersaturation and inhibiting further activation.
+
+To relax the cloud-base-only constraint and extend ARG into a local operator, we propose modifying the supersaturation tendency equation by incorporating the effects of these preexisting hydrometeors.
+When coupling this approach to a full model, additional clipping may be needed: the modified ARG should only be applied when the local supersaturation does not exceed the predicted maximum, and when the number concentration of hydrometeors remains below the potential activation of aerosols.
+
+
+### Modified Supersaturation Budget
+
+We modify the supersaturation tendency equation to include additional vapor sink terms:
+
+```math
+\frac{dS}{dt} = \alpha w - \gamma \frac{d\chi}{dt}_{\text{activation}} - \gamma \frac{d\chi}{dt}_{\text{liquid}} - \gamma_i \frac{d\chi}{dt}_{\text{ice}},
+```
+
+where the thermodynamic coefficient ``\gamma_i`` is defined as
+```math
+\begin{equation}
+\gamma_i = \frac{R_v T}{p_{vap}^{sat}} + \frac{p_{vap}}{p_{vap}^{sat}}\frac{R_m \, L_v L_s}{R_v \, c_{pm} \, T \, p}
+\end{equation}
+```
+where ``L_s`` is the latent heat of sublimation. The new sink terms are:
+
+- Liquid sink:
+```math
+\frac{d\chi}{dt}_{\text{liquid}} = 4 \pi \rho_w N_l r_l G_l S,
+```
+
+- Ice sink:
+```math
+\frac{d\chi}{dt}_{\text{ice}} = 4 \pi \rho_i N_i r_i G_i S_i,
+```
+
+with
+```math
+S_i = \xi (S + 1) - 1, \quad \xi = \frac{e_{s,\ell}}{e_{s,i}}.
+```
+
+Here:
+
+- ``N_l``, ``r_l``, and ``G_l`` are the number concentration, radius, and condensational growth coefficient of liquid droplets,
+- ``N_i``, ``r_i``, and ``G_i`` are the corresponding values for ice crystals,
+- ``\rho_w`` and ``\rho_i`` are the densities of liquid water and ice.
+
+We define:
+```math
+K_l = 4 \pi \rho_w N_l r_l G_l \gamma, \quad K_i = 4 \pi \rho_i N_i r_i G_i \gamma_i,
+```
+
+yielding the modified supersaturation equation at peak supersaturation:
+```math
+0 = \alpha w - \gamma f(S_{max}) - K_l S_{max} - K_i \big[ \xi (S_{max} + 1) - 1 \big].
+```
+
+### Approximate Solution
+
+Let ``S_0`` be the ARG solution obtained by solving:
+```math
+\alpha w = \gamma f(S_0),
+```
+and approximate the nonlinear sink function ``f(S)`` as linear between 0 and ``S_0``:
+```math
+f(S) \approx f(S_0) \cdot \frac{S}{S_0}.
+```
+
+Substituting into the modified budget and solving for ``S_1``:
+```math
+\alpha w = \gamma f(S_0) \cdot \frac{S_1}{S_0} + K_l S_1 + K_i \big[\xi (S_1 + 1) - 1 \big],
+```
+yields:
+```math
+S_1 = \frac{(\alpha w - K_i (\xi - 1)) S_0}{\alpha w + (K_l + K_i \xi) S_0}.
+```
+
+The final expression for the modified supersaturation becomes:
+```math
+S_{\text{max}}^{\text{modified}} = \frac{(\alpha w - K_i (\xi - 1)) S_{\text{max}}^{\text{ARG}}}{\alpha w + (K_l + K_i \xi) S_{\text{max}}^{\text{ARG}}}.
+```
+
+- When ``K_l = K_i = 0``, the formula reduces to the original ARG solution.
+- Larger values of $K_l$ or $K_i$ reduce the supersaturation due to vapor uptake by preexisting hydrometeors.
+- This extension allows the ARG activation scheme to be applied at any vertical level, making it suitable for use as a local operator in models that resolve cloud-scale dynamics.
+
+### Example
+
+To illustrate the behavior of the modified ARG scheme, we compare the maximum supersaturation it predicts against that from a detailed parcel model in four idealized scenarios:
+
+- Varying the liquid droplet number concentration with no ice particles present,
+- Varying the ice crystal number concentration with no liquid droplets present,
+- Varying the initial droplet radius with no ice particles present,
+- Varying the initial ice particle radius with no liquid droplets present.
+
+In all cases, the parcel model assumes an initial zero supersaturation and includes activation, condensational growth, and depositional growth during adiabatic ascent. A background population of sulfate aerosol with an initial concentration of 500 cm⁻³ is present in all scenarios. For simplicity, ice particles are assumed to be monodisperse. Liquid droplets formed through activation and those present initially are represented as two separate monodisperse populations.
+
+The figure below shows the ratio of the maximum supersaturation reached in the parcel model to that predicted by the modified ARG scheme. This comparison highlights the conditions under which the modified ARG provides a good approximation to the more detailed parcel-model results.
+
+```@example
+include("plots/ARGmodified_plots.jl")
+```
+![](parcel_vs_modifiedARG_aerosol_activation.svg)
+

docs/src/IceNucleationParcel0D.md

@@ -122,9 +122,7 @@ The crux of the problem is modeling the ``\frac{dq_l}{dt}`` and ``\frac{dq_i}{dt
   for different homogeneous and heterogeneous ice nucleation paths.
 
 ## Supported size distributions
-Currently the parcel model supports monodisperse and gamma size distributions of cloud droplets and ice crystals,
-  and solves prognostic equations for the cloud water and cloud ice specific contents
-  ``q_l``, ``q_i`` and number concentrations ``N_l``, ``N_i``.
+Currently, the parcel model supports monodisperse and gamma size distributions of cloud droplets and ice crystals, and solves prognostic equations for the cloud water and cloud ice specific contents (`q_l`, `q_i`) and number concentrations (`N_l`, `N_i`). Additionally, a `monodisperseMix` option is now supported for cloud liquid droplets. This allows representing the droplet population as a mixture of two monodisperse modes: one corresponding to an initial set of preexisting droplets, and the other representing droplets formed through activation during the simulation. This feature enables more realistic treatment of scenarios involving preexisting hydrometeors.
 
 For a monodisperse size distribution of cloud droplets or ice crystals
 ```math
@@ -167,15 +165,10 @@ As a result
 ## Supported source terms
 ### Aerosol Activation
 Aerosol activation is described by ([see discussion](https://clima.github.io/CloudMicrophysics.jl/dev/AerosolActivation/#Number-and-mass-of-activated-particles)).
-  It is inherently assumed that the aerosols have a lognormal size distribution.
-  For simplicity, the parcel accepts one mode and one aerosol type at a time, therefore,
-  internal mixing is not needed. The maxiumum supersaturation as described in the above
-  mentioned documentation is replaced by the liquid supersaturation in the parcel as
-  it evolves over time.
+It is inherently assumed that the aerosols have a lognormal size distribution. For simplicity, the parcel accepts one mode and one aerosol type at a time, therefore, internal mixing is not needed. The maxiumum supersaturation as described in the above mentioned documentation is replaced by the liquid supersaturation in the parcel as it evolves over time.
 
 !!! note
-    Standard deviation and r_{mean} of the aerosol size distribution may change as aerosols activate.
-    For now, we will neglect these effects.
+    The standard deviation and $r_{\text{mean}}$ of the aerosol size distribution represent the initial (pre-activation) aerosol population and are assumed to remain fixed throughout the simulation. They do not evolve as aerosols activate.
 
 ### Condensation growth
 

docs/src/plots/ARGmodified_plots.jl

@@ -0,0 +1,137 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+
+import ClimaParams as CP
+import CloudMicrophysics as CM
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.ThermodynamicsInterface as TDI
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+
+function run_parcel_model(Nₐ, Nₗ, Nᵢ, rₗ, rᵢ, w, FT)
+
+    # Get free parameters
+    tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
+    wps = CMP.WaterProperties(FT)
+
+    # Initial conditions
+    T₀ = FT(230)
+    cᵥ₀ = FT(5 * 1e-5)
+    ln_INPC = FT(0)
+
+    # Constants
+    ρₗ = wps.ρw
+    ρᵢ = wps.ρi
+    ϵₘ = TDI.Rd_over_Rv(tps)
+    eₛ = TDI.saturation_vapor_pressure_over_liquid(tps, T₀)
+    qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀)
+    Sₗ = FT(1.0)
+    e = Sₗ * eₛ
+    p₀ = e / cᵥ₀
+    ρ_air = TDI.air_density(tps, T₀, p₀, qᵥ, FT(0), FT(0))
+    qₗ = Nₗ * FT(4 / 3 * π) * (rₗ)^3 * ρₗ / ρ_air
+    qᵢ = Nᵢ * FT(4 / 3 * π) * (rᵢ)^3 * ρᵢ / ρ_air
+    IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, ln_INPC, qₗ, Nₗ]
+
+    # Simulation parameters passed into ODE solver
+    const_dt = FT(1e-2)    # model timestep
+    t_max = FT(60)     # total time
+    aerosol = CMP.Sulfate(FT)
+
+    condensation_growth = "Condensation"
+    deposition_growth = "Deposition"
+    aerosol_act = "AeroAct"     # turn on aerosol activation
+    aero_σ_g = FT(2.3)
+    r_nuc = FT(4e-8)
+
+    params = parcel_params{FT}(
+        w = w,
+        const_dt = const_dt,
+        aerosol_act = aerosol_act,
+        aerosol = aerosol,
+        aero_σ_g = aero_σ_g,
+        r_nuc = r_nuc,
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        Nₐ = Nₐ,
+        liq_size_distribution = "MonodisperseMix",
+    )
+
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, params)
+
+    S_max = maximum(sol[1, :]) - FT(1)
+
+    # ARG results
+    ad = AM.Mode_κ(
+        params.r_nuc,
+        params.aero_σ_g,
+        params.Nₐ,
+        (FT(1.0),),
+        (FT(1.0),),
+        (params.aerosol.M,),
+        (params.aerosol.κ,),
+    )
+    all_ad = AM.AerosolDistribution((ad,))
+    S_max_ARG = AA.max_supersaturation(params.aap, all_ad, params.aps, params.tps, T₀, p₀, w, qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+    S_max_mod =
+        AA.max_supersaturation(params.aap, all_ad, params.aps, params.tps, T₀, p₀, w, qᵥ + qₗ + qᵢ, qₗ, qᵢ, Nₗ, Nᵢ)
+
+    return S_max / S_max_ARG, S_max_mod / S_max_ARG
+end
+
+
+FT = Float32
+Nₐ = FT(5e8)
+Nₗ = FT(1e8)
+Nᵢ = FT(1e6)
+rₗ = FT(20e-6)
+rᵢ = FT(20e-6)
+w = FT(1.2)   # updraft speed
+
+n_points = 10
+Nₗ_range = range(FT(0), stop = FT(1e8), length = n_points)
+Nᵢ_range = range(FT(0), stop = FT(1e6), length = n_points)
+rₗ_range = range(FT(0e-6), stop = FT(25e-6), length = n_points)
+rᵢ_range = range(FT(0e-6), stop = FT(25e-6), length = n_points)
+
+S_max_parcel_Nₗ = zeros(n_points)
+S_max_ARGmod_Nₗ = zeros(n_points)
+S_max_parcel_Nᵢ = zeros(n_points)
+S_max_ARGmod_Nᵢ = zeros(n_points)
+S_max_parcel_rₗ = zeros(n_points)
+S_max_ARGmod_rₗ = zeros(n_points)
+S_max_parcel_rᵢ = zeros(n_points)
+S_max_ARGmod_rᵢ = zeros(n_points)
+for i in 1:n_points
+    S_max_parcel_Nₗ[i], S_max_ARGmod_Nₗ[i] = run_parcel_model.(Nₐ, Nₗ_range[i], FT(0), rₗ, rᵢ, w, FT)
+    S_max_parcel_Nᵢ[i], S_max_ARGmod_Nᵢ[i] = run_parcel_model.(Nₐ, FT(0), Nᵢ_range[i], rₗ, rᵢ, w, FT)
+    S_max_parcel_rₗ[i], S_max_ARGmod_rₗ[i] = run_parcel_model.(Nₐ, Nₗ, FT(0), rₗ_range[i], rᵢ, w, FT)
+    S_max_parcel_rᵢ[i], S_max_ARGmod_rᵢ[i] = run_parcel_model.(Nₐ, FT(0), Nᵢ, rₗ, rᵢ_range[i], w, FT)
+end
+
+# Plot results
+fig = MK.Figure(size = (800, 600), fontsize = 20)
+ax1 = MK.Axis(fig[1, 1], ylabel = "S_max / S_max_ARG", xlabel = "Nₗ [cm⁻³]", title = "No ice particle, rₗ=20 μm")
+ax2 = MK.Axis(fig[1, 2], ylabel = "S_max / S_max_ARG", xlabel = "Nᵢ [cm⁻³]", title = "No liquid droplets, rᵢ=20 μm")
+ax3 = MK.Axis(fig[2, 1], ylabel = "S_max / S_max_ARG", xlabel = "rₗ [μm]", title = "No ice particle, Nₗ=100 cm⁻³")
+ax4 = MK.Axis(fig[2, 2], ylabel = "S_max / S_max_ARG", xlabel = "rᵢ [μm]", title = "No liquid droplets, Nᵢ=1 cm⁻³")
+
+MK.lines!(ax1, Nₗ_range * 1e-6, S_max_parcel_Nₗ, label = "Parcel", linewidth = 2, color = :blue)
+MK.lines!(ax1, Nₗ_range * 1e-6, S_max_ARGmod_Nₗ, label = "Modified ARG", linewidth = 2, color = :red)
+MK.lines!(ax2, Nᵢ_range * 1e-6, S_max_parcel_Nᵢ, label = "Parcel", linewidth = 2, color = :blue)
+MK.lines!(ax2, Nᵢ_range * 1e-6, S_max_ARGmod_Nᵢ, label = "Modified ARG", linewidth = 2, color = :red)
+MK.lines!(ax3, rₗ_range * 1e6, S_max_parcel_rₗ, label = "Parcel", linewidth = 2, color = :blue)
+MK.lines!(ax3, rₗ_range * 1e6, S_max_ARGmod_rₗ, label = "Modified ARG", linewidth = 2, color = :red)
+MK.lines!(ax4, rᵢ_range * 1e6, S_max_parcel_rᵢ, label = "Parcel", linewidth = 2, color = :blue)
+MK.lines!(ax4, rᵢ_range * 1e6, S_max_ARGmod_rᵢ, label = "Modified ARG", linewidth = 2, color = :red)
+
+MK.ylims!(ax1, -0.05, 1.1)
+MK.ylims!(ax3, -0.05, 1.1)
+MK.axislegend(ax1, framevisible = false, labelsize = 16, position = :rc)
+MK.axislegend(ax2, framevisible = false, labelsize = 16, position = :rc)
+MK.axislegend(ax3, framevisible = false, labelsize = 16, position = :rc)
+MK.axislegend(ax4, framevisible = false, labelsize = 16, position = :rc)
+
+MK.save("parcel_vs_modifiedARG_aerosol_activation.svg", fig)

papers/ice_nucleation_2024/calibration_setup.jl

@@ -135,7 +135,7 @@ function perf_model_IC(FT, IN_mode)
         e = eᵥ(qᵥ, p₀, Rₐ, R_v)
         Sₗ = FT(e / eₛ)
     end
-    return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
+    return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0), FT(0)]
 end
 
 function perf_model_pseudo_data(FT, IN_mode, params, IC)

parcel/Example_AerosolActivation.jl

@@ -58,6 +58,7 @@ params = parcel_params{FT}(
     aero_σ_g = aero_σ_g,
     r_nuc = r_nuc,
     condensation_growth = condensation_growth,
+    Nₐ = Nₐ,
 )
 
 # solve ODE

parcel/ParcelCommon.jl

@@ -21,3 +21,13 @@ function AIDA_rate(model_t, data_t, data)
         return 0
     end
 end
+
+# Return activated particle radius. See below eq. 19 in [Abdul-Razzaketal1998](@cite)
+function get_particle_activation_radius(
+    ap::CMP.AerosolActivationParameters,
+    T::FT,
+    S::FT,
+) where {FT}
+    A::FT = AA.coeff_of_curvature(ap, T)
+    return FT(2 / 3) * A / S
+end