Tier I — Gaussian family

Forward model: closed-form analytical solution for steady/transient point-source dispersion.

This is the first tier to build to completion (all six steps), because:

1. The forward model is microseconds, so MCMC is feasible end-to-end.
2. Every downstream tier validates against Tier I in the limit of weak turbulence and stationary winds — Tier I is the analytical reference.
3. The amortized predictor at Step 5 doubles as the working prototype for the inference UX (input/output shapes, posterior visualisation, uncertainty budget).

(1) Simple model¶

Gaussian plume (steady-state, continuous source)¶

The primed coordinates $(x', y')$ are the source-aligned frame: $x'$ is downwind along $\theta_\text{wind}$, $y'$ is crosswind. The image-source sum runs over $n \in \{\dots,-1, 0, 1, \dots\}$ to enforce no-flux at the ground ($z=0$) and at the capping inversion ($z=L = \text{PBL height}$). For $L \to \infty$ only the $n=0$ pair survives — that’s the unbounded-domain case in most textbooks. For typical PBL $L \approx 1\,\text{km}$ and $\sigma_z$ growing past $\sim L/3$, the upper image is non-negligible.

Parameters:

• $Q$ — source strength (kg/s)
• $x_0 = (\text{lat}_0, \text{lon}_0)$ — source location
• $H_\text{stack}$ — physical stack height
• $\Delta h(F_b, F_m, \bar{u}, \text{stability})$ — Briggs plume rise from buoyancy/momentum fluxes (Briggs, 1973); $H_\text{eff} = H_\text{stack} + \Delta h$
• $\bar{u}$, $\theta_\text{wind}$ — wind speed and direction (from met, with uncertainty — see Inference)
• $\sigma_y(x'), \sigma_z(x')$ — crosswind/vertical spread; PG-class lookup (Pasquill, 1961Turner, 1970) for v1, MO-similarity (Monin & Obukhov, 1954Stull, 1988) for v2 (see prereqs)
• $L$ — PBL capping height (from met)

For methane super-emitter inversions where stack height is poorly known, $H_\text{eff}$ is part of what you infer; for known facilities $H_\text{stack}$ is fixed and only $\Delta h$ is computed.

Gaussian puff (time-varying, episodic source)¶

Each puff $k$ advects with the wind and diffuses independently — handles intermittent / burst-mode releases that violate the steady-state assumption.

From $c(x,y,z)$ to a satellite-comparable observation¶

The forward used by inference is not $c(x,y,z)$ directly. It’s:

• Column integrate: the satellite sees a column-integrated enhancement, not a single altitude.
• Background $c_\text{bg}(x,y)$: the regional background (typically ~1900 ppb for CH₄). Not optional — without it the model predicts the enhancement above background, which is what $c$ actually is, but the satellite delivers absolute column densities. Either subtract $c_\text{bg}$ from $y_\text{obs}$ upstream, or jointly model it. Either way, it’s day-1 infrastructure, not a Step-6 upgrade.
• Averaging kernel $\mathbf{A}$: the satellite-product AK from the prereqs; see ((1)). Required when comparing to L2 XCH₄ products. Skip it only when working with a flat-AK assumption (rare and worth flagging).

Extended Gaussian (AERMOD-style)¶

Adds terrain corrections (plume rise over hills, beyond Briggs), building downwash (Huber–Snyder), and receptor-grid output. See Cimorelli et al., 2005. Useful as a sanity benchmark against operational tools, not a primary research target.

(2) Model-based inference¶

The plume is analytical, so the full forward $H : (Q, x_0, H_\text{eff}, \bar{u}, \theta_\text{wind}, \dots) \to y_\text{model}(x,y)$ is differentiable via JAX end-to-end (column integration + AK included).

Likelihood model¶

• Default: heteroscedastic Gaussian on column-XCH₄ enhancement, with per-pixel σ from the L2 retrieval-error map.
• Heavy-tail variant: Student-$t$ with $\nu \approx 5$ for retrievals near the detection floor where outliers dominate (TROPOMI single-pass Veefkind et al., 2012, EMIT off-axis Green & others, 2022).
• Mask: quality-flag mask from the L2 product is passed through as an indicator weight on the likelihood — flagged pixels contribute zero log-likelihood.

Priors¶

MAP / MCMC¶

