Assessing the capability of different satellite observing configurations to resolve the distribution of methane emissions at kilometer scales

Anthropogenic methane emissions originate from a large number of fine-scale and often transient point sources. Satellite observations of atmospheric methane columns are an attractive approach for monitoring these emissions but have limitations from instrument precision, pixel resolution, and measurement frequency. Dense observations will soon be available in both low-Earth and geostationary orbits, but the extent to which they can provide fine-scale information on methane sources has yet to be explored. Here we present an observation system simulation experiment (OSSE) to assess the capabilities of different satellite observing system configurations. We conduct a 1-week WRF-STILT simulation to generate methane column footprints at 1.3 × 1.3 km² spatial resolution and hourly temporal resolution over a 290 × 235 km² domain in the Barnett Shale, a major oil and gas field in Texas with a large number of point sources. We sub-sample these footprints to match the observing characteristics of the recently launched TROPOMI instrument (7 × 7 km² pixels, 11 ppb precision, daily frequency), the planned GeoCARB instrument (2.7 × 3.0 km² pixels, 4 ppb precision, nominal twice-daily frequency), and other proposed observing configurations. The information content of the various observing systems is evaluated using the Fisher information matrix and its eigenvalues. We find that a week of TROPOMI observations should provide information on temporally invariant emissions at ∼ 30 km spatial resolution. GeoCARB should provide information available on temporally invariant emissions ∼ 2–7 km spatial resolution depending on sampling frequency (hourly to daily). Improvements to the instrument precision yield greater increases in information content than improved sampling frequency. A precision better than 6 ppb is critical for GeoCARB to achieve fine resolution of emissions. Transient emissions would be missed with either TROPOMI or GeoCARB. An aspirational high-resolution geostationary instrument with 1.3 × 1.3 km² pixel resolution, hourly return time, and 1 ppb precision would effectively constrain the temporally invariant emissions in the Barnett Shale at the kilometer scale and provide some information on hourly variability of sources.

1 Introduction

Methane is a greenhouse gas emitted by a range of natural and anthropogenic sources (Kirschke et al., 2013; Saunois et al., 2016; Turner et al., 2017). Anthropogenic methane emissions are difficult to quantify because they tend to originate from a large number of potentially transient point sources such as livestock operations, oil or gas leaks, landfills, and coal mine ventilation. Atmospheric methane observations from surface and aircraft have been used to quantify emissions (e.g., Miller et al., 2013; Caulton et al., 2014; Karion et al., 2013, 2015; Lavoie et al., 2015; Conley et al., 2016; Peischl et al., 2015, 2016; Houweling et al., 2016) but are limited in spatial and temporal coverage. Satellite measurements have dense and continuous coverage but limitations from observational errors and pixel resolution need to be understood. Here we perform an observing system simulation experiment (OSSE) to investigate the information content of different configurations of satellite instruments for observing fine-scale and transient methane sources, taking as a test case the oil and gas production sector.

Low-Earth orbit satellite observations of methane by solar backscatter in the shortwave infrared (SWIR) have been available since 2003 from the SCIAMACHY instrument (2003–2012; Frankenberg et al., 2005) and from the GOSAT instrument (2009–present; Kuze et al., 2009, 2016). SWIR instruments measure the atmospheric column of methane with near-unit sensitivity throughout the troposphere. SCIAMACHY and GOSAT demonstrated the capability for high-precision (<1 %) measurements of methane from space (Buchwitz et al., 2015), but SCIAMACHY had coarse pixels (30 × 60 km² in nadir) and GOSAT has sparse coverage (10 km diameter pixels separated by 250 km). Inverse analyses have used observations from these satellite-based instruments to estimate methane emissions at ∼ 100–1000 km spatial resolution  (e.g., Bergamaschi et al., 2009, 2013; Fraser et al., 2013; Monteil et al., 2013; Wecht et al., 2014a; Cressot et al., 2014; Kort et al., 2014; Turner et al., 2015, 2016a; Alexe et al., 2015; Tan et al., 2016; Buchwitz et al., 2017; Sheng et al., 2018a, b). But such coarse resolution makes it difficult to resolve individual source types because of spatial overlap (Maasakkers et al., 2016).

