Dry Deposition

Model

DryDepositionGasFractional calculates dry deposition velocities and first-order loss rates for gas-phase species. It uses fractional (mosaic) land use: rather than one land-use class per grid cell, the deposition velocity is the area-weighted sum over the 11 Wesely classes,

Vd = Σᵢ fᵢ / (Ra + Rb + Rcᵢ)

where the aerodynamic resistance Ra and quasi-laminar sublayer resistance Rb follow Seinfeld and Pandis (2006) eqs. 19.13/19.14 and 19.17 and are the same for every class, and the surface resistance Rcᵢ of class i follows Wesely (1989).

Deposition is applied only in the lowest model level; the rate is k = Vd · g · ρA / ΔP.

Building the model

using AtmosphericDeposition
using ModelingToolkit
using ModelingToolkit: t

model = DryDepositionGasFractional()
sys = mtkcompile(model)

The system is purely algebraic — it produces deposition velocities and rates, not concentrations. After mtkcompile there are no unknowns. The observed quantities are a velocity v_X(t) [m/s] and a rate k_X(t) [1/s] for each of the 132 supported species. They also include the quantities that do not depend on the species, which are computed once and shared by every species: the aerodynamic resistance Ra(t), and the parts of the Wesely surface resistance that are the same for all species (the stomatal resistance and table entries of each land-use class, e.g. rs_urban(t) and rac_water(t), and the class-independent rdc(t) and rcold(t)).

observed(sys)   # Ra(t), rdc(t), ..., v_ACET(t), ..., k_RCOOH(t)

Parameters

parameters(sys)

Meteorology and surface state:

parameterdescription
Tssurface air temperature [K]
zheight to the mid-point of level 1 [m]
del_Ppressure thickness of level 1 [Pa]
z₀roughness length [m]
u_starfriction velocity [m/s]
LMonin-Obukhov length [m]
ρAair density [kg/m³]
Gsolar irradiation [W/m²]
θslope of the local terrain [radians]
levlevel index; deposition is applied only where lev == 1
seasonWesely (1989) season index, 1–5

Land-use area fractions, one per Wesely (1989) class, expected to sum to ≈ 1:

f_urban, f_agricultural, f_range, f_deciduous, f_coniferous, f_mixedforest, f_water, f_barren, f_wetland, f_rangeag, f_rockyshrubs.

By default f_mixedforest = 1.0 and the rest are zero.

The remaining parameters (g, κ, Rc_unit, unit_m, …) are physical constants and unit-handling scalars; they are not intended to be varied.

Coupling to chemistry

On its own the model only reports v_X and k_X. To see concentrations change, couple it to a chemical mechanism — the k_X are then applied as loss terms to the matching species:

using AtmosphericDeposition, GasChem, EarthSciMLBase, ModelingToolkit

model = couple(SuperFast(), DryDepositionGasFractional())
sys = convert(System, model)

GEOSChemGasPhase and Pollu are wired the same way; see ext/GasChemExt.jl for the species maps.

Coupling to GEOS-FP meteorology

Coupling with EarthSciData.GEOSFP binds the surface meteorology, drives season from the simulation date, and looks the 11 land-use fractions up from the bundled MODIS-derived CONUS dataset at each cell's longitude and latitude:

using AtmosphericDeposition, EarthSciData, EarthSciMLBase, ModelingToolkit, Dates

domain = DomainInfo(
    DateTime(2022, 5, 1), DateTime(2022, 5, 2);
    latrange = deg2rad(24.0f0):deg2rad(2):deg2rad(50.0f0),
    lonrange = deg2rad(-125.0f0):deg2rad(2.5):deg2rad(-66.0f0),
    levrange = 1:10,
)

model = couple(GEOSFP("4x5", domain), DryDepositionGasFractional(domain))
sys = convert(System, model)

Passing domain to DryDepositionGasFractional declares the simulation grid, and the land-use source is conservatively area-averaged onto it. Without it the system still builds, but the first land-use lookup raises.

Note the bundled land-use dataset covers CONUS only (-130…-60 °E, 20…55 °N); a domain extending outside it raises rather than being clamped to the boundary. season assumes Northern-Hemisphere mid-latitudes.