The SBMP, MBSP and MBMP family

The radiative transfer equation in the SWIR¶

TL;DR. Two terms of the radiative transfer equation vanish in the shortwave infrared, and what is left integrates to Beer–Lambert.

Problem. For the intensity $I_\lambda$ in a plane-parallel atmosphere, the radiative transfer equation is

where $\mu = \cos\theta$ is the cosine of the zenith angle $\theta$, $\tau_\lambda$ is the optical depth and $J_\lambda$ is the source function (emission plus scattering into the beam), and all of them are scalars at the wavelength $\lambda$. Between 1.6 and 2.3 μm both source terms are small. For thermal emission, a 300 K surface emits $2 \times 10^{-4}$ of the sunlight it reflects at 2.2 μm (with an albedo of 0.2 and $\mathrm{SZA} = 35°$). For Rayleigh scattering, the optical depth falls as $\lambda^{-4}$, so it is 16 times weaker at 2 μm than at 1 μm, and it is $3.7 \times 10^{-4}$ in band 12. So with $J_\lambda \to 0$, (1) reduces along the path $s$ to

where $\sigma_\lambda$ is the absorption cross-section and $n(s)$ is the absorber’s density. We integrate (2) down to a Lambertian surface of albedo $A_\lambda$ and back up, through the column $\Omega = \int n \mathrm{d}z$ along the slant paths, which gives the TOA reflectance

with the geometric air-mass factor

where $L_\lambda$ is the radiance at the sensor, $E_{\lambda,\odot}$ is the solar irradiance, and $\mathrm{SZA}$ and $\mathrm{VZA}$ are the sun and view zenith angles, as in radtrax · 1. As in radtrax · 1, $\Omega$ is the methane column in kg/m² (not a spatial domain), and every symbol here is a scalar per pixel and wavelength. So every pixel has two unknowns per wavelength, $A_\lambda$ and $\Omega$, and the rest of the page is about removing $A_\lambda$.

1
2
3
4
5
6
# Shapes: scalars; λ in μm
# B_λ(T) = 2hc² / λ⁵ / (exp(hc / λkT) − 1)
planck(2.2, 300.0) / L_reflected      # → 2.0 × 10⁻⁴           user
# τ_R(λ) ≈ 0.008569 λ⁻⁴ (1 + 0.0113 λ⁻² + 0.00013 λ⁻⁴)
tau_rayleigh(1.0) / tau_rayleigh(2.0) # → 16.1                 user
tau_rayleigh(2.2)                     # → 3.7 × 10⁻⁴           user

The band model and its linearisation¶

TL;DR. Average over the band, split the column into background and plume, and see how far a straight line carries.

Problem. Band $b$ integrates over its spectral response $f_b(\lambda)$. With the albedo flat across the band, $R_b = A_b T_b(\Omega)$, where the band transmittance $T_b(\Omega) \in [0, 1]$ is the band average of $e^{-\mathrm{AMF} \sigma_\lambda \Omega}$ at this page’s fixed $\mathrm{AMF}$ (the exact band model of radtrax · 1). We split the column as $\Omega = \Omega_\text{bg} + \Delta\Omega$, with $\Omega_\text{bg} = 0.0109$ kg/m² (1.9 ppm). The exact plume optical depth of band $b$ is $\tau_b(\Delta\Omega) = -\ln[T_b(\Omega_\text{bg} + \Delta\Omega) / T_b(\Omega_\text{bg})]$, and its tangent at the background is the sensitivity

so that, to first order in $\Delta\Omega$,

where $\alpha_b \in \mathbb{R}_{>0}$ is in m²/kg, like $\bar{\sigma}_b$, and $R_b^0 \in \mathbb{R}_{\ge 0}$ is the plume-free reflectance.

1
2
3
4
5
6
7
8
9
10
11
# Shapes: N_λ = 100,000 wavelengths (0.01 nm); AMF = 2.225
srf_s2 = read_srf("S2A")         # 1 nm table, one column per band  today
# f_b(λ) on the 0.01 nm grid, zero outside the table        (N_λ,) each
f11 = np.interp(WL, srf_s2.index, srf_s2["B11"], left=0, right=0)  # user
f12 = np.interp(WL, srf_s2.index, srf_s2["B12"], left=0, right=0)  # user
# σ̄_b = ⟨σ⟩_b: 0.040 and 0.400 m²/kg (chosen)
# α_b = −(1/AMF) ∂ ln T_b / ∂Ω at Ω_bg     → 0.0398, 0.3953  user
# AMF α₁₂ ΔΩ / τ₁₂(ΔΩ) at ΔΩ = 0.3 kg/m²   → 1.16           user
# nB_b(ΔΩ) = ⟨exp(−AMF σ ΔΩ)⟩_b, no background       (181,)
nb = build_nb_lut(ds, srf, T_K=280, p_atm=0.8, amf=m)       # today
# −ln nB₁₂ / τ₁₂ at 0.3 kg/m²                → 1.011          user

