Emission rate from first principles

The advection–diffusion equation¶

The concentration $c(\vect{s}, t)$ of a trace gas obeys

where $\vect{s} = (x, y, z) \in \mathbb{R}^3$ is the position (along-wind, crosswind and height), $\vect{u} \in \mathbb{R}^3$ is the wind [m/s], $\mat{K} \in \mathbb{R}^{3 \times 3}$ is the turbulent diffusivity [m²/s], $\lambda$ is a first-order loss rate [1/s] and $S$ is the source. Every other symbol on this page is a scalar in $\mathbb{R}$, apart from the column flux $\vect{F}$ and the normal $\vect{n}$ of the divergence method. For incompressible flow ($\nabla \cdot \vect{u} = 0$), the advection term is $\vect{u} \cdot \nabla c$.

Methane. Methane has a lifetime of about nine years, so $\lambda \approx 0$ on a plume’s timescale of minutes to hours. We align $x$ with a mean wind $(U, 0, 0)$, keep a diagonal $\mat{K}$ with entries $K_x, K_y, K_z$, and drop the along-wind diffusion when the Péclet number $\mathrm{Pe} = U L / K_x$ is large:

The column. A satellite measures the column $\Omega(x, y, t) = \int_0^\infty c \mathrm{d}z$. With no flux through the ground and $c \to 0$ aloft, the vertical diffusion integrates to zero, and for a point source $q = Q \delta(x) \delta(y)$ we get

Pulling $U$ out of the $z$ integral assumes that the wind does not vary with height, and stage 2 measures what that costs. Each method below comes either from (2) or from this column equation (Figure 1).

The Gaussian plume (plumax)¶

TL;DR. Solve (2) at steady state, integrate it over height, and invert the column for $Q$.

Problem. At steady state, $x$ plays the role of time (the travel time $\tau_x = x / U$), and (2) is a two-dimensional heat equation in $(y, z)$, $\partial_{\tau_x} c = K_y \partial_{yy} c + K_z \partial_{zz} c$. Integrating (2) across $x = 0$ turns the source into its initial condition, $U c(0^+, y, z) = Q \delta(y) \delta(z - H)$, so the source injects $Q / U$ per metre of travel. With a reflecting ground and a release height $H$,

where $\sigma_y(x)$ and $\sigma_z(x)$ are the crosswind and vertical dispersion widths. The column of (4) is

which we invert for the emission rate on the plume’s axis:

1
2
3
4
5
6
7
8
9
10
11
# Shapes: Q = 2 kg/s, U = 3 m/s, class D; x = 500, 1500 m
# σ = a x (1 + b x)^c, Briggs, for y and z
D = get_dispersion_params("D")                         # today
# c(x, y, z) by (er1-eq-plume3d), with the image source
c = plume_concentration(x, y, z, 0, 0, H, U, 0, Q, D)  # today
# Ω = ∫₀^∞ c dz, quadrature, H = 0 and 30 m          (·) → ()
# → Q / (√(2π) U σ_y) e^{−y²/2σ_y²} to 10⁻⁶         user
# σ_y: 39.0 m at 500 m, 146.1 m at 2 km
# d ln σ_y / d ln x = 0.99 between 100 and 400 m
# U ∬ c dy dz on a 1 m grid, H = 30 m              (Y, Z) → ()
# → Q to 10⁻³ at 500 m and 1.5 km                   user

The column formula (5) holds for any release height. The two image terms together integrate to $\int_{-\infty}^{\infty} e^{-(z-H)^2/2\sigma_z^2} \mathrm{d}z = \sqrt{2\pi} \sigma_z$ whatever $H$ is, so plumax’s plume released at 30 m has the same column as one released at the ground, to quadrature accuracy. In addition, the flux through every cross-section, $U \iint c \mathrm{d}y \mathrm{d}z$, is $Q$, which is simply the initial condition carried downwind (and the cross-sectional flux in three dimensions).

The Gaussian puff (plumax)¶

TL;DR. Release mass in puffs, keep along-wind diffusion, and recover the plume by superposition.

