The Beer–Lambert column retrieval

The forward model¶

The reflectance a sensor measures in band $b$ at pixel $(x, y)$ is the surface reflectance times the transmittance of the two-way path through the atmosphere, averaged over the band’s spectral response $f_b(\lambda)$:

Here $\sigma(\lambda)$ is the methane cross-section per unit mass, $\Omega_{\text{bg}}$ the background column and $\Delta\Omega$ the plume’s enhancement. When every line is optically thin, the band average of the exponential is the exponential of the band-averaged cross-section, and (1) becomes the source note’s model:

so that the plume’s optical depth is linear in the column:

(2) is the right picture for small plumes and the wrong one for strong lines. The running example measures how wrong in stage 2.

The variables¶

Table 1 gives each term of (1), its units and the piece of the stack that supplies it.

Where each term comes from¶

Each term has one source in a satellite product (Figure 1).

Reflectance. For Sentinel-2 L1C the granule’s JPEG 2000 bands are top-of-atmosphere reflectance scaled to integers:

with a quantification value of 10,000 and, since processing baseline 04.00 (January 2022), an offset of −1,000. Hyperspectral L1 products are radiance $L_b$, converted with the band’s solar irradiance $E_b$ and the Earth–Sun distance $d$ in astronomical units:

Air-mass factor. The geometric two-way path through a plane-parallel layer is

with $\theta_s$ and $\theta_v$ the sun and view zenith angles.

Cross-section. A line list gives $\sigma(\lambda)$ at the layer’s pressure and temperature; the band’s response function averages it, (2), or, for the exact model, averages the transmittance, (1).

1
2
3
4
5
6
7
8
9
10
11
# Shapes: S2 L1C tile (13, H, W); EMIT L1B (285, H, W); SRF table (877, 13)
s2 = S2ImageL1C(safe_path)                                              # today
s2.quantification_value(), s2.radio_add_offsets()  # 10000, {B12: -1000}  today
# R_b = (DN_b + offset_b) / Q
R = s2.load()                                     # reflectance           today
# m = 1/cos θ_s + 1/cos θ_v, from the 5 km angle grids
sza, vza = s2.sza, s2.vza["B12"]                  # GeoTensors            today
m = 1 / np.cos(np.radians(sza)) + 1 / np.cos(np.radians(vza))           # user
# R_b = π L_b d² / (E_b cos θ_s)
emit = EMITImage(l1b_path); R_h = emit.load(as_reflectance=True)        # today
srf = read_srf("S2A")              # B11, B12 at 1 nm, offline            today

The running example: one plume, two sensors¶

A 4 km × 4 km Sentinel-2 scene at 20 m, 200 × 200 pixels. The surface has a smooth band 11 reflectance between 0.20 and 0.40 and two cover types, with band 12 at 0.84 and 0.80 of band 11 and a 1 % per-pixel scatter. A Gaussian plume, elongated and rotated by 30°, peaks at $\Delta\Omega = 0.3$ kg/m². The background is 1.9 ppm, a column of 0.0109 kg/m², and the sun and view zenith angles are 35° and 5°, so (6) gives $m = 2.225$. Each band adds noise with a standard deviation of 0.0015 in reflectance. Stage 5 adds an EMIT-like cube.

Stage 1 · The forward model (georeader → radtrax)¶

TL;DR. Build the exact model (1) with Sentinel-2’s own response functions.

1
2
3
4
5
6
7
8
9
10
# Shapes: G = 100,000 wavelengths (0.01 nm); H = W = 200
srf = read_srf("S2A")                     # (877, 13) at 1 nm         today
f11, f12 = np.interp(WL, srf.index, srf["B11"]), ...        # (G,)     user
# σ̄_b = ∫ f_b σ dλ / ∫ f_b dλ
s11, s12 = band_mean(sigma, f11), band_mean(sigma, f12)     # 0.04, 0.40  user
# T_b(m, Ω) = ∫ f_b exp(−m σ Ω) dλ / ∫ f_b dλ
T = radtrax.band_transmittance(sigma, f12)          # (m, Ω) → T   proposed (radtrax)
# R_b = ρ_b T_b(m, Ω_bg + ΔΩ) + ε_b
R12 = rho12 * T(m, OMEGA_BG + d_omega) + eps        # (H, W)       user
# T₁₂(m, Ω_bg) = 0.990: the background takes 1 % of band 12

Sentinel-2A’s band 12 is above half its peak from 2113 to 2286 nm, and band 11 reaches 1668 nm. The background methane already removes 1 % of band 12’s light along this path.