$\alpha_b$ sits just below $\bar{\sigma}_b$, because the background has already darkened the strongest line cores. At 0.3 kg/m² the straight line is 16 % above the exact optical depth (the curve of growth). Band 11, at the edge of the 1.66 μm band, absorbs a tenth as much as band 12 on this spectrum, which is a modelling choice (though it is not zero on real spectra either). plumax’s band table, nB, integrates the plume without the background column, which overstates band 12’s plume optical depth by 1.1 %.

Three ways to remove the surface¶

TL;DR. Ratio across passes (SBMP), across bands (MBSP), or both (MBMP); each ratio cancels what it assumes is constant.

Problem. We write (6) for bands 11 and 12, on a plume pass $(p)$ and a plume-free reference pass $(r)$ over the same pixels. Each of the three estimators turns the reflectances into a log signal $\Delta\tau \in \mathbb{R}$ per pixel. SBMP compares band 12 across the two passes,

MBSP compares the two bands within the plume pass,

with the fractional band difference

and MBMP compares the band difference across the two passes,

where $r = A_{12} / A_{11} \in \mathbb{R}_{>0}$ is the band ratio. [CHECK: taking $-\ln$ of the ratios gives the leftover log terms with a minus sign, e.g. $\tau_{12} - \ln(c_s A_{12}^{(p)} / A_{12}^{(r)})$; the sign does not change what each method assumes away.] Each scale factor is a scalar least-squares fit through the origin over background pixels, $c_s$ from $R_{12}^{(r)} \approx c_s R_{12}^{(p)}$ and $c$ per pass from $R_{11} \approx c R_{12}$. The leftover log is what each method assumes away. SBMP needs the band-12 albedo unchanged between passes, up to one scene-wide factor. MBSP needs the band ratio $r$ to be the same everywhere in the scene, up to $c$. MBMP only needs $r$ unchanged at each pixel between passes, up to one scene-wide factor per pass, so it is a difference in differences (across bands and then across dates). The column comes from inverting the exact table, $\Delta\Omega = \tau_{12}^{-1}(\Delta\tau_\text{SBMP})$ for SBMP and $\Delta\Omega = g^{-1}(\Delta\tau)$ with $g = \tau_{12} - \tau_{11}$ for MBSP and MBMP, or, to first order, from $\Delta\tau_\text{SBMP} / (\mathrm{AMF} \alpha_{12})$ and $\Delta\tau / (\mathrm{AMF} \Delta\alpha)$ with $\Delta\alpha = \alpha_{12} - \alpha_{11}$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Shapes: obs, ref (2, 200, 200) = (B11, B12) × (H, W)
# c = Σ x y / Σ x² over background pixels                 () per fit
c_s = fit_c(R12r, R12p, bg)                                # user
c_p, c_r = fit_c(R11p, R12p, bg), fit_c(R11r, R12r, bg)    # user
# δ = c R₁₂ / R₁₁ − 1, per pass               (200, 200) → (200, 200)
d_p, d_r = c_p * R12p / R11p - 1, c_r * R12r / R11r - 1    # user
# ΔΩ = τ₁₂⁻¹(−ln(c_s R₁₂ / R₁₂′))              (200, 200) → (200, 200)
sbmp = invert(-np.log(c_s * R12p / R12r), tau12_tab)       # user
# ΔΩ = g⁻¹(−ln(1 + δ))                          (200, 200) → (200, 200)
mbsp = invert(-np.log1p(d_p), g_mb)                        # user
# ΔΩ = g⁻¹(−ln((1 + δ) / (1 + δ′)))             (200, 200) → (200, 200)
mbmp = invert(-(np.log1p(d_p) - np.log1p(d_r)), g_mb)      # user
# ln(R₁₁/R₁₂) − ln(R₁₁′/R₁₂′), no scale factors   (2, 200, 200) → ()
gz = gz.plume.SBMP(swir1=0, swir2=1, axis=0)(obs, ref)     # today

The running example: one scene, every nuisance¶

TL;DR. One plume, two passes, and every surface and atmospheric change the derivations warn about, each in its own patch.

Problem. A 200 × 200 scene with 20 m pixels holds a plume with a peak of 0.3 kg/m² (148,700 kg) on the plume pass only. Both passes share a textured surface: the albedo varies by 15 %, the band ratio $r$ by 4 %, and a mineral patch has $r$ 15 % low. Between the passes, the whole scene brightens by 5 %, a field is wetted (both bands × 0.8), one patch’s ratio drops by 3 %, and band 12’s water-vapour transmission drops by 1 % everywhere. The pixels enter through plumax’s inject_plume, with the band table built against the background column.