Improved observations of methane from space are expected in the near future (Jacob et al., 2016). The TROPOMI instrument (Veefkind et al., 2012; Butz et al., 2012; Hu et al., 2016, 2018), launched in October 2017, will provide global mapping at 7 × 7 km² nadir resolution once per day. The GeoCARB geostationary instrument (Polonsky et al., 2014; O'Brien et al., 2016) will be launched in the early 2020s with current design values of 3 × 3 km² pixel resolution and twice-daily return time. Additional instruments are presently in the proposal stage with improved combinations of pixel resolution, return time, and instrument precision (Fishman et al., 2012; Butz et al., 2015; Xi et al., 2015).

An OSSE simulates the atmosphere as it would be observed by an instrument with a given observing configuration and error specification. Several OSSEs have been conducted to evaluate the potential of satellite observations to quantify methane sources, but they have been conducted at either coarse (∼ 50 × 50 km²) spatial resolution (Wecht et al., 2014b; Bousserez et al., 2016) or assumed idealized flow conditions (Bovensmann et al., 2010; Rayner et al., 2014). Here we use a 1-week simulation of atmospheric methane with 1.3 × 1.3 km² resolution over a 290 × 235 km² domain to simulate continuous and transient emissions in the Barnett Shale region of Texas, and from there we quantify the capability of different satellite instrument configurations to resolve and quantify these sources at the kilometer and hourly scales. Our choice of scales is guided by the resolution of the planned satellite observations, and our choice of the Barnett Shale is guided by the availability of a high-resolution emission inventory for the region (Lyon et al., 2015). The pattern and density of methane emissions in the Barnett Shale is typical of other source regions in the US (Maasakkers et al., 2016).

2 High-resolution OSSE environment

Figure 1 High-resolution OSSE domain. Left panel shows the successive nested WRF domains at 36, 12, 4, and 1.3 km spatial resolutions, with the coarser domains providing initial and boundary conditions for the finer domains. Black shaded region is the Barnett Shale region in Texas. Right panel shows the domain for the OSSE. Green box is the innermost 1.3 km WRF domain, dashed orange box is the observation domain, and solid orange box is the domain over which the footprints are computed. Light blue lines indicate the counties in the Barnett Shale.

We simulate atmospheric methane concentrations over the Barnett Shale in Texas at 1.3 × 1.3 km² horizontal resolution for the period of 19–25 October 2013 using a framework similar to that of Turner et al. (2016b). The simulation uses version 3.5 of the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) over a succession of nested domains (left panel in Fig. 1) with 1.3 × 1.3 km² spatial resolution in the innermost domain covering 290 × 235 km². There are 50 vertical layers up to 100 hPa. Boundary-layer physics are represented with the Mellor–Yamada–Janjíc scheme and the land surface is represented with the five-layer slab model (Skamarock et al., 2008). The simulation is initialized with assimilated meteorological observations from the North American Regional Reanalysis (https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-regional-reanalysis-narr, last access: 4 May 2018). Overlapping 30 h forecasts were initialized every 24 h at 00:00 UTC and the first 6 h of each forecast were discarded to allow for model spinup. Grid nudging was used in the outermost domain.