Stage 2 · The effective cross-section against the exact band (radtrax)¶

TL;DR. Compare the linear optical depth (3) with the exact one, the curve of growth.

Problem. The plume’s exact band optical depth, above the background, is

and the linear model claims $\tau_b = m \bar\sigma_b \Delta\Omega$. A line’s core saturates long before its wings do, so the band absorbs less than the linear model says.

The exact optical depth (7) falls short of the linear one (Figure 2) by 3.8 % at 0.05 kg/m², 15 % at the plume’s peak (0.3 kg/m², where band 12 dips by 20 %, not the 23 % the linear model gives), 22 % at 0.5 and 36 % at 1 kg/m². Band 11 absorbs about a tenth as much, which is the ratio this spectrum was built with.

Stage 3 · The multispectral proxy baseline (geotoolz, radtrax)¶

TL;DR. Estimate the plume-free band 12 from band 11 two ways, invert linearly, and integrate the plume’s mass.

Problem. Under the plume, $\rho_{12} \tau_{\text{bg},12}$ is unknown. Band 11 sees the same surface. (a) The source note fits $R_{12} \approx k R_{11} + b$ on pixels outside a suspected-plume mask and uses $k R_{11} + b$ as the baseline. (b) Varon et al.'s multi-band single-pass (MBSP) retrieval Varon et al., 2021, whose estimator radtrax · 2 derives, fits $R_{11} \approx c R_{12}$ through the origin over the whole scene and uses the fractional change

to first order, which keeps band 11’s own absorption. Both are scored by the plume’s mass above the median of a ring of clean pixels around it, as the integrated mass enhancement method does.

1
2
3
4
5
6
7
8
9
10
# Shapes: H = W = 200; plume mask 3,777 pixels; true mass 108,979 kg
# (a) R₁₂ ≈ k R₁₁ + b on clean pixels;  τ = −ln(R₁₂ / (k R₁₁ + b))
k, b = np.polyfit(R11[~suspect], R12[~suspect], 1)  # 0.816, 0.0016  user
omega_a = -np.log(R12 / (k * R11 + b)) / (m * s12)  # B11 "blind"    user
# → 80 % of the mass
# (b) R₁₁ ≈ c R₁₂ through the origin;  ΔΩ = −ln(c R₁₂ / R₁₁) / (m (σ̄₁₂ − σ̄₁₁))
c = (R11 * R12).sum() / (R12 * R12).sum()           # 1.220          user
omega_b = -np.log(c * R12 / R11) / (m * (s12 - s11))                  # user
# → 89 % of the mass; noise off the plume 0.029 (a), 0.032 (b) kg/m²
gz.plume.SBMP()(scene, reference_scene)   # the multi-pass ratio      today

The note’s recipe loses a fifth of the plume. Two causes add: band 11 absorbs under the plume, so the baseline is too low there by about $\bar\sigma_{11} / \bar\sigma_{12}$ of the signal, a tenth here; and saturation (stage 2) takes most of the rest. MBSP removes the first by subtracting $\bar\sigma_{11}$, and recovers 89 %. Neither baseline can model the second surface cover: its patches, where band 12 is 0.80 of band 11 rather than 0.84, come out as false enhancements of 0.04 kg/m² (Figure 3). That surface term, not the sensor noise, sets the noise floor.

Stage 4 · A look-up table of the band observable (radtrax)¶

TL;DR. Tabulate the MBSP observable exactly over air-mass factor and enhancement, and invert each pixel by interpolation.

Problem. Replace the first-order relation in (8) by the exact band observable

tabulated once, and solve $g(m, \Delta\Omega) = -\ln(c R_{12} / R_{11})$ for $\Delta\Omega$ in each pixel.

1
2
3
4
5
6
7
8
# Shapes: table (7, 81) over m ∈ [2, 3.5] and ΔΩ ∈ [0, 0.4] kg/m²
# g(m, ΔΩ) from the exact band transmittances, both bands
lut = radtrax.Table(observable="mbsp", sensor="S2A",
                    m=np.linspace(2, 3.5, 7),
                    omega=np.linspace(0, 0.4, 81))                 proposed (radtrax)
# ΔΩ = g⁻¹(−ln(c R₁₂ / R₁₁); m), odd-extended for negative noise
omega_lut = lut.invert(-np.log(c * R12 / R11), m)  # (H, W)    proposed (radtrax)
# → 96 % of the mass; RMSE over the plume 0.0175 kg/m² (linear: 0.0206)

