Emission rate from first principles

The advection–diffusion equation¶

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

with $\mathbf{u}$ the wind [m/s], $\mathbf{K}$ the turbulent diffusivity [m²/s], $\lambda$ a first-order loss rate [1/s] and $S$ the source. For incompressible flow, $\nabla \cdot \mathbf{u} = 0$, the advection term is $\mathbf{u} \cdot \nabla c$.

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

The column. A satellite measures $\Omega(x, y, t) = \int_0^\infty c\, dz$. 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)$:

Pulling $U$ out of the $z$ integral assumes the wind does not vary with height; stage 2 measures what that costs. Each method below comes 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 ($\tau = x / U$) and (2) is a two-dimensional heat equation in $(y, z)$, $\partial_\tau 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)$: the source injects $Q / U$ per metre of travel. With a reflecting ground,

and the column of (4) is

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}\,dz = \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. And the flux through every cross-section, $U \iint c\, dy\, dz$, is $Q$: 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(\mathbf{x}, 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}$, and with $\sigma_i^2 = 2 K_i t$ and an image at $z = -H$ for a reflecting ground,

which is plumax’s puff_concentration. Its column is

again for any release height. A steady release is the superposition $\Omega = \int_0^t Q\, \Omega_1(x, y, t - t')\, dt'$ of unit puffs, each (7) 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 differ by 4 % at 5, where along-wind diffusion, which the plume drops, matters. 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. 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 to $Q\, \delta(x)$:

with $\Theta$ the step function in (8) and nothing upwind. So the flux is constant downwind and equals $Q$:

When the wind varies with height, the flux through a transect is $\iint u(z)\, c\, dy\, dz$, 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 (9) says (Figure 2). 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: a plume hugging the ground moves slower than the 10 m wind (2.55 against 3 m/s at 100 m), and a deeper one 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. Write the column equation in flux form with $\mathbf{F} = \int \mathbf{u}_H\, c\, dz \approx \mathbf{u}_H\, \Omega$ and integrate it over a region $\mathcal{A}$ holding the source. Gauss’s theorem turns both divergences into boundary integrals:

Averaged over a long period and dropping the diffusive term, the source density 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$. A box narrower than the plume returns 0.81 kg/s. The rest is the diffusive term of (11), which the method drops: plumax’s Briggs column solves the column equation with an effective diffusivity $K_{\text{eff}}(x) = \tfrac{U}{2}\, d\sigma_y^2 / dx$, and with it the diffusive flux through the narrow box’s side edges is 1.18 kg/s. Advective (0.82 exactly, 0.81 on the grid) plus diffusive is $Q$ to 10^-6. Dropping diffusion is safe only for a boundary that 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. That 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 (12) on a model grid; 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 > t$. On its boundary $\Omega = t$, so the advective boundary integral is

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

The IME sidesteps the boundary with a residence time, $\tau = \text{IME} / Q$, parameterised as $\tau = 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 $\text{IME} \sim \Omega A$, gives $Q \sim U\, \text{IME} / L$. Balancing it against what really leaves, crosswind diffusion out of a plume of width $W$ and length $L$, gives $Q \sim K_y\, \text{IME} / W^2$, and diffusive growth, $W^2 \sim K_y L / U$, turns that into the same $U\, \text{IME} / L$. The order-one factor the scaling leaves out is what the calibration fixes:

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 t [f(x_end) − f(0)]  → 0         user
# −∮ K_y ∂Ω/∂y n_y dl = 2 K_y t ∫ f/σ² dx  → 2.000 = Q  user
# IME = Σ Ω_j A_j over Ω > 2 g/m², plumax: 796 kg; τ = 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, (13) and (14), 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. It replaces the one boundary term that carries the whole source, diffusion across the plume’s edge, with a calibrated wind. That 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 (15), offers three length methods (max_axis by default, convex_hull, skeleton), none of them $\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 integral form. The IME is the cross-sectional flux’s area analogue: it trades the wind direction for a calibration. plumax also carries forward models 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 use only the column equation, 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 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: 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 Varon et al. made with WRF; 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.