We find the hierarchy the derivations describe (Figure 1, Table 1). SBMP, with or without a scale factor, reads the wetted field as 0.28 to 0.29 kg/m² of methane, which is nearly the plume’s peak. MBSP reads the stable mineral patch as 0.21 kg/m² and the 4 % ratio texture as 0.054 kg/m² of clutter. MBMP removes both, and the scene-wide brightening and the water-vapour change with them. What it keeps is the patch whose band ratio changed, where 3 % reads as 0.039 kg/m². The plume’s mass then depends on two choices which the derivations make differently. A scale factor fitted over every pixel takes in the plume, moves by 0.31 % and leaves a −0.009 kg/m² offset that costs 11 % of the mass. In addition, the slope that converts the signal to a column matters: the exact table recovers 99 %, a linear inversion 91 %, the $\delta$ difference 88 %, and $\bar{\sigma}_{12}$ with $\sigma_{11} = 0$ 78 %. geotoolz’s SBMP keeps the 1 % water-vapour change as a +0.013 kg/m² offset over the whole scene, which is 16 % extra mass on this plume (why).

Noise¶

TL;DR. Each ratio cancels a surface term and adds the noise of the bands it uses.

Problem. With independent noise of relative size $\varepsilon \in \mathbb{R}_{>0}$ in each band and pass, the log ratios add variances. SBMP uses two band-12 images, MBSP one image of each band and MBMP four images, and the two-band methods divide by the smaller $\Delta\alpha$. So, to first order, the standard deviations of the retrieved columns are

1
2
3
# Shapes: obs, ref (2, 200, 200); ε = 0.01 per band and pass
# sd(ΔΩ̂_noisy − ΔΩ̂_clean) off the plume       (200, 200) → ()
# → SBMP 0.016, MBSP 0.018, MBMP 0.025 kg/m²; theory within 3 %

As (11) and (13) predict, MBMP pays for its robustness with 1.6 times SBMP’s noise. On a surface that is stable in time, SBMP with a scale factor is the more precise retrieval, and on a surface that is uniform in its band ratio, MBSP needs no reference pass. MBMP is the default because real surfaces are neither, and its noise falls with averaging whereas a surface artefact does not.

When MBMP is not the right choice¶

Derivation (1) lists three cases, and the run bears on each:

1. No good reference pass. MBMP and SBMP both need a plume-free date with the same surface. If the only reference is from another season, the local ratio changes grow, and MBSP (which needs no reference) may do better.

2. A strongly wavelength-dependent surface. Where $r$ varies widely across the scene, one $c$ fits it poorly. MBMP still cancels the stable part of $r$ (the mineral patch is 0.000 kg/m² above the offset), but the fitted $c$ is less well defined, and SBMP avoids the band ratio.

3. Water vapour that changes between passes. Band 12 has water-vapour lines. A scene-wide change is a scene-wide factor, which the per-pass $c$ absorbs. However, a change that varies across the scene (a front, terrain) is not absorbed, and it reads like a changed band ratio, where 1 % is 0.013 kg/m².

What each method cancels¶

Table 2 collects the run: what each method removes, what it leaves, and its noise.

Where it lives in the stack¶

The package boundaries inventory lists the pieces: plumax’s radtran.nb_lut (build_nb_lut, lookup_nb, inject_plume) for the band tables, georeader’s read_srf for the real response functions (B11 from 1569 to 1658 nm and B12 from 2113 to 2286 nm at half maximum), and geotoolz’s plume.SBMP for the two-pass, two-band log signal. The scale factors, the background masks and the exact inversion exist nowhere in the stack, so we write them in the notebook.

Open questions¶

1. What should the operators be called? An MBMP operator in radtrax, with SBMP and MBSP beside it, each with the scale factors and a background mask, would match the literature, and the current SBMP could become MBMP’s scale_factors=False mode.

2. Where does the background mask come from? In the run, we use a first-pass detection grown by a Gaussian of 2 pixels. The alternatives are an iterated fit or the plume mask from geotoolz’s PlumeMask, and the choice moves the offset.

3. Should the band table include the background? plumax’s nB omits it. Including it is one extra integral and changes band 12 by 1 % here, and more over a high background or long air-mass paths.

Upcoming work¶

Table 3 lists the follow-ups.

References¶

• The two derivations: “Deriving the MBMP method family from first principles” Johnson, 2026 and the step-by-step optical derivation Johnson, 2026.

• The SBMP, MBSP and MBMP retrievals Varon et al., 2021; the IME and the fine-scale plume retrievals they feed Varon et al., 2018; a ground-truth calibration of Sentinel-2 retrievals Chen & others, 2022; transient emissions from geostationary observations Watine-Guiu et al., 2023.

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