Quantifying methane point sources from fine-scale satellite observations of atmospheric methane plumes

Anthropogenic methane emissions originate from a large number of relatively small point sources. The planned GHGSat satellite fleet aims to quantify emissions from individual point sources by measuring methane column plumes over selected $\sim \mathrm{10}\times \mathrm{10}$ km² domains with $\le \mathrm{50}\times \mathrm{50}$ m² pixel resolution and 1 %–5 % measurement precision. Here we develop algorithms for retrieving point source rates from such measurements. We simulate a large ensemble of instantaneous methane column plumes at 50×50 m² pixel resolution for a range of atmospheric conditions using the Weather Research and Forecasting model (WRF) in large eddy simulation (LES) mode and adding instrument noise. We show that standard methods to infer source rates by Gaussian plume inversion or source pixel mass balance are prone to large errors because the turbulence cannot be properly parameterized on the small scale of instantaneous methane plumes. The integrated mass enhancement (IME) method, which relates total plume mass to source rate, and the cross-sectional flux method, which infers source rate from fluxes across plume transects, are better adapted to the problem. We show that the IME method with local measurements of the 10 m wind speed can infer source rates with an error of 0.07–0.17 t h${}^{-\mathrm{1}}+\mathrm{5}$ %–12 % depending on instrument precision (1 %–5 %). The cross-sectional flux method has slightly larger errors (0.07–0.26 t h${}^{-\mathrm{1}}+\mathrm{8}$ %–12 %) but a simpler physical basis. For comparison, point sources larger than 0.3 t h⁻¹ contribute more than 75 % of methane emissions reported to the US Greenhouse Gas Reporting Program. Additional error applies if local wind speed measurements are not available and may dominate the overall error at low wind speeds. Low winds are beneficial for source detection but detrimental for source quantification.

1 Introduction

Satellite instruments can measure atmospheric methane columns from solar backscatter in the shortwave infrared (SWIR) with near-uniform sensitivity down to the surface (Frankenberg et al., 2005). There is considerable interest in using these measurements to quantify methane emissions (Jacob et al., 2016). Most current and planned instruments have pixel resolutions of 1–10 km and column precisions of 0.1 %–1 % (Bovensmann et al., 1999; Butz et al., 2011; Veefkind et al., 2012; Polonsky et al., 2014; Kuze et al., 2016). Jacob et al. (2016) show that these measurements can successfully map regional methane emissions but have limited ability to resolve individual methane point sources, even with imaging capabilities, because the sources tend to be relatively small and spatially clustered (e.g., oil/gas fields, livestock operations, landfills, coal mine vents). The GHGSat microsatellite fleet (Germain et al., 2017; McKeever et al., 2017) aims to address this gap by observing methane columns over selected scenes of order 10×10 km² with $\le \mathrm{50}\times \mathrm{50}$ m² effective pixel resolution and moderate precision (1 %–5 %). Here we present algorithms for interpreting the instantaneous plumes observed by this instrument in terms of the implied point source (facility-level) emissions and estimate the associated errors and detection limits as a function of instrument precision.

Aircraft remote sensing of methane columns over oil/gas and coal mining facilities shows that the instantaneous plumes have irregular shapes and detectable sizes of order 0.1–1 km (Thorpe et al., 2016; Thompson et al., 2015, 2016; Frankenberg et al., 2016). A standard method used to retrieve source rates from plume observations is to assume Gaussian plume behaviour, as expected from statistically averaged turbulence (Bovensmann et al., 2010; Krings et al., 2011, 2013; Rayner et al., 2014; Fioletov et al., 2015; Nassar et al., 2017; Schwandner et al., 2017). This method may induce large errors for small instantaneous plumes, which generally do not follow steady-state Gaussian behaviour. Several authors have addressed this difficulty. Krings et al. (2011, 2013) proposed a cross-sectional flux method to derive the source rate as the product of the local wind and the concentration integrated over a plume cross section, expanding on a similar method used for in situ plume measurements (White et al., 1976; Cambaliza et al., 2014; Conley et al., 2016). Jacob et al. (2016) described a mass balance method for inferring the source rate solely based on the enhancement in the source pixel. Frankenberg et al. (2016) inferred the source rate empirically from the total detectable mass of methane in the plume (integrated mass enhancement or IME).

A common feature of all these methods for retrieving the point source rate Q from plume observations is their need for independent knowledge of the wind speed U driving transport of the plume. In the cross-sectional flux method applied to in situ aircraft observations, methane and local wind speed are measured concurrently (Conley et al., 2016). In remote sensing, however, the wind speed for the instantaneous column plume is not directly measured and may be variable both vertically and horizontally across the plume.

Here we use observing system simulation experiments to develop algorithms for retrieving individual point source rates from fine-scale satellite observations of instantaneous methane plumes. We review previously used plume inversion methods and show with large eddy simulations (LES) that the IME and cross-sectional flux methods are best suited to the problem. We further develop the IME method to provide a physical basis for its general application. We consider different combinations of instrument precision, meteorological environment, and wind information to test the methods and quantify errors. Our work is motivated by GHGSat but is more generally applicable to any fine-scale plume observations from space.

2 Review of methods for retrieving point sources from observations of column plumes

