Air Refreshing and Cloud Ice Uptake Limitations (Luo & Yu, 2023)

Overview

This module implements novel parameterizations for air refreshing limitation and cloud ice uptake limitation on wet scavenging, as described by Luo and Yu (2023). These parameterizations address two key unresolved processes in global atmospheric models:

1. Air refreshing limitation (Section 2.1): Accounts for the finite time needed for species in subgrid cloud-free (rain-free) air to be transported and mixed into cloudy (rainy) air before being subject to wet scavenging. The standard well-mixed assumption overestimates wet scavenging when the subgrid air mixing time scale is comparable to the model time step.

2. Cloud ice uptake limitation (Section 2.2): Limits cold cloud rainout efficiency by the rate at which interstitial water-soluble aerosols are captured by ice crystals via coagulation. The original GEOS-Chem assumption of 100% rainout efficiency for cold clouds is replaced by a physically-based efficiency that depends on ice crystal properties and temperature.

Reference: Luo, G., & Yu, F. (2023). Impact of air refreshing and cloud ice uptake limitations on vertical profiles and wet depositions of nitrate, ammonium, and sulfate. Geophysical Research Letters, 50, e2023GL104258. https://doi.org/10.1029/2023GL104258

AirRefreshingLimitation(; name = :AirRefreshingLimitation)

CloudIceUptakeLimitation(; name = :CloudIceUptakeLimitation)

WetScavengingLimitations(; name = :WetScavengingLimitations)

Implementation

Air Refreshing Limitation

The air refreshing limited removal rate $R_A$ (Eq. 2) is:

\[R_A = \frac{1}{\frac{1}{f R_i} + \tau_A}\]

where $f$ is the cloud/rainfall fraction, $R_i$ is the in-cloud/under-rain removing rate, and $\tau_A$ is the grid refreshing time (Eq. 5):

\[\tau_A = \frac{1 - f}{f K_i}\]

The cloudy air refreshing rate $K_i$ (Eq. 11) depends on turbulence kinetic energy and grid spacing:

\[K_i = \sqrt{\frac{2}{3}\text{TKE}} \left[\frac{1}{\sqrt{f}\,\Delta x} + \frac{1}{\sqrt{f}\,\Delta y} + \frac{1}{\Delta z}\right]\]

Cloud Ice Uptake Limitation

The cold cloud rainout efficiency $F_I$ (Eq. 12) is:

\[F_I = 1 - e^{-R_{A,U}\,\Delta t}\]

where $R_{A,U}$ is the air refreshing limited cloud ice uptake rate (Eq. 13):

\[R_{A,U} = \frac{f\,R_U\,K_i}{K_i + (1-f)\,R_U}\]

The cloud ice uptake rate $R_U$ for HNO₃ (Eq. 14, based on Jacob, 2000) is given in the original CGS form as:

\[R_{U,\text{HNO}_3} = \frac{N_I \cdot S_I}{\frac{r}{D_g} + \frac{2.749064 \times 10^{-4}\,\sqrt{M}}{\gamma\,\sqrt{T}}}\]

where $N_I$ is ice crystal number concentration (cm⁻³), $S_I$ is ice crystal surface area (cm²), $r$ is ice crystal radius (cm), $D_g$ is gas-phase diffusion coefficient (cm² s⁻¹), $M$ is molar mass (g mol⁻¹), $T$ is temperature (K), and $\gamma$ is uptake efficiency. The prefactor $2.749064 \times 10^{-4}$ has units s cm⁻¹ and equals $4\sqrt{\pi/(8R)}$ from kinetic gas theory, where $R = 8.314 \times 10^7$ erg mol⁻¹ K⁻¹. The implementation uses fully SI units with $R = 8.314$ J mol⁻¹ K⁻¹, giving a prefactor of $4\sqrt{\pi/(8R)} \approx 0.8693$ s m⁻¹ when $M$ is in kg mol⁻¹. $M$ and $T$ are non-dimensionalized by reference values of 1 kg mol⁻¹ and 1 K respectively to avoid fractional unit powers.

The uptake efficiency for HNO₃ on ice (Eq. 15, based on Hudson et al., 2002) is:

\[\gamma = 3 \times 10^{-3} + 4 \times 10^{-3} \max\!\left(0,\, \min\!\left(1,\, \frac{220 - T}{220 - 209}\right)\right)\]

State Variables

using DataFrames, ModelingToolkit, DynamicQuantities
using AtmosphericDeposition

sys = WetScavengingLimitations()
vars = unknowns(sys)
DataFrame(
    :Name => [string(ModelingToolkit.Symbolics.tosymbol(v, escape = false)) for v in vars],
    :Units => [dimension(ModelingToolkit.get_unit(v)) for v in vars],
    :Description => [ModelingToolkit.getdescription(v) for v in vars]
)