• MAP estimation: jax.grad(log_posterior) → L-BFGS or Adam. Convergence in <1 s for a single overpass.
• MCMC: NumPyro NUTS over $(Q, x_0, H_\text{eff}, \bar{u}, \theta_\text{wind}, c_\text{bg})$ jointly. Forward pass is ~µs, so 10k samples is seconds. Keep $\bar{u}$, $\theta_\text{wind}$ as inference variables (with tight met priors) so their posterior contracts can be diagnosed downstream.
• Linear-Gaussian special case: if observations are linear in $c$ (column-integrated XCH₄ with a fixed AK) and only $Q$ is free with all geometry fixed, the posterior over $Q$ is analytically tractable — exact Bayesian inversion via gaussx. Useful as a sanity check on the NUTS implementation.

$Q / \bar{u}$ identifiability — corrected¶

The classic statement is “$Q$ and $\bar{u}$ enter only as $Q/\bar{u}$, so they’re degenerate in a single transect” Varon et al., 2018. This is true only if $\bar{u}$ has a flat prior. In production, the tight met-derived prior on $\bar{u}$ resolves the degeneracy: the posterior $p(Q \mid y, \bar{u}_\text{met}, \sigma_{\bar{u},\text{met}})$ contracts as long as $\sigma_{\bar{u},\text{met}} / \bar{u}_\text{met} \ll \sigma_Q / Q_\text{prior}$. The “need two crosswind transects” claim is a special case for fully-free $\bar{u}$ and should not be the default reading.

When the wind prior is not tight (rare — typically remote regions without good reanalysis coverage), then yes, multiple overpasses with different wind directions break the ratio. But the first-line fix is the prior, not the geometry.

This is the first working end-to-end inverse pipeline. Build it here, validate against synthetic releases with known truth, then move to real controlled-release data (see Validation).

(3) Model emulator¶

• Input: $(Q, x_0, H_\text{eff}, \bar{u}, \theta_\text{wind}, \text{stability\_class}, L)$ — the raw inputs, not the evaluated σ profiles. The emulator should learn the σ functional too; otherwise you’re benchmarking interpolation, not modelling.
• Output: $y_\text{model}(x,y)$ on the satellite-pixel grid (post column-integration, post AK).
• Architecture: MLP for low-dim parameters → small transposed-conv decoder for spatial output, or lightweight DeepONet.
• Value: validates emulator training pipeline (data generation, training loop, residual checks) before applying to expensive PDE models.

(4) Emulator-based inference¶

Same inference loop as Step 2, but the forward pass is the neural network. Required validation:

• Posterior mean from emulator-MCMC ≈ posterior mean from analytical-MCMC, within $1\sigma$.
• Posterior covariance from emulator-MCMC ≈ analytical posterior covariance, in operator norm.

If those pass for Tier I, the same diagnostics apply unchanged at Tier III.

(5) Amortized inference (predictor)¶

Train a summary network + posterior network mapping observation patches directly to the source posterior:

Input shape — commit to a 2D patch¶

$y_\text{patch}$ is a fixed-size 2D image patch of column XCH₄ enhancement (background-subtracted) centred on the facility candidate, plus per-pixel uncertainty and quality-mask channels. Concretely:

• Shape: $(C, H_\text{patch}, W_\text{patch})$ with $C = 3$ channels (enhancement, $\sigma_\text{retr}$, mask).
• Resolution: instrument-native (e.g. 30 m for Tanager Carbon Mapper, 2024, 60 m for EMIT Green & others, 2022, ~5 km for TROPOMI Veefkind et al., 2012). Train one predictor per instrument family — don’t try to share across resolutions.
• Footprint: large enough to capture the plume’s downwind extent at expected wind speeds (typical: 2–10 km square for point-source instruments).

This commitment matters because it pins the architecture (CNN-style backbones over set-transformers or graph nets).

Conditioning on context¶

The predictor is not observation-only. At inference time we know $\text{context} = (\text{facility}_\text{lat,lon}, \bar{u}_\text{met}, \theta_\text{wind,met}, \text{stability\_class}, L_\text{met})$. The clean way to wire this in is feature-wise modulation of the summary network — the FiLM / hypernet conditioning primitives being built in pyrox.nn.

summary   = CNN(y_patch)              # observation features
modulated = FiLM(summary, context)    # per-feature γ(context)·summary + β(context)
posterior = NPE_head(modulated)       # conditional flow over (Q, x₀, H_eff)