A methane point source produces a turbulent plume of atmospheric methane with characteristics determined by the strength of the source, the wind field, and turbulence that depends on atmospheric stability and surface roughness. Four different methods have been proposed to quantify point source rates from plume observations: (1) the Gaussian plume inversion method (Bovensmann et al., 2010; Krings et al., 2011, 2013; Rayner et al., 2014; Fioletov et al., 2015; Nassar et al., 2017; Schwandner et al., 2017), (2) the source pixel method (Jacob et al., 2016; Buchwitz et al., 2017), (3) the cross-sectional flux method (White et al., 1976; Conley et al., 2016; Krings et al., 2011, 2013; Tratt et al., 2011, 2014; Frankenberg et al., 2016), and (4) the IME method (Thompson et al., 2016; Frankenberg et al., 2016). Here we discuss these methods for remote sensing observations of column plumes. This is a somewhat different problem than for in situ observations of plumes. In situ observations benefit from a stronger signal but require characterization of the plume in the vertical dimension, which is integrated in a column measurement.

Satellite remote sensing of methane plumes retrieves column concentrations with vertical sensitivity that depends on atmospheric scattering and absorption. Clear-sky observations in the SWIR have near-unit sensitivity throughout the tropospheric column, while observations in the thermal infrared (TIR) have strong vertical dependence determined by temperature contrast with the surface (Worden et al., 2013). Here we focus on SWIR observations, where we can ignore vertical dependence in sensitivity. TIR remote sensing has been used effectively to detect methane plumes from low-flying aircraft (Tratt et al., 2014; Frankenberg et al., 2016) but is not practical from space because of interference from the background methane column above the plume (Jacob et al., 2016).

Methane column concentrations retrieved from remote sensing are commonly expressed in the literature as column-average dry molar mixing ratio X [ppb]. The plume is then characterized by an enhancement $\mathrm{\Delta }X=X-X_{\mathrm{b}}$ relative to the local background X_b. For our purposes of relating plume observations to the source rate Q [kg s⁻¹], a more useful measure of plume concentration is the column mass enhancement ΔΩ with units [kg m⁻²]. ΔΩ is related to ΔX by

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }\mathrm{\Omega }={\displaystyle \frac{M_{{\mathrm{CH}}_{\mathrm{4}}}}{M_{\mathrm{a}}}}{\mathrm{\Omega }}_{\mathrm{a}}\mathrm{\Delta }X,\end{array}$$

where $M_{{\mathrm{CH}}_{\mathrm{4}}}$ and M_a are the molar masses of methane and dry air [kg mol⁻¹] and Ω_a is the column of dry air [kg m⁻²].

2.1 Gaussian plume inversion method

The Gaussian plume inversion method fits a Gaussian plume model to the measured columns. Assuming a steady wind U oriented along the x axis and integrating the three-dimensional Gaussian plume equation vertically, one obtains an expression for Q in terms of the vertical column enhancement ΔΩ(x,y) downwind of a point source located at the origin (Bovensmann et al., 2010):

$$\begin{array}{}\text{(2)}& {\displaystyle }Q=U\mathrm{\Delta }\mathrm{\Omega }\left(x,y\right)\left(\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_y\left(x\right)e^{\frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y(x{)}^{\mathrm{2}}}}\right).\end{array}$$

The empirical dispersion parameter σ_y(x) [m] describes the horizontal spread of the plume along the y axis orthogonal to the wind direction. It is commonly parameterized as (Martin, 1976)

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\sigma }}_y\left(x\right)=a{\left({\displaystyle \frac{x}{x_{\mathrm{0}}}}\right)}^{\mathrm{0.894}},\end{array}$$

where x₀=1000 m, and the dispersion coefficient a [m] depends on atmospheric stability as defined by the Pasquill–Gifford atmospheric stability categories (Pasquill, 1961). The solution to Eq. (2) may involve non-linear optimal estimation fitting of a to the observed plume (Krings et al., 2011). The fit may not be successful if the instantaneous plume shows large departure from steady-state Gaussian behaviour.

2.2 Source pixel method

In the source pixel method used by Jacob et al. (2016) to compare different satellite-observing configurations, emissions are inferred solely from methane enhancements in the source pixel relative to the local background. For an observation pixel of dimension W [m] containing a methane point source ventilated by a uniform wind speed U [m s⁻¹], the source rate Q [kg s⁻¹] is calculated from the mean source pixel enhancement ΔΩ [kg m⁻²]:

$$\begin{array}{}\text{(4)}& Q={\displaystyle \frac{UWp}{g{\mathrm{\Omega }}_{\mathrm{a}}}}\mathrm{\Delta }\mathrm{\Omega },\end{array}$$

where p is the surface pressure and g is the acceleration of gravity. The source pixel method ignores additional information from the plume downwind and is therefore not optimal. In addition, the instantaneous wind U may have large uncertainty for small pixels because of turbulence. The method is also vulnerable to systematic errors in the column enhancement retrieved over the source pixel (e.g., due to different reflectance properties of the emitter compared to the surrounding area) and errors in the local background estimate.

2.3 Cross-sectional flux method

In the cross-sectional flux method, the source rate is estimated by computing the flux through one or more plume cross sections orthogonal to the plume axis. This approach is commonly used for aircraft in situ observations (White et al., 1976; Mays et al., 2009; Cambaliza et al., 2014, 2015; Conley et al., 2016). Krings et al. (2011, 2013) and Tratt et al. (2011, 2014) extended it to methane columns observed by aircraft remote sensing. By mass balance, the source rate Q must be equal to the product of the wind speed and the column plume transect along the y axis perpendicular to the wind:

$$\begin{array}{}\text{(5)}& Q={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}U(x,y)\mathrm{\Delta }\mathrm{\Omega }(x,y)\mathrm{d}y,\end{array}$$