Parameters

pars = parameters(sys)
DataFrame(
    :Name => [string(ModelingToolkit.Symbolics.tosymbol(p, escape = false)) for p in pars],
    :Units => [dimension(ModelingToolkit.get_unit(p)) for p in pars],
    :Description => [ModelingToolkit.getdescription(p) for p in pars]
)

Equations

eqs = equations(sys)

\[ \begin{align} \mathtt{K_i}\left( t \right) &= \sqrt{\frac{2}{3} \mathtt{TKE}} \left( \frac{1}{\mathtt{{\Delta}z}} + \frac{1}{\sqrt{f} \mathtt{{\Delta}x}} + \frac{1}{\sqrt{f} \mathtt{{\Delta}y}} \right) \\ \mathtt{\tau\_A}\left( t \right) &= \frac{1 - f}{f \mathtt{K_i}\left( t \right)} \\ \mathtt{R\_A}\left( t \right) &= \frac{1}{\mathtt{\tau\_A}\left( t \right) + \frac{1}{\mathtt{R_i} f}} \\ \gamma\left( t \right) &= \mathtt{\gamma\_base} + max\left( \mathtt{zero\_dimless}, min\left( \mathtt{one\_dimless}, \frac{ - T + \mathtt{T\_upper\_luo}}{ - \mathtt{T\_lower\_luo} + \mathtt{T\_upper\_luo}} \right) \right) \mathtt{\gamma\_delta} \\ \mathtt{R\_U}\left( t \right) &= \frac{\mathtt{N\_I} \mathtt{S\_I}}{\frac{\mathtt{r\_ice}}{\mathtt{D\_g}} + \frac{\mathtt{kinetic\_prefactor\_si} \sqrt{\frac{M}{\mathtt{M\_ref}}}}{\sqrt{\frac{T}{\mathtt{T\_ref}}} \gamma\left( t \right)}} \\ \mathtt{R\_AU}\left( t \right) &= \frac{f \mathtt{R\_U}\left( t \right) \mathtt{K_i}\left( t \right)}{\mathtt{K_i}\left( t \right) + \left( 1 - f \right) \mathtt{R\_U}\left( t \right)} \\ \mathtt{F\_I}\left( t \right) &= 1 - e^{ - \mathtt{R\_AU}\left( t \right) \mathtt{{\Delta}t}} \end{align} \]

Analysis

Reported Impacts (Table 1, Luo & Yu, 2023)

The paper reports the following relative changes (%) in concentrations of HNO₃, inorganic nitrate (InOrgNIT), NH₄, and SO₄ at different pressure levels when applying the air refreshing (AR) limitation versus the baseline (BSL), and the combined air refreshing + ice uptake (ARIU) limitation versus BSL:

Key findings:

• Air refreshing limitation has the largest impact below 800 hPa where turbulence in the planetary boundary layer is strong
• Cloud ice uptake limitation is most important above 500 hPa where cold cloud scavenging dominates
• The combined effect (ARIU) shows the largest changes in the upper troposphere (<500 hPa)

These results were obtained using GEOS-Chem v14.1.1 and cannot be reproduced from the parameterization alone. The following analysis plots demonstrate the behavior of the individual parameterization components.

Effect of TKE on Air Refreshing Limited Removal Rate

This figure shows how $R_A / (f \cdot R_i)$ varies with turbulence kinetic energy (TKE) for different cloud fractions. When TKE is large (strong mixing), $R_A \approx f R_i$ (ratio approaches 1). When TKE is small (weak mixing), the air refreshing limitation significantly reduces wet scavenging.

using AtmosphericDeposition: air_refreshing_limited_rate, grid_refreshing_time,
                             cloudy_air_refreshing_rate
using Plots

@parameters f_p Rᵢ_p [unit = u"s^-1"] TKE_p [unit = u"m^2/s^2"]
@parameters Δx_p [unit = u"m"] Δy_p [unit = u"m"] Δz_p [unit = u"m"]

R_A_expr = air_refreshing_limited_rate(f_p, Rᵢ_p,
    grid_refreshing_time(f_p, cloudy_air_refreshing_rate(f_p, TKE_p, Δx_p, Δy_p, Δz_p)))

# GEOS-Chem-like grid: ~200 km horizontal, ~1 km vertical
base = Dict(Rᵢ_p => 0.001, Δx_p => 200000.0, Δy_p => 200000.0, Δz_p => 1000.0)

tke_range = 10 .^ range(-2, 2, length = 100)
cloud_fracs = [0.1, 0.3, 0.5, 0.8]

p = plot(xlabel = "TKE (m²/s²)", ylabel = "R_A / (f·Rᵢ)",
    title = "Air Refreshing Limitation vs. TKE",
    xscale = :log10, ylim = (0, 1.05), legend = :bottomright)