The table recovers 96 % of the mass. The remaining 4 % comes from the ring’s median, which the surface patches pull below zero. Because $g$ of (9) is increasing in $\Delta\Omega$ at every $m$ in the table, the inversion is a one-dimensional interpolation. The note’s table has water vapour as a second index. That matters where water and methane lines overlap in band 12; this synthetic spectrum has no water, so it does not test it.

Stage 5 · Hyperspectral: a continuum against a matched filter (georeader, geotoolz)¶

TL;DR. On an EMIT-like cube, fit the note’s quadratic continuum and compare it with geotoolz’s matched filter, over a surface with a carbonate absorption.

Problem. A 120 × 120 cube with 68 bands from 2000 to 2500 nm, 7.4 nm apart with an 8.5 nm full width at half maximum. Surfaces follow a quadratic in wavelength, and 15 % of pixels carry calcite’s absorption near 2340 nm, inside the methane band. The plume peaks at 0.15 kg/m². (a) The note’s continuum: a quadratic per pixel through the 2120–2150 and 2430–2470 nm windows, then a least-squares column over the 27 bands from 2200 to 2400 nm. (b) The matched filter Thompson et al., 2015, with the linearised target $t = -\mu \odot m \bar\sigma$ for the background mean $\mu$:

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: cube (C, H, W) = (68, 120, 120); 9 window bands, 27 fit bands
srf_h = reflectance.srf(centres, fwhm, WLH)         # (G, C) Gaussian   today
# (a) R / (a λ² + b λ + c) fitted on the windows; ΔΩ by least squares
omega_cont = continuum_column(cube, windows, fit)   # (H, W)            user
# (b) Δ̂Ω = (x − μ)ᵀ Σ⁻¹ t / tᵀ Σ⁻¹ t,  t = −μ ⊙ m σ̄
mf = gz.matched_filter.MatchedFilter(target=t, fit_on_call=True, axis=0)  # today
omega_mf = mf(cube)                 # background includes the plume     today
masked = cube.copy(); masked[:, suspect] = np.nan   # NaN: left out of the fit
mf_bg = gz.matched_filter.MatchedFilter(target=t, axis=0).fit(masked)    # today
omega_mf_bg = mf_bg.transform(cube)                                     # today
# mass: continuum 95 %, matched filter 72 % (plume in Σ), 91 % (masked)
# over carbonate: continuum +0.0069, matched filter +0.0002 kg/m²

The continuum fit recovers 95 % of the mass but reads carbonate as methane (Figure 4): pixels with calcite come out 0.0069 kg/m² high, against 0.0029 kg/m² low over the rest, and its noise off the plume is 0.0095 kg/m². The matched filter (10) learns the carbonate from the scene covariance and leaves no bias there (0.0002 kg/m²), with a noise of 0.0016 kg/m². Its one trap is the plume itself: with the plume inside the background statistics, the filter recovers only 72 % of the mass; with the suspected plume masked out of the fit, 91 %. The rest is saturation, as in stage 2, because the target is linearised. At the peak, 0.15 kg/m² mixed over a 500 m layer is a 442 ppm enhancement (xrtoolz’s column_mass_to_delta_vmr).

What radtrax would own¶

Where each piece is today is in the package boundaries inventory; Table 2 keeps what radtrax would hold.

Written in JAX, the band transmittance is differentiable, so the table becomes optional: a Gauss–Newton step on (1) inverts each pixel directly, as Report 17’s Laplace step does on the five-band sensor, and a learned retrieval can be trained through it (Report 5).

Open questions¶

1. Are the note’s continuum windows clean? 2120–2150 and 2430–2470 nm need checking against HITRAN’s methane and water lines at EMIT’s resolution; the long-wave window sits at the edge of methane’s 2.3 μm band.

2. Which band 11 to band 12 ratio does Sentinel-2 really have? Stage 3’s bias scales with $\bar\sigma_{11} / \bar\sigma_{12}$; the tenth used here is a choice, and the real ratio follows from HITRAN and the two response functions.

3. A table or a solver? A table is fast and auditable; a differentiable forward model is exact at every pixel and composes with learning. radtrax could ship both, the table built from the forward model.

4. Which atmosphere? The single layer with a geometric air-mass factor ignores scattering and the vertical profile. When does a layered model pay for itself on Sentinel-2’s bands?

Upcoming work¶

References¶

• The source note on the forward model, provenance, proxies and look-up table Johnson, 2026.

• Varon et al. on Sentinel-2 methane retrievals Varon et al., 2021, and on the integrated mass enhancement Varon et al., 2018.

• Thompson et al. on the matched filter for imaging spectroscopy Thompson et al., 2015.

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