WRF meteorology is used to drive the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003). STILT is a Lagrangian particle dispersion model. It advects an ensemble of particles backward in time from selected receptor locations, using the archived hourly WRF wind fields and boundary-layer heights. STILT calculates the footprint for the receptors; a spatiotemporal map of the sensitivity of observations to emissions contributing to the concentration at each selected receptor location and time. We use STILT to calculate 10-day footprints for hourly column concentrations at 1.3 × 1.3 km² resolution over a 70 × 70 km² domain in the innermost WRF nest, tracking the resulting footprints over a 290 × 235 km² domain (right panel in Fig. 1). With this system we examine the constraints on emissions over the 290 × 235 km² domain provided by dense SWIR satellite observations (over the 70 × 70 km² domain) that have up to 1.3 km pixel resolution and hourly daytime frequency. Footprints for each column are obtained by releasing 100 STILT particles from vertical levels centered at 28, 97, 190, and 300 m above the surface and 8 additional levels up to 14 km altitude spaced evenly on a pressure grid. The column footprints are then constructed by summing the pressure-weighted contributions from individual levels, using a typical SWIR averaging kernel taken from Worden et al. (2015) with near-uniformity in the troposphere and correcting for water vapor (see Appendix A in O'Dell et al., 2012).

The footprint for the i^th receptor location and time can be expressed as a vector ${\mathit{h}}_i={\left(\partial y_i/\partial \mathit{x}\right)}^T$ describing the sensitivity of the column concentration y at that receptor location and time to the emission fluxes x over the 290 × 235 km² domain and previous times extending up to 10 days. Here x is arranged as a vector of length n assembling all the emission grid cells and hours, allowing the emissions to vary on an hourly basis. The column concentration is expressed as the dry air column-average mixing ratio (ppb) following common practice (Jacob et al., 2016). The emissions x have units of nmol m⁻² s⁻¹, so that the footprint has units of ppb nmol⁻¹ m² s. The column concentration for the i^th observation (y_i) can be reconstructed from its footprint as

$$\begin{array}{}\text{(1)}& y_i={\mathit{h}}_i\mathit{x}+b_i,\end{array}$$

where b_i is the background column concentration upwind of the 290 × 235 km² domain. We can then write the full set of observations as a vector y of length m and reshape the set of m footprint vectors h into an m×n sparse matrix $\mathbf{H}=\partial \mathit{y}/\partial \mathit{x}$ (where m is the number of observations and n is the number of state vector elements):

$$\begin{array}{}\text{(2)}& \mathit{y}=\mathbf{H}\mathit{x}+\mathit{b},\end{array}$$

where b is the background vector with elements b_i and H is the Jacobian matrix that maps emissions to concentration enhancements due to emissions within our domain.

Figure 2 Summed methane column footprints for all 1.3 × 1.3 km² grid cells in the 70 × 70 km² observation domain defined by the dashed orange box. The footprints are calculated from 08:00 to 17:00 LT over the 290 × 235 km² domain defined by the solid orange box. Bottom right panel shows the summed footprint for the full week.

Figure 2 shows the sum of all column footprints produced on individual days for the 70 × 70 km² observation domain. Computing these high-resolution footprints was a non-trivial computational task and ultimately yielded more than 4 Tb of footprints for the week of pseudo-satellite observations in the Barnett Shale. The footprints show large variability from day to day over the course of the week, reflecting meteorological variability. For example, winds are from the north on 19 October and from the south on 20 October. The winds are weak on 24 October, resulting in a strong local contribution to the footprint. Summing the footprints over the course of the week (bottom right panel of Fig. 2), we find that the observations are mainly sensitive to the core 70 × 70 km² domain where they are made, with a diffuse sensitivity over the outer 290 × 235 km² domain. Additional observations within the outer domain would need to be considered to constrain emissions in that domain. However, information on emissions in the 70 × 70 km² core domain is mainly contributed by observations within the domain. Thus our focus will be to determine the capability of the observations in the 70 × 70 km² domain to constrain emissions within that same domain, but we include the outer 290 × 235 km² domain in our footprint analysis for completeness in accounting of information. Previous work Turner et al. (2016b, Supplement Sect. 6.1) investigated the impact of domain size on error reduction for WRF-STILT inversions in California's Bay Area and found that it had a negligible impact

Figure 3 Gridded Environmental Defense Fund (EDF) methane emission inventory for the Barnett Shale in Texas in October 2013 (Lyon et al., 2015). Spatial resolution is 4 × 4 km². White areas are outside the inventory domain.

Figure 4 Simulated methane concentration enhancements relative to background (ΔXCH₄=Hx) in the 70 × 70 km² observation domain of the Barnett Shale (dashed orange box), as derived from the downscaled EDF methane inventory (x) and the WRF-STILT footprints (H) within the 290 × 235 km² OSSE domain (solid orange box). Values are for 23 October at 09:00 LT. Zeros are due to missing data because of unfinished computations.

The footprint information can be combined with an emission inventory for the 290 × 235 km² domain to generate a field of column concentrations over the 70 × 70 km² domain as would be observed from satellite. For this purpose we use the Environmental Defense Fund (EDF) inventory for the Barnett Shale in October 2013 at 4 × 4 km² resolution compiled by Lyon et al. (2015). We downscale the EDF inventory by uniform attribution from 4 × 4 km² to 1.3 × 1.3 km² spatial resolution. The inventory is shown in Fig. 3 and includes contributions from oil and gas production, livestock operations, landfills, and urban emissions from the Dallas–Fort Worth area. It provides mean monthly values with no temporal resolution, but presumes that some sources will behave as sporadic large transients (Zavala-Araiza et al., 2015). Figure 4 shows an example of the methane column enhancements above background (Hx) computed at 09:00 local time on 23 October. We find enhancements in the range of 0–10 ppb due to emissions within the 290 × 235 km² OSSE footprint domain. In what follows we will examine the potential of different satellite observing systems to detect these enhancements relative to the background and interpret them in terms of local sources.

3 Information content of different satellite observing systems

Table 1 Satellite observing systems considered in this work.

^a Hourly observations are 10 times per day at 08:00–17:00 LT, twice-daily observations are at 10:00 and 14:00 LT, and daily observations are at 13:00 LT.
^b Aspirational instrument with the highest observation frequency and pixel resolution that can be simulated within our OSSE framework.

We aim to determine the information content from different satellite-based observing systems regarding the spatial and temporal distribution of emissions in the Barnett Shale. We consider both steady and potentially transient emissions with five different satellite observing configurations (Table 1). TROPOMI (global daily mapping, 7 × 7 km² nadir pixel resolution, 11 ppb precision; Veefkind et al., 2012) was launched in October 2017 and is expected to provide an operational data stream by the end of 2018. GeoCARB (geostationary, 2.7 × 3.0 km² pixel resolution, 4 ppb precision; O'Brien et al., 2016) is planned for launch in the early 2020s and its observation schedule is still under discussion with a tentative design for observations twice daily; here we examine different return frequencies of hourly, twice daily, and daily. Finally, the hypothetical “hi-res” configuration assumes geostationary hourly observations at the 1.3 × 1.3 km² pixel resolution of our WRF simulation and with 1 ppb precision; it represents an aspirational system that combines the frequent return time, fine pixel resolution, and high precision of instruments presently at the proposal stage (Bovensmann et al., 2010; Fishman et al., 2012; Xi et al., 2015). All configurations are filtered for cloudy scenes.

The various satellite observing configurations of Table 1 differ in their return frequency, pixel resolution, and instrument precision. The benefit of improving any of these attributes may be limited by error in the forward model used in the inverse analysis (i.e., the Jacobian matrix H) and by spatial or temporal correlation of the errors. These limitations are described by the model–data mismatch error covariance matrix (R) including summed contributions from the instrument, forward model, and representation errors (Turner and Jacob, 2015; Brasseur and Jacob, 2017). Representation errors are negligible here because the instrument pixels are commensurate or coarser than the model grid resolution. Instrument error (i.e., precision) is listed in Table 1. Forward model error is estimated by computing STILT footprints for a subset of the meteorological period using the Global Data Assimilation System (GDAS; https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-data-assimilation-system-gdas, last access: 4 May 2018) applying the two sets of footprints to either the EDF methane inventory (Fig. 3; Lyon et al., 2015) or the gridded EPA inventory (Maasakkers et al., 2016), and computing semivariograms of differences in column concentrations. From this we obtain a forward model error standard deviation of 4 ppb with an error correlation length scale of 40 km. We assume a temporal model error correlation length of 2 h. Sheng et al. (2018b) previously derived a temporal model error correlation length of 5 h in simulation of TCCON methane column observations at 25 km resolution, and we expect our correlation length to be shorter because of the finer resolution.

Bayesian inference is commonly used when estimating methane emissions with atmospheric observations, allowing for errors in the observations and in the prior estimates:

$$\begin{array}{}\text{(3)}& P\left(\mathit{x}\right|\mathit{y})\propto P(\mathit{y}\left|\mathit{x}\right)P\left(\mathit{x}\right),\end{array}$$

where P(x|y) is the posterior probability density function (pdf) of the state vector (x) given the observations (y), P(y|x) is the conditional pdf of y given x, and P(x) is the prior pdf of x. A common assumption is that P(y|x) and P(x) are normally distributed which allows us to write the posterior pdf as

$$\begin{array}{ll}\text{(4)}& {\displaystyle }P\left(\mathit{x}\right|\mathit{y})\propto & {\displaystyle }\exp \left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{y}-\mathbf{H}\mathit{x}\right)}^T{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{H}\mathit{x}\right)\right.{\displaystyle }& {\displaystyle }\left.-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}_a\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_a\right)\right\},\end{array}$$