for f_val in cloud_fracs
    ratios = Float64[]
    for tke_val in tke_range
        vals = merge(base, Dict(f_p => f_val, TKE_p => tke_val))
        R_A_val = Float64(ModelingToolkit.Symbolics.value(substitute(R_A_expr, vals)))
        push!(ratios, R_A_val / (f_val * 0.001))
    end
    plot!(p, tke_range, ratios, label = "f = $f_val")
end
p

Effect of Grid Spacing on Air Refreshing Limitation

This figure shows how the air refreshing limitation depends on horizontal grid spacing (Δx = Δy) for a typical TKE value. Coarser grids experience stronger limitation because the subgrid air mixing takes longer.

R_A_expr2 = air_refreshing_limited_rate(f_p, Rᵢ_p,
    grid_refreshing_time(f_p, cloudy_air_refreshing_rate(f_p, TKE_p, Δx_p, Δy_p, Δz_p)))

base2 = Dict(f_p => 0.3, Rᵢ_p => 0.001, Δz_p => 1000.0)
dx_range = 10 .^ range(3, 6, length = 100)  # 1 km to 1000 km
tke_vals_plot = [0.1, 0.5, 1.0, 5.0]

p2 = plot(xlabel = "Horizontal grid spacing (km)", ylabel = "R_A / (f·Rᵢ)",
    title = "Air Refreshing Limitation vs. Grid Spacing (f = 0.3)",
    xscale = :log10, ylim = (0, 1.05), legend = :topright)

for tke_val in tke_vals_plot
    ratios = Float64[]
    for dx_val in dx_range
        vals = merge(base2, Dict(TKE_p => tke_val, Δx_p => dx_val, Δy_p => dx_val))
        R_A_val = Float64(ModelingToolkit.Symbolics.value(substitute(R_A_expr2, vals)))
        push!(ratios, R_A_val / (0.3 * 0.001))
    end
    plot!(p2, dx_range ./ 1000, ratios, label = "TKE = $tke_val m²/s²")
end
p2

HNO₃ Uptake Efficiency on Ice (Eq. 15)

The uptake efficiency $\gamma$ of nitric acid on ice crystals as a function of temperature, based on Hudson et al. (2002). The efficiency transitions from 0.003 at $T \geq 220$ K to 0.007 at $T \leq 209$ K.

using AtmosphericDeposition: hno3_uptake_efficiency, T_upper_luo, T_lower_luo,
                             γ_base, γ_delta, zero_dimless, one_dimless

@parameters T_plot [unit = u"K"]
γ_expr = hno3_uptake_efficiency(T_plot)

const_defaults = Dict(
    T_upper_luo => 220.0, T_lower_luo => 209.0,
    γ_base => 3e-3, γ_delta => 4e-3,
    zero_dimless => 0.0, one_dimless => 1.0
)

T_range = 190.0:0.5:240.0
γ_vals = [Float64(ModelingToolkit.Symbolics.value(
              substitute(γ_expr, merge(const_defaults, Dict(T_plot => T_val)))))
          for T_val in T_range]

p3 = plot(T_range, γ_vals .* 1000,
    xlabel = "Temperature (K)", ylabel = "γ (×10⁻³)",
    title = "HNO₃ Uptake Efficiency on Ice (Eq. 15, Hudson et al., 2002)",
    legend = false, linewidth = 2)
vline!([209.0, 220.0], linestyle = :dash, color = :gray, label = "")
p3

Cold Cloud Rainout Efficiency

The cold cloud rainout efficiency $F_I$ as a function of the air refreshing limited cloud ice uptake rate $R_{A,U}$ for different time steps. Longer time steps and higher uptake rates lead to more complete scavenging.

using AtmosphericDeposition: cold_cloud_rainout_efficiency

@parameters R_AU_plot [unit = u"s^-1"]
@parameters Δt_plot [unit = u"s"]
F_I_expr = cold_cloud_rainout_efficiency(R_AU_plot, Δt_plot)

R_AU_range = 10 .^ range(-5, -1, length = 100)
dt_vals = [60.0, 300.0, 600.0, 1800.0, 3600.0]

p4 = plot(xlabel = "R_A,U (s⁻¹)", ylabel = "F_I",
    title = "Cold Cloud Rainout Efficiency (Eq. 12)",
    xscale = :log10, ylim = (0, 1.05), legend = :topleft)

for dt_val in dt_vals
    F_I_vals = [Float64(ModelingToolkit.Symbolics.value(
                    substitute(F_I_expr, Dict(R_AU_plot => rau, Δt_plot => dt_val))))
                for rau in R_AU_range]
    label_str = dt_val < 3600 ? "Δt = $(Int(dt_val)) s" : "Δt = $(Int(dt_val/3600)) h"
    plot!(p4, R_AU_range, F_I_vals, label = label_str)
end
p4