Architecture options¶

• Conditional normalizing flow (gauss_flows) for the posterior head — preferred for low-dim posteriors over $(Q, x_0, H_\text{eff})$.
• BNN posterior head (pyrox) when posterior multimodality is unlikely.
• NPE / SBI-style — drop-in for either, gives a clean SBC validation interface.

Training dataset is free: simulate millions of plume configurations in seconds, sampling met context from the WRF/ERA5 climatology. Validation: posterior calibration via simulation-based calibration (SBC) — uniform rank statistics across 1k held-out simulations, stratified by met regime.

(6) Improve¶

• PG → MO σ swap. Replace PG lookup tables for $\sigma_y, \sigma_z$ with MO-similarity-derived functions parameterised by $(u_*, L_\text{Obukhov}, z_0)$ from the prereqs. Closer to physics, continuous in stability, and the right pre-step for tier-II validation.
• Multi-source. Real basins host 5–50 simultaneous emitters Frankenberg et al., 2016Jacob et al., 2022. This is not an array-shape change — the inference becomes a mixture model with unknown component count, requiring reversible-jump MCMC (RJMCMC) or a Dirichlet-process prior over source count. Plan for this from the start: the posterior interface that Tier V.A consumes must handle variable-$K$ per overpass.
• Learned σ from LES. Train a small NN against LES output to learn stability-dependent σ functions that go beyond MO. Slot in as a swap-in for either PG or MO σ. Useful for super-emitter regimes that LES has resolved better than empirical fits.
• Distributed-source field $Q(\mathbf{x})$. Replace point source with a spatial source field — opens the door to Tier II/III for spatially extended emissions. At Tier I this is just a sum of point sources sharing met context; the multi-source mixture (above) is the same code path with finer support.

Module layout¶

Validation strategy¶

• Forward model — mass flux conservation. Integrate mass flux $\int\!\!\int c \, u_\perp \, \mathrm{d}A$ through a transverse plane downwind of the source, compare to $Q$. Should be exact in the no-deposition, unbounded-domain limit. (Common bug: integrating $\int\!\!\int c \, \mathrm{d}A$ instead — that’s mass per unit wind speed, dimensionally wrong.)
• Ground reflection consistency. Set $H = 0$; the two image-source terms must collapse to a single Gaussian with doubled prefactor. Catches sign/index bugs in the image sum.
• PBL capping. With $L \to \infty$ the upper-image series must vanish; with finite $L$ and $\sigma_z$ growing past $L/2$, concentrations must converge to a vertically well-mixed limit:

Tests both branches of the $n$-sum.

• MAP recovery. Synthetic release with known $(Q^*, x_0^*, H^*)$. Recovery target: $|Q_\text{MAP} - Q^*| < \sigma_\text{post}(Q)$ — i.e. recovery within the posterior’s own claimed uncertainty. (Not a fixed % — that conflates SNR with method quality.)
• MCMC calibration. SBC on 1000 simulated releases — rank histograms uniform across all parameters, stratified by met regime.
• $Q/\bar{u}$ identifiability — empirical. Run MAP over a sweep of $\sigma_{\bar{u},\text{met}} / \bar{u}_\text{met} \in [0.05, 0.5]$. Posterior $\operatorname{CV}(Q)$ should grow monotonically with met-wind uncertainty; flat CV signals a wiring bug.
• Real-data benchmark. Invert a published controlled-release flight (e.g. Stanford / Sherwin et al. 2024 controlled releases over Tanager / GHGSat / EMIT). Posterior median for $Q$ must contain the metered release rate within 95% credible interval. This is the only test that proves the inference works on real radiances — synthetic validation is necessary but not sufficient.
• Emulator agreement. Step 4 posterior matches Step 2 posterior in mean and covariance to within Monte Carlo error.

Aggregating across overpasses¶

The MAP / MCMC posterior $p(Q \mid \text{overpass})$ produced here is the per-event evidence consumed by the population layer. See Tier V.A — Instantaneous emission estimation for the formal interface that turns this posterior into a mark likelihood for the TMTPP fit, and Tier V.D — Total emission estimation for why per-overpass averages systematically misrepresent regional totals.

The multi-source extension (Step 6) makes this interface variable-$K$ per overpass — the V.A adapter must handle that.

Open questions¶