where B is the n×n prior error covariance matrix and x_a is the n×1 vector of prior fluxes. The most probable solution is obtained by minimizing the cost function:

$$\begin{array}{ll}\text{(5)}& {\displaystyle }\mathcal{J}\left(\mathit{x}\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{y}-\mathbf{H}\mathit{x}\right)}^T{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{H}\mathit{x}\right){\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}_a\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_a\right),\end{array}$$

yielding the posterior estimate ($\widehat{\mathit{x}}$):

$$\begin{array}{}\text{(6)}& \widehat{\mathit{x}}={\mathit{x}}_a+\underset{\text{posterior error covariance matrix}}{\underbrace{{\left({\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}+{\mathbf{B}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}}}{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{H}\mathit{x}\right)\end{array}$$

with an n×n posterior error covariance matrix:

$$\begin{array}{}\text{(7)}& \mathbf{Q}=(\underset{\text{observations}}{\underbrace{{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}}}+\underset{\text{prior}}{\underbrace{{\mathbf{B}}^{-\mathrm{1}}}}{)}^{-\mathrm{1}}\end{array}$$

that characterizes the uncertainty in the solution. The first term in the posterior covariance matrix is known as the Fisher information matrix: $\mathcal{F}={\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}$ (see, for example, Rodgers, 2000; Tarantola, 2004).

Comparison between 𝓕 and B⁻¹ identifies the extent to which the observations reduce the uncertainty in the fluxes. Specifically, the number of pieces of information on emissions acquired to better than measurement error is the number of eigenvalues of ${\mathbf{B}}^{\mathrm{1}/\mathrm{2}}\mathcal{F}{\mathbf{B}}^{\mathrm{1}/\mathrm{2}}$ that are greater than unity (Rodgers, 2000). As such, the Fisher information matrix and prior error covariance matrix can quantify the effective rank of the observing system.

A drawback with this formulation of the information content is that it relies on the assumption of a Gaussian prior pdf. A number of papers have suggested that the pdf of methane emissions from a given source may be skewed, with a “fat tail” of transient high emissions (e.g., Brandt et al., 2014; Zavala-Araiza et al., 2015; Frankenberg et al., 2016). Alternate formulations for the cost function to be minimized may include no prior information (least-squares regression), a prior constraint that promotes a sparse solution (e.g., Candes and Wakin, 2008), a prior constraint based on frequentist regularization approaches (such as LASSO regression or Tikhonov regularization), or a prior constraint based on the spatial patterns of emissions rather than their magnitudes (geostatistical inversion). Table 2 lists the corresponding formulations. From Table 2 we see that the observation term is the same in all cases. Thus the Fisher information matrix provides a general measure of the information content provided by an observing system, independent of the form of the prior constraint, and we use it in what follows as a measure of the information content.

Table 2 Cost functions for different formulations of the inverse problem^a.

^a γ is the regularization parameter for LASSO regression and Tikhonov regularization. G is a matrix with columns corresponding to different spatial datasets and β is a vector of drift coefficients for the spatial datasets. Other variables defined in the text.

The Fisher information matrix is an n×n matrix. Each of its n eigenvectors represent an independent normalized emission flux pattern and the corresponding eigenvalues are the inverses of the error variances associated with that pattern. A more useful way of stating this is that the inverse square root of the i^th eigenvalue of 𝓕 represents the flux threshold f_i needed for the observations to be able to constrain the emission flux pattern represented by the i^th eigenvector. Whether that flux threshold is useful depends on the magnitude of the emissions, and this can be assessed for the problem at hand. Thus the eigenanalysis of the Fisher information matrix gives us a general estimate of the capability of an observing system to quantify emissions, which can then be applied to any actual n×n emission field.

For a given emission field, we may expect that some of the n emission flux patterns will be usefully constrained by the observing system while others are not. The number of patterns that are usefully constrained represents the number ℐ≤n pieces of information on emissions provided by the observing system. We will equivalently refer to it as the rank of the Fisher information matrix. This is determined by comparing the eigenvalues of an emission inventory (e_i) to the flux thresholds. The number of e_i larger than the corresponding f_i provides a cut-off to estimate ℐ:

$$\begin{array}{}\text{(8)}& \mathcal{I}=\sum _i\left\{\begin{array}{ll}\mathrm{1},& e_i>f_i\\ \mathrm{0},& e_i\le f_i\end{array}\right.\end{array}$$

In the case of Bayesian inference, this is roughly equivalent to the degrees of freedom for signal with a diagonal prior error covariance matrix and a relative uncertainty of 100 %. But the eigenanalysis of the Fisher information matrix provides a more general approach of the capability of an observing system that can be confronted to any prior constraint and allows intercomparison of different observing system configurations.

There is an inconsistency in this formulation of ℐ: 𝓕 and B⁻¹ have different eigenspaces. In this work we have chosen to treat these matrices separately because, in practice, it is computationally infeasible to directly compute the eigenvalues of the matrix product if n is large, as in the case here of constraining hourly emissions of the spatially distributed inventory. This inconsistency results in our estimate of ℐ likely being an upper bound on the information content (see Appendix for details).

4 Comparing different satellite configurations

The eigenanalysis of Sect. 3 allows us to intercompare the value of different satellite configurations for resolving the fine-scale patterns of methane emissions within a given domain. Here we apply it to the Barnett Shale domain of Sect. 2. We consider two limiting cases: Case 1 assumes the emissions to be temporally invariant and Case 2 assumes the emissions to vary hourly with no temporal correlation. In Case 1 the problem is typically overdetermined (m>n), depending on the satellite configuration, and the maximum rank of 𝓕 is n (the number of emission grid cells). In Case 2 the problem is underdetermined (m<n) and the maximum rank of 𝓕 is m (the number of observations).

In both Case 1 and 2, the observations only provide useful information (as defined by Eq. 8) if the signal is larger than the noise, as diagnosed by the e_i>f_i criterion of Eq. 8. Here the emissions are the downscaled EDF inventory, which includes 40 140 grid cells in the 290 × 235 km² inversion domain (n= 40 140 in Case 1 with temporally invariant emissions) but only 2 601 of those grid cells are within the 70 × 70 km² observation domain (dashed orange box in Fig. 1), where we might expect the observations to provide the strongest constraints. In Case 2 with temporally variable emissions we have $n=\mathrm{40}\mathrm{140}\times \mathrm{24}=\mathrm{963}\mathrm{360}$ grid cells for a single day.

Figure 5 Capability of different configurations for satellite observations of atmospheric methane (Table 1) to resolve the fine-scale (1.3 × 1.3 km²) patterns of variability of temporally invariant emissions in a 290 × 235 km² domain and for a 1-week observation period. The colored lines show the flux thresholds for the different emission patterns of variability in the domain, as given by the ordered inverse square roots of the eigenvalues of the Fisher information matrix. Solid black line is the eigenvalues of the emissions from the EDF Barnett Shale methane inventory (Lyon et al., 2015) and the solid gray line is the gridded EPA inventory. The region above the black line is where the noise is larger than the signal. Filled circles indicate the information content of the observing system (ℐ) for a given satellite configuration at 1.3 × 1.3 km² spatial resolution. Inset table lists the information contents for the five configurations.

Figure 6 Capability of GeoCARB-like satellite configurations to resolve the fine-scale (1.3 × 1.3 km²) patterns of variability of temporally invariant emissions in a 290 × 235 km² domain and for a 1-week observation period. Left panel shows the results for a configuration with 10 returns per day (hourly observations) where the instrument precision is varied from 0 to 14 ppb. Right panel shows the results for a configuration with 4 ppb instrument precision and the return frequency per day is varied from 1 to 10. Solid black line shows eigenvalues of the EDF Barnett Shale methane emission inventory (Lyon et al., 2015). The region above the black line is where the noise is larger than the signal. The change in flux threshold as the sampling frequency increases in right panel is not necessarily monotonic, because some of the cases use different subsets of observation (e.g., daily observations are at 13:00 LT while twice daily are at 10 and 14).

Figure 5 shows the ensemble of flux thresholds for the five satellite configurations, assuming temporally invariant emissions. The ranked flux patterns are on the abscissa; leading flux patterns correspond to larger patterns of variability (e.g., regional-scale emissions), and the trailing flux patterns correspond to fine-scale variability. The corresponding flux thresholds are on the ordinate. The flux threshold is lowest for the leading flux patterns and largest for the trailing flux patterns. This means that the regional-scale emissions are easiest to quantify and the finer-scale emissions are increasingly difficult to quantify. The information content (ℐ) is obtained from the intersection of the flux thresholds (colored lines) with the eigenvalues from the emission inventory (black line). A higher information content means that finer scales of emission variability can be detected.

From Fig. 5, we see that a week of TROPOMI observations provides five pieces of information on emissions for the 70 × 70 km² core domain out of a possible 2601 pieces of information describing the emissions on the 1.3 × 1.3 km² grid. The actual pieces of information are the eigenvectors of the Fisher information matrix, and the ranked eigenvectors describe gradually finer patterns of variability from 70 × 70 to 1.3 × 1.3 km². The k^th ranked eigenvector may be assumed to describe an emission pattern of dimension 70$/\sqrt{k}$, implying that TROPOMI can resolve emissions on a 30 km scale.

The three GeoCARB configurations provide 98–961 pieces of information dependent on whether the observations are daily, twice daily, or hourly. Following the above assumption, this corresponds to resolving emissions on a ∼ 2–7 km scale. Hourly observations provide 10 times more information (as defined by Eq. 8) on emission patterns than daily observations, and 3 times more than twice-daily observations (the default configuration of GeoCARB). Remarkably, more is gained by going from daily to twice daily (factor of 3.4) than going from twice daily to hourly (factor of 2.9) because of the temporal error correlation in the transport model. The aspirational hi-res satellite configuration provides 2 221 pieces of information on temporally invariant sources, corresponding to 85 % of the flux patterns in the 70 × 70 km² observation region, which means that much of the spatial variability in the 1.3 × 1.3 km² emissions in the Barnett Shale is resolved.

Figure 6 further quantifies the importance of instrument precision and return frequency for the GeoCARB pixel resolution of 2.7 × 3.0 km². It shows the flux thresholds for a set of configurations where the instrument precision is varied from 0 to 14 ppb and the return frequency is varied from 1 to 10 returns per day. We find that instrument precision is more important than return frequency for increasing the information content from the observations.

Figure 7 Same as Fig. 5 but for temporally variable sources on 21 October 2013.

Figure 8 Same as Fig. 6 but for temporally variable sources on 21 October 2013.

In Case 2 we assume that the methane sources in individual pixels vary in time on an hourly basis with no correlation from one hour to the next, making the problem generally underdetermined (m<n) for all satellite configurations. Here we aim to determine the ability of the satellite observations to quantify the hourly emissions over the spatial patterns defined by the eigenvectors of 𝓕 and making no assumption as to the persistence of those emissions. We treat each day independently and compute the eigenvalues of the Fisher information matrix for each day. Figure 7 shows the flux thresholds for the five satellite configurations on a representative day. From Fig. 7, we see that TROPOMI is unable to provide any information on hourly emissions in the Barnett Shale. The three GeoCARB configurations provide 2–54 pieces of information. Figure 8 evaluates the impact of sampling frequency and instrument precision for the GeoCARB configurations. As with the temporally invariant case, we find that instrument precision is more important for increasing the information content. The aspirational “hi-res” configuration (shown in Fig. 7) is the only configuration that is able to provide substantial information (458 pieces of information) on temporally variable emissions.

Figure 9 Information content ℐ as a function of the instrument precision (left column) and the sampling frequency per day (right column) for a satellite with a pixel resolution of 2.7 × 3.0 km². Top row is for Case 1 where the sources are assumed to be temporally invariant and bottom row is for Case 2 where the sources are temporally variable. Solid black line is the median information content. A 4 ppb model error is included; see Sect. 3. Uncertainty is from randomly sampling e_i from the eigenvalues of the EDF inventory.

Figure 9 summarizes the findings from Figs. 6 and 8. It compares the information content ℐ from configurations with 2.7 × 3.0 km² spatial resolution (GeoCARB) as the instrument precision and return frequency are varied from 0 to 14 ppb and 1 to 10 returns per day, respectively, for both temporally variable and constant sources. Uncertainty of ℐ is estimated by randomly sampling e_i from the ensemble of emission inventory eigenvalues and comparing to f_i in Eq. 8. For the temporally invariant sources (Case 1), we find considerable increases in information content for instrument precisions better than 6 ppb (top left panel in Fig. 9) and an approximately linear relationship between information content and return frequency (top right panel in Fig. 9). The satellite configurations provide considerably less information for the temporally variable sources (Case 2). We find that satellite configurations with instrument precision worse than 6 ppb provide no information on temporally variable sources (bottom left panel in Fig. 9). As with the temporally invariant case, we find an approximately linear relationship between information content and return frequency (bottom right panel in Fig. 9). From this, we conclude that a GeoCARB-like instrument would greatly benefit from having an instrument precision better than 6 ppb.

5 Conclusions

We conducted an observing system simulation experiment (OSSE) to evaluate the potential of different satellite observation systems for atmospheric methane to quantify methane emissions at kilometer scale. This involved a 1-week WRF-STILT simulation of atmospheric methane columns with 1.3 × 1.3 km² spatial resolution over the 290 × 235 km² Barnett Shale domain to quantify the information content of different satellite instrument configurations for resolving the kilometer-scale distribution of methane emissions within that domain. We evaluated the information content of the different satellite observing systems through an eigenanalysis of the Fisher information matrix 𝓕, which characterizes the capability of an observing system independently of the form of the prior information. The eigenvalues of 𝓕 define the emission flux thresholds for detection of emission patterns down to 1.3 km in scale as defined by the eigenvectors. Here we put these flux thresholds in context of the high-resolution EDF emission inventory for the Barnett Shale to quantify the information content from different satellite observing configurations. The same approach could be readily used for different observation domains and different prior inventories.

We find from this analysis that the recently launched TROPOMI satellite instrument (low-Earth orbit, 7 × 7 km² pixels, daily return time, 11 ppb precision) should be able to constrain the mean emissions in the Barnett Shale and provide some coarse-resolution information on the distribution of temporally invariant emissions at ∼ 30 km scales. The planned GeoCARB instrument (geostationary orbit, 2.7 × 3.0 km² pixels, twice-daily return time, 4 ppb precision) will provide 50 times more information than TROPOMI. The observing frequency of GeoCARB is still under discussion; we find that twice-daily observations triple the information content relative to daily observations, while hourly observations allow another tripling. The 4  ppb precision of GeoCARB is well adapted to the magnitude of methane sources; we find that a precision larger than 6 ppb would considerably decrease the information content. An aspirational “hi-res” instrument using attributes of currently proposed instruments (geostationary orbit, 1.3 × 1.3 km² pixels, hourly return time, 1 ppb precision) can resolve much of the kilometer-scale spatial distribution in the EDF inventory. This assumes that the emissions are constant in time or that their temporal variability is known. Resolving hourly variable emissions at the kilometer scale will be very limited even with the aspirational “hi-res” instrument.

Appendix A: Computing the information content

We treat 𝓕 and B⁻¹ separately because it is computationally infeasible to compute the eigenvalues of the matrix product when we attempt to resolve hourly emissions as n>10⁶ and both 𝓕 and B⁻¹ are n×n matrices. This separation of 𝓕 and B⁻¹ results in our estimate of ℐ likely being an upper bound on the information content. This follows from Bhatia (1997) who prove that $\mathit{\lambda }\left(\mathbf{CD}\right){\prec }_w{\mathit{\lambda }}^{\downarrow }\left(\mathbf{C}\right)\cdot {\mathit{\lambda }}^{\downarrow }\left(\mathbf{D}\right)$, where C and D are Hermitian positive definite matrices, λ^↓(X) denotes the vector of eigenvalues of X in decreasing order, ≺_w is the weak majorization preorder, and $\mathbf{p}\times \mathbf{q}=\left(p_{\mathrm{1}}{q}_{\mathrm{1}},\mathrm{\dots },p_n{q}_n\right)$. Therefore, directly computing the eigenvalues of ${\mathbf{B}}^{\mathrm{1}/\mathrm{2}}\mathcal{F}{\mathbf{B}}^{\mathrm{1}/\mathrm{2}}$, as Rodgers (2000) suggests for the Bayesian inference case with Gaussian errors, would likely yield fewer eigenvalues larger than unity than our estimate.

In the case of temporally variable emissions, the system is generally underdetermined (m<n) and we can use a singular value decomposition to efficiently compute the eigenvalues of 𝓕. For an m×n real matrix A, the non-zero singular values of A^TA and AA^T are identical even though the singular vectors are different (see, for example, Rodgers, 2000) but the dimensions of these two matrices are n×n and m×m, respectively, and the eigenvalues can be computed from the square root of the non-zero singular values. We can write $\mathcal{F}={\widehat{\mathbf{H}}}^T\widehat{\mathbf{H}}$, where $\widehat{\mathbf{H}}={\mathbf{L}}^{-\mathrm{1}}\mathbf{H}$ is the pre-whitened Jacobian and L is a lower triangular matrix from a Cholesky decomposition of R (such that R=LL^T). Thus, the eigenvalues of 𝓕 can be obtained by analysis of either ${\widehat{\mathbf{H}}}^T\widehat{\mathbf{H}}$ (an n×n matrix) or $\widehat{\mathbf{H}}{\widehat{\mathbf{H}}}^T$ (an m×m matrix).