where the integral is approximated in the observations as a discrete summation of the product $U(x,y)\mathrm{\Delta }\mathrm{\Omega }(x,y)$ over the detectable width of the plume.

Compared to in situ aircraft measurements, an advantage of remote sensing is that the full vertical extent of the plume is covered by the measurement. A disadvantage is that the wind U(x,y) is not as well characterized: it must describe some vertical average over the plume extent and there is generally no information on its horizontal variability over the scale of the plume. This may require estimation of an effective wind speed U_eff applied to the cross-plume integral C [kg m⁻¹] of the column along the y axis:

$$\begin{array}{}\text{(6)}& {\displaystyle }Q=C{U}_{\mathrm{eff}},\mathrm{with}C={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{\Delta }\mathrm{\Omega }(x,y)\mathrm{d}y.\end{array}$$

If U_eff is assumed uniform with distance x downwind of the source, then the integral C is independent of x and can be computed for different values of x to improve accuracy through averaging. We show in Sect. 6 how to estimate U_eff for use in Eq. (6).

2.4 Integrated mass enhancement (IME) method

The IME method relates the source rate to the total plume mass detected downwind of the source. The IME of an observed column plume comprising N pixels of area A_j (j=1…N) is

$$\begin{array}{}\text{(7)}& \mathrm{IME}={\sum }_{j=\mathrm{1}}^N\mathrm{\Delta }{\mathrm{\Omega }}_j{A}_j.\end{array}$$

Frankenberg et al. (2016) defined an empirical linear relationship between IME and Q for their methane plumes detected from aircraft by using independent estimates of a few sources from the cross-sectional flux method. They then applied this linear relationship to all their observed plumes.

More fundamentally, the relationship between IME and Q is defined by the residence time τ of methane in the detectable plume: $Q=\mathrm{IME}/\mathit{\tau }$. One can express τ dimensionally in terms of an effective wind speed U_eff [m s⁻¹] and a plume size L [m]:

$$\begin{array}{}\text{(8)}& Q={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\mathrm{IME}={\displaystyle \frac{U_{\mathrm{eff}}}{L}}\mathrm{IME}={\displaystyle \frac{U_{\mathrm{eff}}}{L}}{\sum }_{j=\mathrm{1}}^N\mathrm{\Delta }{\mathrm{\Omega }}_j{A}_j.\end{array}$$

U_eff and L would have simple physical meanings of wind speed and plume length if dissipation of the plume occurred by uniform transport to a terminal distance downwind of the source. But the actual mechanism for plume dissipation is turbulent diffusion, which takes place in all directions. U_eff and L must therefore be viewed as operational parameters to be related to observations of wind speed U and plume extent. In Sect. 5 we derive these relationships from synthetic plumes generated by LES. The detectable plume size L depends on Q and U, introducing non-linearity in Eq. (8).

3 Synthetic GHGSat observations of methane plumes

We generate synthetic GHGSat plumes with the Weather Research and Forecasting model in LES mode (WRF-LES; Moeng et al., 2007; Skamarock et al., 2008) to evaluate the ability of the methods described in Sect. 2 to retrieve methane point source rates from satellite observations. The WRF-LES simulations are conducted at 50×50 m² resolution and are sampled virtually with the GHGSat instrument by column integration and with consideration of instrument precision. In this section, we briefly describe the GHGSat instrument and the application of LES to produce synthetic plumes.

3.1 The GHGSat instrument

GHGSat is a lightweight satellite instrument (∼15 kg) designed by GHGSat, Inc. to detect atmospheric methane plumes from individual point sources. A demonstration instrument (GHGSat-D) was launched in June 2016 into sun-synchronous orbit (local solar viewing time of 09:30 on the descending node) to test the instrument performance and column retrieval algorithm. The launch of the first operational satellite is scheduled for early 2019. The long-term plan is for a constellation of sun-synchronous GHGSat microsatellites in low Earth orbit, providing frequent observations of different sources of interest.

GHGSat measures backscattered solar SWIR radiation at 1635–1670 nm (0.1 nm spectral resolution) over $\sim \mathrm{10}\times \mathrm{10}$ km² targeted domains (12×12 km² for GHGSat-D) with $\le \mathrm{50}\times \mathrm{50}$ m² pixels. The design precision for the methane column retrieval is 1 %–5 %. This is coarser precision than that of other satellite instruments (e.g., 0.6 % for TROPOMI, 0.7 % for GOSAT), but for observations of point sources it is more than offset by the finer pixel resolution (Jacob et al., 2016).

3.2 Large eddy simulations (LES)

We apply WRF-LES to simulate turbulent plume transport at 50×50 m² horizontal resolution and 30 m vertical resolution over a 6×6 km² domain. We use a modified version of the WRF v3.8 default LES case with cloud-free conditions and no topography (WRF User's Guide, 2018; Moeng et al., 2007). Simulations are performed with one-way nesting from an external simulation over a 7.5×7.5 km² domain with 150×150 m² resolution and periodic boundary conditions. A uniform sensible heat flux H=100 W m⁻² is applied at the surface to drive buoyant turbulence. Mechanical turbulence is generated by surface drag with an aerodynamic roughness height of 0.1 m. Forcing from a large-scale pressure gradient maintains momentum across the domain.

Each LES simulates 5 h of atmospheric transport. The first 3 h spin up realistic turbulence and the final 2 h are used for analysis. We use a range of initial mixing depths and wind speed soundings to produce different simulations. The potential temperature soundings are uniform at 290 K from the surface to a mixing depth set at either 500, 800, or 1100 m altitude, with an inversion above that altitude and the model top set 700 m above the inversion. For each of these three mixing depths, we conduct simulations using five initially uniform southerly wind profiles with speeds of 2–8 m s⁻¹. The resulting LES ensemble of 15 simulations is broadly representative of the range of meteorological conditions that could be sampled with a SWIR daytime instrument.

We use the WRF-LES passive tracer transport capability (Nottrott et al., 2014; Nunalee et al., 2014) to generate a plume from a single constant point source in the WRF-LES meteorological environment. From there we integrate the plume over vertical columns and add random noise to produce GHGSat pseudo-observations. We archive the tracer column field every 30 s as an independent realization of the instantaneous plume. From the 15 WRF-LES simulations we thus archive a collection of 3600 scenes, representing our statistical ensemble for different possible realizations of turbulence.

The WRF-LES point source in the archived ensemble has a normalized source rate, which we subsequently scale from the output to simulate source rates Q in the range 0.05–2.25 t CH₄ h⁻¹ (0.5–20 kt a⁻¹). This range covers the top 25 % of sources reporting to the US Greenhouse Gas Reporting Program (GHGRP) and contributing 75 % of total GHGRP methane emissions (Jacob et al., 2016). A uniform background methane column of 0.01 kg m⁻² (roughly 1850 ppb) is added to the tracer column. Uncorrelated measurement noise is then added as a random increment of the background sampled from a normal distribution with zero mean bias and standard deviation σ=1 %–5 %, corresponding to the range of expected instrument precision. The column enhancement ΔΩ is then inferred by subtracting the 0.01 kg m⁻² background, which is therefore assumed to be known.

Figure 1 Examples of column plume pseudo-observations generated by an LES on a 50×50 m² grid. The circle identifies the location of the source. Each panel shows the same synthetic plume observation for a source Q=1 t h⁻¹, with instrument precision σ varying from 1 % to 5 %.

Figure 1 shows examples of synthetic plume observations produced in this manner for a source Q=1 t h⁻¹, assuming different levels of instrument precision. As the precision decreases, the plume is increasingly difficult to detect.

4 Inadequacy of the Gaussian plume and source pixel methods

Previous studies of CO₂ column observations from the OCO-2 satellite instrument with $\sim \mathrm{1.3}\times \mathrm{2.25}$ km² nadir pixel resolution (Crisp et al., 2017) have shown that the Gaussian plume inversion method can be effective for quantifying very large CO₂ emissions from power plants (Nassar et al., 2017) and volcanoes (Schwandner et al., 2017). We find here that the approach fails when applied to fine-scale methane plumes, because the plumes depart too much from Gaussian behaviour. CO₂ point sources can be considerably larger relative to background concentrations and instrument precision levels, and the resulting plumes can then be observed over distances of tens of kilometres. Such a large size allows for statistical averaging of eddies and hence better approximation of Gaussian behaviour, even for an instantaneous plume. To demonstrate this, Fig. 2 shows an LES snapshot of a large power plant emitting 3.75 kt CO₂ h⁻¹ in a 72 km domain with 300×300 m² pixel resolution. Fitting a Gaussian plume to the 300 m pixel enhancements yields a coefficient of determination R² of only 0.24, but R² increases to 0.86 when the LES image is averaged over 3×3 km² pixels. Spatial averaging of turbulence over kilometre-scale pixels thus greatly improves the accuracy of Gaussian plume models, but this requires sufficiently large plumes. Methane plumes are generally too small to allow such averaging (Frankenberg et al., 2016).

Figure 2 CO₂ column enhancements relative to background for a 3.75 kt CO₂ h⁻¹ (33 Mt CO₂ a⁻¹) power plant plume simulated by LES at 300 m resolution. (a) shows the plume with 300 m pixel resolution and (b) shows the same plume but with pixel resolution degraded to 3000 m. The coefficient of determination (R²) inset measures the ability to fit each LES plume to a Gaussian form (Eqs. 2–3). (c) shows how the coefficient of determination and the root-mean-square error (expressed as a percentage of the median pixel enhancement) vary with pixel resolution.

The source pixel retrieval method only considers the column enhancement over the point source pixel, thus inferring the source rate from ventilation of that pixel by the local wind. It assumes in effect that the near-field plume is diluted over the source pixel and neglects information from the plume downwind. This can be an effective method when pixel resolution is coarse, so that most of the information is in the source pixel and the mean wind across the pixel can be well defined (Buchwitz et al., 2017). However, it has three major shortcomings when applied to GHGSat 50×50 m² pixels: (1) it does not exploit the information from downwind pixels, where most of the plume mass typically resides; (2) small-scale turbulence generates strong variability in the wind; and (3) source pixel ventilation may take place by turbulent horizontal diffusion rather than advection by the mean wind, leading to negative bias. With regard to (2), the residence time in a GHGSat pixel is only ∼30 s, and there is large variability in the wind on such a short timescale that cannot be described deterministically. For example, in a typical LES under moderately unstable conditions we find a 10 m wind speed of 2.45±0.8 m s⁻¹, where the standard deviation is for the 30 s data. This 30 s variability in wind speed alone thus introduces a factor of 30 % uncertainty in the source estimate. With regard to (3), the relative importance of turbulent diffusion and advection is diagnosed by the Péclet number $\mathit{Pe}=UL/K_{\mathrm{H}}$, where K_H is the turbulent horizontal diffusion coefficient (Brasseur and Jacob, 2017). For a typical K_H=50 m² s⁻¹ (d'Isidoro et al., 2010) with U=2 m s⁻¹ and L=50 m we find Pe∼1, so that turbulent diffusion and advection are of comparable importance.

5 Computing the source rate by the IME method

We showed in Sect. 2.4 how the IME method for retrieving the point source rate Q from the measured IME hinges on knowledge of the residence time of methane in the detectable plume. We refer to this residence time as the plume lifetime $\mathit{\tau }=\mathrm{IME}/Q$, which in turn is related to two parameters: an effective wind speed U_eff and a characteristic plume size L. IME and L can be inferred from the plume observations, while U_eff can be inferred from the observable 10 m wind speed U₁₀ at the point of emission.

5.1 Inferring the plume mass (IME) and size (L)

Inferring IME and L from the plume observations requires that we define the horizontal extent of the plume through a pixel selection procedure that separates signal from noise. Careful selection is important. Consider an array of N pixels of equal area and with retrieved column enhancements ΔΩ_j (j=1…N). If each pixel enhancement includes a contribution s_j from signal (actual plume enhancement) and ε_j from random noise, then as per Eq. (7),

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{IME}}{A}}={\sum }_{j=\mathrm{1}}^N\mathrm{\Delta }{\mathrm{\Omega }}_j={\sum }_{j=\mathrm{1}}^N(s_j+{\mathit{\epsilon }}_j)={\mathit{\epsilon }}_{\mathrm{a}}+{\sum }_{j=\mathrm{1}}^N{s}_j,\end{array}$$

where ε_a is the total (absolute) measurement error. The relative error ε_r is then ${\mathit{\epsilon }}_{\mathrm{r}}={\mathit{\epsilon }}_{\mathrm{a}}/{\sum }_j^N{s}_j$. If the noise is normally distributed and uncorrelated, then the error standard deviation is proportional to $\sqrt{N}$, so that the standard deviation σ_r of the relative error scales as

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{r}}\propto {\displaystyle \frac{\sqrt{N}}{{\sum }_j^N{s}_j}}.\end{array}$$

Now consider two extreme cases: (1) all pixels contain the same signal s₀, and (2) only one pixel contains signal s₀ and the other pixels contain only noise. In case (1), the total signal ${\sum }_j^N{s}_j$ is proportional to N, meaning ${\mathit{\sigma }}_{\mathrm{r}}\propto \mathrm{1}/\sqrt{N}$. By contrast, in case (2), the total signal is equal to s₀, so ${\mathit{\sigma }}_{\mathrm{r}}\propto \sqrt{N}$. Thus, we see that aggregating plume pixels can either decrease or increase the error on the IME depending on whether these pixels have significant signal or not.

Figure 3 Illustration of the procedure for constructing plume masks in the IME method. (a) Satellite pseudo-observation generated by LES for a point source Q=1 t h⁻¹, with instrument precision σ=1 % (same as in Fig. 1). (b) The output of the t test procedure for significant signal. (c) The final plume mask after application of median and Gaussian filters.

Figure 3 illustrates how we construct a plume mask to select plume pixels with significant signal-to-noise ratios. The background distribution (mean ± standard deviation) is first characterized by an upwind sample of the measured columns, mimicking what one would do with actual observations. Next, we sample the 5×5 pixels neighbourhood centred on each pixel in the viewing domain and compare the sample distributions to the background distribution by means of a Student's t test. Pixels with 5×5 neighbourhoods that follow a distribution significantly different than the background at a confidence level of 95 % or higher are assigned to the plume and others to the background. The resulting Boolean plume mask contains some random classification errors, so we smooth it with a median filter followed by a Gaussian filter and thresholding. The median filter replaces each pixel in the mask by the median value of its 3×3 neighbourhood. The Gaussian filter convolves the mask with a two-dimensional Gaussian of standard deviation 2–5, with larger values for higher noise levels.

We compute the IME by summing pixel enhancements within the plume mask following Eq. (7). A simple measure of the plume size L can be taken as

$$\begin{array}{}\text{(11)}& L=\sqrt{A_{\mathrm{M}}},\end{array}$$

where A_M [m²] is the area of the plume mask. Another possible estimate of L would be the mask's perimeter, which can be obtained by contour tracing. The definition of L is not critical as long as it has some physical basis relating it to the observed plume geometry. A different definition would imply a different calculation of U_eff.

5.2 Inferring U_eff from the 10 m wind speed U₁₀

The effective wind speed U_eff is a parameter of the IME method that should be related to the measurable 10 m wind speed at the location of the point source, and here we use the LES to derive the U_eff=f(U₁₀) relationship. If U₁₀ is not actually measured at the site, it can be estimated from an operational meteorological database at the cost of some representation error. We discuss that error in Sect. 7.

We derive the U_eff=f(U₁₀) relationship from a training set of column plumes comprising two-thirds of the LES ensemble selected at random (i.e., 2400 plume instances). The remaining plumes serve as a test set for evaluating the retrieval algorithm. For each plume in the training set, U_eff is computed from Eq. (8) as $U_{\mathrm{eff}}=QL/\mathrm{IME}$, based on the known source rate Q and with L and IME determined from the plume masks. The corresponding U₁₀ time series at the location of the source is obtained from the LES, averaged over the plume lifetime $\mathit{\tau }=\mathrm{IME}/Q$. Values of τ in our ensemble range from 1 to 60 min depending on instrument precision, source rate, and wind speed. In practice, τ is unknown a priori and must be inferred from the plume observations and local wind speed information. We discuss this in Sect. 5.3.

Figure 4 shows the relationship between U_eff and U₁₀ inferred from the LES ensemble. We find that we can fit the data to a logarithmic dependence:

$$\begin{array}{}\text{(12)}& U_{\mathrm{eff}}={\mathit{\alpha }}_{\mathrm{1}}\log U_{\mathrm{10}}+{\mathit{\alpha }}_{\mathrm{2}},\end{array}$$

where α₁=0.9–1.1, α₂=0.6 m s⁻¹, and the range on the first coefficient is for the 1 %–5 % range of instrument precision. For 1 % instrument precision, the logarithmic function plotted in Fig. 4 captures 75 % of the variance (R²=0.75). This decreases to 35 % of the variance for 5 % instrument precision. The convexity of the relationship is an important result, as it implies that error in U_eff is smaller than error in U₁₀. One might expect U_eff from the IME method to be proportional to U₁₀, such that IME∕L would be inversely proportional to U₁₀ as per Eq. (8). However, even though that inverse relationship holds for plume concentrations (see Eq. 6), it is much weaker for the IME because the concentrated plume in the near-field of the source remains in the signal even at high wind speeds. Thus, the plume observations themselves interpreted with the IME method contain some information on U₁₀ that slackens the dependence of U_eff on U₁₀.

Figure 4 Relationship between the effective and local 10 m wind speeds in the IME method, characterized with LES training plumes assuming 1 % instrument precision. Each point represents a different LES plume pseudo-observation from the training set. The red line fits the data to a logarithmic dependence. The 1:1 line is shown in black. See text for similar results with 3 % or 5 % instrument precision.

Table 1 Error standard deviations for retrieving point source rates from column plume observations.

^a Sum of absolute and relative errors when local measurements of 10 m wind speed U₁₀ are available (see text). ^b Additional relative error when local wind speed data are not available. The values given here are inferred from a sample of the GEOS-FP global database over the US and should only be viewed as illustrative. The range is for GEOS-FP wind speeds of 2–7 m s⁻¹, with the larger errors for the smaller wind speeds.

Figure 5 Flow chart describing the IME retrieval algorithm. Algorithm inputs are shown in green, operations in grey, and output in blue. There are two possible paths depending on the availability of 10 m wind speed data: (a) local high-frequency wind speed measurements at the location of the source (right branch), and (b) a temporally averaged meteorological database (left branch).

5.3 Computing the source rate

Figure 5 summarizes the algorithm for retrieving source rates with the IME method. The algorithm accepts two inputs: (1) a map of plume enhancements ΔΩ(x,y) over the plume mask, and (2) the 10 m wind speed U₁₀ from either local high-frequency measurements or an operational meteorological database. If local high-frequency measurements of U₁₀ are available, then there is an opportunity to iteratively refine the plume lifetime τ over which U₁₀ should be averaged, and for this we make a first guess τ₀=5 min. If only coarse-resolution wind speed data are available, then we assume that these are representative of the local value averaged over the plume lifetime and add the associated error to the overall error budget (see Sect. 7).

Figure 6 Evaluation of the IME method for retrieving source rates Q using the LES test set with three different instrument precisions (1 %, 3 %, 5 %). The inset gives the coefficient of determination, R².

Figure 7 Relationship between the effective and local 10 m wind speeds in the cross-sectional flux method, characterized with LES training plumes assuming 1 % instrument precision. Each point represents a different LES plume pseudo-observation from the training set. The red line fits the data to a linear function, excluding the lowest wind speed population (U₁₀<2 m s⁻¹ and U_eff<2 m s⁻¹). See text for corresponding results with 3 % and 5 % instrument precision. U_eff for the cross-sectional flux method is different than for the IME method (Fig. 4).

Figure 6 shows the instantaneous source rates retrieved from our IME algorithm when applied to the test set of LES plumes in different instrument precision scenarios, compared to the true source rates. There is good agreement with the 1:1 line in all cases (R²≥0.86). Retrieval uncertainty (expressed as absolute and relative contributions and defined by the standard deviation of departure from the 1:1 line) increases from 0.07 t h${}^{-\mathrm{1}}+\mathrm{5}$ % for σ=1 %, to 0.13 t h${}^{-\mathrm{1}}+\mathrm{7}$ % for σ=3 %, and 0.17 t h${}^{-\mathrm{1}}+\mathrm{12}$ % for σ=4 % (Table 1). For sources 1.5 t h⁻¹ or larger, retrieval error is less than 25 % even with instrument precision up to 5 %. For Q=1 t h⁻¹, instrument precision up to σ=3 % yields uncertainty less than 20 % of the true source rate. Source rates larger than 0.3 t h⁻¹ contribute more than 75 % of total methane emitted from point sources reporting to the US Greenhouse Gas Reporting Program (Jacob et al., 2016). An instrument with σ=1 % measurement uncertainty can quantify these emissions to within 30 % of the true source rate.

Figure 8 Evaluation of the cross-sectional flux method for retrieving source rates Q using the LES test set with three different instrument precisions (1 %, 3 %, 5 %). The inset gives the coefficient of determination, R².

Figure 9 Evaluation of 10 m wind speeds from the 3 h GEOS-FP global database when used as estimate of local wind speed for source rate calculations in the IME and cross-sectional flux methods. The figure compares 3 h average 10 m wind speeds from the MesoWest database measured at 10 US airports (ABQ, ATL, BOS, DFW, LAX, MCI, MSP, PDX, PHL, and PHX) to corresponding values from the GEOS-FP database. The GEOS-FP data have been corrected for a local roughness height $z_{\mathrm{0},\mathrm{m}}=\mathrm{0.025}$ m (see text). The data are for daytime June 2017 (15:00–21:00 UTC). The fit to a reduced major axis (RMA) regression line is also shown, which closely overlaps the 1:1 line.

6 Computing the source rate by the cross-sectional flux method

Much of our analysis of the IME method in Sect. 5 can be applied to the cross-sectional flux method commonly used for in situ aircraft observations and extended by Krings et al. (2011, 2013) and Tratt et al. (2011, 2014) for remote sensing observations. We compute the plume mask as described in Sect. 5.1, and infer the wind direction from the axis of the plume, based on a weighted average of plume pixel coordinates using the column enhancements as weights. From there, we obtain the mean cross-plume integral C of the column enhancements at different distances downwind of the source (see Eq. 6).

We again use the LES training set to characterize the relationship between the effective wind speed U_eff in Eq. (6) and the local 10 m wind speed U₁₀. U_eff in the cross-sectional flux method is different than U_eff in the IME method. For each plume in the training set, U_eff is computed from Eq. (6) based on C and the known source rate Q. The plume lifetime over which to average local high-frequency U₁₀ measurements for comparison with U_eff is computed as $\mathit{\tau }=L/U_{\mathrm{eff}}$, where the plume size parameter L now has a specific physical meaning as the maximum along-wind distance from the source over which transects can be computed (as defined by the plume mask).

Figure 7 shows the resulting relationship between U_eff and U₁₀. The relationship is near-linear, as would be expected, and the fit U_eff=βU₁₀ with β=1.3–1.5 (where the range is for the 1 %–5 % range of instrument precisions) captures 25 %–75 % of the variance ($\mathrm{0.25}\le R^{\mathrm{2}}\le \mathrm{0.75}$) for U₁₀≥2 m s⁻¹, depending on instrument precision. The 30 %–50 % increase relative to U₁₀ reflects the increase in wind speed with the altitude at which the plume is transported. The departure from the linear relationship for U₁₀<2 m s⁻¹ is because low winds are more variable in direction. The cross-sectional flux method should not be used under calm-wind conditions.

Figure 8 shows the results of the cross-sectional flux retrieval algorithm applied to the LES test plumes, excluding those from the plume population with U₁₀<2 m s⁻¹ and U_eff<2 m s⁻¹. In all instrument precision scenarios, the retrieved source rates are consistent with the 1:1 line. However, residuals are slightly larger than in the IME method (see Fig. 6), as indicated by the smaller coefficients of determination. This results primarily from greater uncertainty in the effective wind speed compared to the IME method. Moreover, analyzing orthogonal plume cross sections requires estimation of the wind direction, which introduces an additional source of error. Absolute and relative retrieval errors estimated in the same way as for the IME method are listed in Table 1. While retrieval uncertainty is slightly higher (0.07–0.26 t h${}^{-\mathrm{1}}+\mathrm{8}$ %–12 %, depending on instrument precision), an advantage of the cross-sectional flux method is that there is a simpler physical basis for relating U₁₀, C, and Q.

7 Inferring the effective wind speed from meteorological databases

Both the IME and cross-sectional flux methods require knowledge of the local wind speed. In the absence of local wind speed measurements, the 10 m wind speed U₁₀ at the time of observation must be estimated from some meteorological database. Here we examine the option of using the GEOS-FP operational reanalysis produced by the NASA Global Modelling and Assimilation Office, available globally as 3 h averages with ${\mathrm{0.25}}^{{}^{\circ }}\times {\mathrm{0.3125}}^{{}^{\circ }}$ resolution ($\approx \mathrm{25}\times \mathrm{25}$ km²) at a lowest grid point level of 60 m above the surface (Molod et al., 2012; https://gmao.gsfc.nasa.gov/GMAO_products/, last access: 24 October 2017). The 10 m wind speed can be obtained from the 60 m wind speed by

$$\begin{array}{}\text{(13)}& U_{\mathrm{10}}=\left[{\displaystyle \frac{\ln \left(\frac{z_{\mathrm{10}}}{z_{\mathrm{0},\mathrm{m}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}}{\ln \left(\frac{z_{\mathrm{60}}}{z_{\mathrm{0},\mathrm{m}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}}}\right]U_{\mathrm{60}},\end{array}$$

where z_0,m [m] is the surface roughness length for momentum, z₁₀=10 m, z₆₀=60 m, and ${\mathrm{\Psi }}_{\mathrm{m}}=f(z/l)$ is a stability correction parameter dependent on the Monin–Obukhov length l (Brasseur and Jacob, 2017). The GEOS-FP data include values for z_0,m and l, but one can use local estimates of these variables if better information is available. Better databases than GEOS-FP may be available to the user depending on region, but an advantage of GEOS-FP is that it can be used as a global default.

Figure 9 evaluates the GEOS-FP U₁₀ data by comparison to 3 h average daytime measurements in June 2017 at 10 US airports obtained from the University of Utah MesoWest database (Horel et al., 2002; MesoWest database, 2018). Here we use $z_{\mathrm{0},\mathrm{m}}=\mathrm{0.025}$ m as input to Eq. (13) to account for the relatively smooth airport terrain. There is no bias in the GEOS-FP data relative to MesoWest. The error standard deviation derived from the difference between the 3 h GEOS-FP and MesoWest 10 m wind speeds is 1.6 m s⁻¹, largely independent of wind speed. Since wind speed is a positive variable, errors at low wind speeds (<2 m s⁻¹) tend to be systematic. There is additional error from using 3 h wind data when the plume lifetime τ is much shorter. From the 5 min resolution of the MesoWest data we find an additional error standard deviation of 2.0 m s⁻¹ for τ=5 min and 1.3 m s⁻¹ for τ=1 h when 3 h average wind speed data are used. Adding these errors in quadrature, we conclude that using GEOS-FP wind data incurs an error standard deviation on the 10 m wind speed of 2.5 m s⁻¹ for small plumes (τ=5 min) and 2.0 m s⁻¹ for large plumes (τ=1 h).

Substitution into the U_eff=f(U₁₀) relations of the IME and cross-sectional flux methods implies an additional error in inferring Q of 15 %–50 % for the IME method and 30 %–65 % for the cross-sectional flux method over the 10 m wind speed range 2–7 m s⁻¹, with largest errors at low wind speeds. The error is larger for the cross-sectional flux method where the dependence of U_eff on U₁₀ is linear rather than logarithmic. Comparison to the other retrieval errors for each method is given in Table 1. At low wind speeds, the error from using GEOS-FP wind data may dominate the overall error budget for inferring source rates. However, our estimate of the error from using operational meteorological databases is intended only to be illustrative. Different errors may apply for other regions or seasons, or when using other meteorological databases than GEOS-FP.

8 Conclusions

We have developed new algorithms for quantifying methane point sources from fine-scale satellite observations of atmospheric column plumes, motivated by the planned fleet of GHGSat instruments ($\le \mathrm{50}\times \mathrm{50}$ m² pixel resolution, 1 %–5 % precision). A challenge is that individual point sources of methane are relatively weak, so that the detectable instantaneous plumes are relatively small (∼1 km) and short lived (<1 h). Using a large ensemble of WRF large eddy simulations (LES) of methane plumes from point sources, we showed that Gaussian plume inversions are unsuccessful because the instantaneous plumes are too small to follow Gaussian behaviour. We also showed how a simple source pixel mass balance method is inappropriate because of wind variability and horizontal turbulent diffusion on the scales of relevance.

Two more promising methods for quantifying source rates from methane column plume observations are the integrated mass enhancement (IME) method and the cross-sectional flux method. Both methods require construction of a plume mask to isolate the plume enhancements from the background noise. The IME method requires estimation of the plume lifetime τ, which in turn depends on an effective wind speed U_eff for the plume and a characteristic plume size L. We showed how these quantities can be estimated from knowledge of the plume mask and of the 10 m wind speed U₁₀ at the location of the source. The source rates are then inferred from the plume observations with expected errors of 0.07–0.17 t h${}^{-\mathrm{1}}+\mathrm{5}$ %–12 % depending on instrument precision (1 %–5 %). For reference, source rates larger than 0.3 t h⁻¹ contribute more than 75 % of total point source emissions in the US Greenhouse Gas Reporting Program (GHGRP) database.

The cross-sectional flux method requires an estimate of the wind direction and of an effective wind speed U_eff reflecting the vertical and horizontal spread of the plume. Again, the LES simulations show how these can be reliably estimated from the plume mask and local U₁₀. We find that for U₁₀≥2 m s⁻¹, U_eff=βU₁₀ with β=1.3–1.5 is a good approximation that accounts for vertical plume spreading. The cross-sectional flux method should not be used for U₁₀<2 m s⁻¹. The errors on the source rates are 0.07–0.26 t h${}^{-\mathrm{1}}+\mathrm{8}$ %–12 %, slightly worse than in the IME method. An advantage of the cross-sectional flux method is its simpler physical basis.

Both the IME and the cross-sectional flux methods parameterize their effective wind speeds U_eff as a function of the local wind speed U₁₀. If local measurements of U₁₀ are not available, then U₁₀ must be estimated from an operational meteorological database or from measurements some distance away. Using the global NASA GEOS-FP archive of wind speeds in June 2017 as an illustrative example compared to US airport data, we find that using this archive would incur source rate errors of 15 %–50 % in the IME method and 30 %–65 % in the cross-sectional flux method over the 2–7 m s⁻¹ range of wind speeds. The largest errors are at low wind speeds where they dominate the overall error budget. Low wind speeds facilitate source detection by improving signal to noise, but worsen source quantification by increasing the uncertainty in the inference of U_eff.

Our source rate retrieval algorithms were motivated by the need to interpret GHGSat plume observations but can readily be applied to any fine-pixel remote sensing measurements of methane column plumes from satellite or aircraft. The precision in retrieving point source rates is much better than can be achieved by current or planned imaging satellites from governmental space agencies, which have higher instrument precision but coarser pixel resolution (Jacob et al., 2016).

Several questions remain to be explored. (1) How will correlated errors in the retrieved methane columns, as observed in GHGSat-D columns (Germain et al., 2017; McKeever et al., 2017), affect source quantification by the IME and cross-sectional flux methods? Such errors would complicate plume definition and plume enhancement calculations. More advanced image segmentation techniques based on machine learning experiments may be useful to differentiate plumes from correlated background errors. (2) How reliably can we parameterize the relationship between effective and 10 m wind speeds given the range of topographies, source elevations, and meteorological environments observed in the real world? Targeted LES experiments may be needed to better constrain the U_eff–U₁₀ relationship for sources in complex topography. (3) How will scattering uncertainties in the photon light paths influence mass enhancement estimates? Clouds in the scene can introduce scattering errors while also masking out portions of the plume. The masked pixels could be estimated by interpolation but scattering errors may be too severe.