Problem. An instantaneous release of mass $M$ at $(0, 0, H)$ and $t = 0$ is the initial condition $c(\vect{s}, 0) = M \delta(x) \delta(y) \delta(z - H)$ of (1), with $S = 0$ afterwards. In the frame moving with the wind, $x' = x - U t$, the chain rule cancels the advection term and leaves the heat equation $\partial_t c = K_x \partial_{x'x'} c + K_y \partial_{yy} c + K_z \partial_{zz} c$. Its Green’s function is the product of three one-dimensional heat kernels, $(4\pi K_i t)^{-1/2} e^{-x_i^2 / 4 K_i t}$ (with $x_i$ the $i$-th coordinate of the moving frame), and with $\sigma_i^2 = 2 K_i t$ and an image at $z = -H$ for a reflecting ground, we get

which is plumax’s puff_concentration, and its column is

again for any release height. A steady release is then the superposition $\Omega = \int_0^t Q \Omega_1(x, y, t - t') \mathrm{d}t'$ of unit puffs, where each $\Omega_1$ is (8) with $M = 1$.

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: release ages a ∈ (0, 3600 s], 20,000 of them
# Ω(x, 0) = ∫₀^T Q Ω₁(x, 0; a) da                    (a,) → ()
# against Q / (√(2π) U σ_y), σ_y² = 2 K x / U
superposed(500, u=5) / kplume(500, u=5)  # Pe 250 → 0.999  user
superposed(50, u=5) / kplume(50, u=5)    # Pe 25  → 0.990  user
superposed(500, u=1) / kplume(500, u=1)  # Pe 50  → 0.995  user
superposed(50, u=1) / kplume(50, u=1)    # Pe 5   → 0.958  user
# plumax's puff, σ² = 2 K t, t = 100 s, H = 30 m, 2 m grid
c = puff_concentration(x - U*t, y, z, 0, 0, H, σ, σ, σ, M)   # today
# ∭_{z ≥ 0} c dV = M to 10⁻⁶                  (301, 301, 201) → ()
# ∫ c dz = (er1-eq-puffcol) to 10⁻⁶      (301, 301, 201) → (301, 301)
# ∂c/∂t + U ∂c/∂x − K ∇²c = 0 to 10⁻⁴, central differences   user

The puff sum and the plume agree to 0.1 % at a Péclet number of 250, and they differ by 4 % at 5, where the along-wind diffusion (which the plume drops) matters. So this is the note’s point about low winds made quantitative: the plume formula is the high-Péclet limit of the puff sum, not its exact steady state.

The cross-sectional flux (geotoolz)¶

TL;DR. Integrate the column across the plume at several downwind distances and multiply by the wind.

Problem. We integrate (3) at steady state over all $y$. The diffusion term integrates to $K_y [\partial_y \Omega]_{-\infty}^{\infty} = 0$, because the plume and its gradient vanish far to the sides, and the source term integrates to $Q \delta(x)$:

which, with nothing upwind, integrates to

where $\Theta$ is the step function. So the flux is constant downwind and equals $Q$:

When the wind $u(z)$ varies with height, the flux through a transect is $\iint u(z) c \mathrm{d}y \mathrm{d}z$, so the wind that multiplies the column is the concentration-weighted one,

1
2
3
4
5
6
7
8
9
10
# Shapes: 100 × 160 pixels at 20 m; transects every 400 m
csf = gz.plume.CrossSectionalFlux(
    source=(0, 0), wind_u=3.0, wind_v=0.0,
    n_transects=5, transect_spacing_m=400)             # today
# Q_k = U Σ_j Ω_kj Δy                          (100, 160) → (5,)
csf(column, x > 0)["flux_kg_s"]           # → 2.000 ×5  today
csf(column, column > 0.002)["flux_kg_s"]  # thresholded today
# → 1.90, 1.62, 1.25, 0.77, 0
# u(z) = 3 (z / 10 m)^0.15;  U_eff = ∫ u c dz / ∫ c dz    user
# → 2.55 m/s at 100 m, 3.38 at 1 km, 3.62 at 2 km

On full transects the flux is $Q$ at every distance, as (11) says (Figure 2). However, a detection threshold removes the plume’s flanks, and the far transects lose more: 5 % at 400 m, 61 % at 1.6 km, and everything at 2 km, where the whole transect is below 2 g/m². Shear changes the wind that applies, so a plume hugging the ground moves slower than the 10 m wind (2.55 against 3 m/s at 100 m) whereas a deeper one moves faster (3.62 at 2 km).

The divergence method (xrtoolz budgets)¶

TL;DR. Map the source as the divergence of the column flux, and check Gauss’s theorem on a box.

Problem. We write the column equation in flux form with the column flux $\vect{F} = \int \vect{u}_\text{h} c \mathrm{d}z \approx \vect{u}_\text{h} \Omega$, where $\vect{u}_\text{h}, \vect{F} \in \mathbb{R}^2$ are the horizontal wind and flux, and integrate it over a region $\mathcal{A}$ which holds the source. Gauss’s theorem turns both divergences into boundary integrals, with $\vect{n} \in \mathbb{R}^2$ the outward unit normal and $l$ the arc length along $\partial \mathcal{A}$:

Averaged over a long period, and dropping the diffusive term, the source density is

and the emission rate of the region is

1
2
3
4
5
6
7
8
9
10
# Shapes: F_x = U Ω (100, 160); boxes x ∈ [−100, 1500] m
# q = ∂F_x/∂x, central differences       (100, 160) → (100, 160)
div = np.gradient(U * column, 20.0, axis=1)            # user
# Q_A = ∬_A q dA = ∮ F·n dl                      (100, 160) → ()
area_sum(box(-100, 1500, half=600))  # → 2.000 kg/s    user
area_sum(box(-100, 1500, half=60))   # → 0.81, too narrow
# the dropped term on that box, K_eff = (U/2) dσ_y²/dx (Briggs)
# 2 ∫₀^{x₁} K_eff Ω(x, h) h / σ_y² dx → 1.18; + 0.82 = Q   user
# one overpass, σ_Ω = 0.003 kg/m²: sd = U σ_Ω Δx √n_rows
# → 1.39 kg/s (2,000 draws: 1.39; this draw 1.51); 49 for a 10 % sd

The area sum and the boundary flux agree to 10^-6, as Gauss’s theorem requires, and a box wide enough to hold the plume returns $Q$. However, a box narrower than the plume returns 0.81 kg/s. The rest is the diffusive term of (13), which the method drops: plumax’s Briggs column solves the column equation with an effective diffusivity $K_\text{eff}(x) = \frac{U}{2} \mathrm{d}\sigma_y^2 / \mathrm{d}x$, and with it the diffusive flux through the narrow box’s side edges is 1.18 kg/s. The advective part (0.82 exactly, 0.81 on the grid) plus the diffusive part is $Q$ to 10^-6. So dropping diffusion is only safe for a boundary which crosses the plume where advection carries it (at the downwind edge) and lies outside it everywhere else. With noise, the area sum of a central difference telescopes to four edge columns, each weighted ½, so its uncertainty is theirs: 1.39 kg/s for one overpass at 0.003 kg/m² of noise (70 % of the signal), and 49 overpasses for a 10 % uncertainty. This is why the method averages many overpasses, as the note says it does at TROPOMI’s scale. xrtoolz’s budgets (control_volume_integral, boundary_flux) computes both sides of (15) on a model grid, but nothing in the stack applies them to methane columns yet.

The integrated mass enhancement (geotoolz)¶

TL;DR. Sum the mass above a detection threshold, divide by a residence time, and see what the boundary terms really are.

Problem. Let $\mathcal{A}$ be the region where $\Omega > \Omega_\text{thr}$, for a detection threshold $\Omega_\text{thr}$. On its boundary $\Omega = \Omega_\text{thr}$, so the advective boundary integral is

because the outward normal of a closed curve integrates to zero. Away from the source, the total flux balance of $\mathcal{A}$ then leaves only diffusion, which means everything the source puts in leaves across the contour by crosswind diffusion,

The IME sidesteps the boundary with a residence time, $\tau_\text{res} = \mathrm{IME} / Q$, which is parameterised as $\tau_\text{res} = L / U_\text{eff}$. A scaling argument gives the same form from either balance. Balancing the source against advection, $U \Omega / L \sim Q / A$ with $\mathrm{IME} \sim \Omega A$, gives $Q \sim U \mathrm{IME} / L$. Balancing it against what really leaves (the crosswind diffusion out of a plume of width $W$ and length $L$) gives $Q \sim K_y \mathrm{IME} / W^2$, and the diffusive growth $W^2 \sim K_y L / U$ turns that into the same $U \mathrm{IME} / L$. The order-one factor which the scaling leaves out is what the calibration fixes. So, with $\Omega_j$ and $A_j$ the column and the area of pixel $j$,

and the emission rate is

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: K plume, K_y = 10 m²/s; contour Ω = 4 g/m², 663 m
# ∮ U Ω n_x dl = −2 U Ω_thr [f(x_end) − f(0)]  → 0     user
# −∮ K_y ∂Ω/∂y n_y dl = 2 K_y Ω_thr ∫ f/σ² dx  → 2.000 = Q  user
# IME = Σ Ω_j A_j over Ω > 2 g/m², plumax: 796 kg; τ_res = 398 s
# L, chord of the mask                          (100, 160) → ()
L = plume_length(mask, transform, method="max_axis")   # today
# → 1,740 m; √A = 465 m
# U_eff = Q L / IME: 1.17 m/s (√A), 4.37 m/s (chord)    user
# Q = U L⁻¹ IME, raw wind and chord             (100, 160) → ()
ime = gz.plume.IMEEstimate(wind_speed=3.0,
                           length_method="max_axis")   # today
ime(column, mask)                         # → 1.37 kg/s

Both contour integrals, (16) and (17), are exact: the advective flux through the contour is zero, and the diffusive flux is $Q$ to 10^-9. So the residence-time argument is not a way around small boundary terms. Instead, it replaces the one boundary term that carries the whole source (diffusion across the plume’s edge) with a calibrated wind, which is also why $U_\text{eff}$ depends on the threshold, the noise and the turbulence, as the note says. The calibration is tied to the length scale it was fitted with: the same noise-free plume needs $U_\text{eff} = 1.17$ m/s with $L = \sqrt{A}$ and 4.37 m/s with the chord, a factor of 3.7. geotoolz’s IMEEstimate, which computes (18) and (19), offers three length methods (max_axis by default, convex_hull and skeleton), and none of them is $\sqrt{A}$.

The five methods side by side¶

The divergence method and the cross-sectional flux (Table 1) are the same statement in differential and in integral form. The IME is the area analogue of the cross-sectional flux, which trades the wind direction for a calibration. plumax also carries forward models that the note does not cover (a Lagrangian particle model and an Eulerian large-eddy model on finitevolx), which can generate the synthetic plumes an IME calibration is fitted on. The cross-sectional flux, the divergence method and the IME only use the column equation and not the Gaussian solution, so they hold for any plume shape.

Open questions¶

1. Which length does the stack’s IME calibration use? geotoolz takes $U_\text{eff}$ as an input. The calibration and the length method should travel together as one object, so that a $\sqrt{A}$ coefficient cannot meet a chord length.

2. Where should the divergence method live? It is a budget over a methane column on a grid, so xrtoolz’s budgets compute it, plumax owns the physics, and a TROPOMI pipeline owns the averaging.

3. Can plumax calibrate $U_\text{eff}$ itself? Its large-eddy model can generate the synthetic plumes which Varon et al. made with WRF, and a calibration workflow would make the coefficients reproducible per instrument.

Upcoming work¶

Table 2 lists the follow-ups.

References¶

• The source note, “Emission Rate Estimation from First Principles” Johnson, 2026.

• The IME and its calibration Varon et al., 2018, Sentinel-2 retrievals Varon et al., 2021, the airborne flux distribution Frankenberg & others, 2016, the cross-sectional flux Krings & others, 2011, the divergence method Beirle & others, 2019Liu & others, 2021Koene & others, 2024, the review Jacob & others, 2022 and the plume and puff comparison Jia et al., 2025.

• plumax, geotoolz, georeader and xrtoolz at the revisions in the header.