High-resolution inversion of methane emissions in the Southeast US using SEAC⁴RS aircraft observations of atmospheric methane: anthropogenic and wetland sources

We use observations of boundary layer methane from the SEAC⁴RS aircraft campaign over the Southeast US in August–September 2013 to estimate methane emissions in that region through an inverse analysis with up to $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ (25×25 km²) resolution and with full error characterization. The Southeast US is a major source region for methane including large contributions from oil and gas production and wetlands. Our inversion uses state-of-the-art emission inventories as prior estimates, including a gridded version of the anthropogenic EPA Greenhouse Gas Inventory and the mean of the WetCHARTs ensemble for wetlands. Inversion results are independently verified by comparison with surface (NOAA∕ESRL) and column (TCCON) methane observations. Our posterior estimates for the Southeast US are 12.8±0.9 Tg a⁻¹ for anthropogenic sources (no significant change from the gridded EPA inventory) and 9.4±0.8 Tg a⁻¹ for wetlands (27 % decrease from the mean in the WetCHARTs ensemble). The largest source of error in the WetCHARTs wetlands ensemble is the land cover map specification of wetland areal extent. Our results support the accuracy of the EPA anthropogenic inventory on a regional scale but there are significant local discrepancies for oil and gas production fields, suggesting that emission factors are more variable than assumed in the EPA inventory.

1 Introduction

Methane is an important greenhouse gas (Myhre et al., 2013). Individual countries must report their national emissions to the United Nations Framework Convention on Climate Change (UNFCCC; United Nation, 1992). Observations of atmospheric methane reviewed by Brandt et al. (2014) have implied that the US national inventory reported by the Environmental Protection Agency (EPA) may be greatly underestimated. Here we use aircraft observations from the NASA SEAC⁴RS aircraft campaign over the Southeast US (Toon et al., 2016), together with a newly gridded version of the EPA inventory (Maasakkers et al., 2016), in a fine-resolution inversion with detailed error characterization to better quantify methane emissions over this major source region.

The EPA (2016) reports a national anthropogenic emission total of 29.2 Tg CH₄ a⁻¹ for 2014, with no significant trend over the past decade and less than ±3 % interannual variability. Major contributors are livestock (32 %), the oil and gas industry (32 %), waste (22 %), and coal mining (8 %). The EPA (2016) inventory is consistent with Lyon et al. (2015) for oil and gas systems and Wolf et al. (2017) for livestock, and 8 % higher than the previous versions (EPA, 2013, 2014), largely due to updated oil and gas emissions. There is also a highly uncertain natural source from wetlands, estimated at 4.5–14 Tg a⁻¹ for the contiguous US in the WETCHIMP compilation of inventories (Melton et al., 2013). Inverse analyses of atmospheric methane observations have suggested that the EPA bottom-up inventory (EPA, 2013, 2014) is too low by about 30 % (Miller et al., 2013; Turner et al., 2015; Alexe et al., 2015). However, Turner et al. (2015) and Alexe et al. (2015) relied on prior estimates from the global EDGAR v4.2 inventory (European Commission, 2011) that have large errors in source patterns particularly for oil and gas systems (Maasakkers et al., 2016; Sheng et al., 2017). Errors in source patterns used as prior estimates can greatly bias inversion results (Jacob et al., 2016), though this depends on the constraint from observations and on the uncertainty assigned to the prior estimates. Miller et al. (2013) used a geostatistical inversion that did not rely on any prior estimates, but had little constraints in the Southeast US.

The SEAC⁴RS aircraft campaign conducted in August–September 2013 offers an opportunity for better estimating methane emissions in the Southeast US, a region that accounts for about half of total anthropogenic methane emissions in the US, according to the gridded EPA inventory (Maasakkers et al., 2016), and also has extensive wetlands. The aircraft flights provided extensive boundary layer measurements of methane across the region. We conduct an inverse analysis of the SEAC⁴RS data with the GEOS-Chem chemical transport model (CTM) at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ resolution, using state-of-the-art prior estimates from the gridded EPA inventory of Maasakkers et al. (2016) and the WetCHARTs extended ensemble wetlands inventory of Bloom et al. (2017). This allows us to evaluate the EPA inventory with better accuracy and resolution than has been done before and also to gain better understanding of US wetland emissions.

2 Methods

We derive an optimized estimate of spatially resolved methane emissions in the Southeast US (domain of Fig. 1) through Bayesian inverse analysis of atmospheric methane observations from the SEAC⁴RS aircraft campaign. Let the vector x represent a gridded ensemble of methane emissions in the region (state vector for the inversion).The inversion minimizes the cost function J(x) by solving ∇_xJ(x)=0 as follows:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }J\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{x}-{\mathit{x}}_{\mathbf{A}}{)}^T{\mathbf{S}}_{\mathbf{A}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_{\mathbf{A}}){\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{y}-\mathbf{K}\mathit{x}{)}^T{\mathbf{S}}_{\mathbf{O}}^{-\mathrm{1}}(\mathit{y}-\mathbf{K}\mathit{x}).\end{array}$$

Here the methane observations are assembled as a vector y, x_A is the prior emission estimate, K is the Jacobian matrix describing the sensitivity of concentrations to emissions, and S_A and S_O are the prior and observational error covariance matrices, respectively. The observational error includes contributions from both the instrument error and the model transport error.

Figure 1 Methane emissions in the Southeast US in August–September 2013. Panel (a) shows the prior anthropogenic and wetland methane emissions, and panel (b) shows the inversion results including posterior emissions, scaling factors (posterior ∕ prior emission ratios), and the diagonal elements of the averaging kernel matrix representing the sensitivity of the inversion results to the observations. The sum of these diagonal elements over the domain (trace of the averaging kernel matrix) quantifies the degrees of freedom for signal (DOFS) of the inversion. Numbers inset in the emission panels are the regional totals expressed as annual means for clarity (i.e., assuming that August–September emission rates hold for the rest of the year). Values in parentheses are the totals for the region with averaging kernel sensitivities larger than 0.05 (stippled areas in lower left panel).

Analytical solution of ∇_xJ(x)=0 yields the optimal estimate $\widehat{\mathit{x}}$, the posterior error covariance matrix $\widehat{\mathbf{S}}$, and the associated averaging kernel matrix A as follows (Rodgers, 2000; Brasseur and Jacob, 2017):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\widehat{\mathit{x}}={\mathit{x}}_{\mathbf{A}}+{\mathbf{S}}_{\mathbf{A}\mathbf{K}}^T({\mathbf{KS}}_{\mathbf{A}\mathbf{K}}^T+{\mathbf{S}}_{\mathbf{O}}{)}^{-\mathrm{1}}(\mathit{y}-\mathbf{K}{\mathit{x}}_{\mathbf{A}}),\text{(3)}& {\displaystyle }& {\displaystyle }{\widehat{\mathbf{S}}}^{-\mathrm{1}}={\mathbf{K}}^T{{\mathbf{S}}_{\mathbf{O}}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathbf{A}}^{-\mathrm{1}},\text{(4)}& {\displaystyle }& {\displaystyle }\mathbf{A}={\mathbf{I}}_n-\widehat{\mathbf{S}}{\mathbf{S}}_{\mathbf{A}}^{-\mathrm{1}},\end{array}$$

where I_n is the identity matrix with n being the dimension of the state vector x. Inversions of atmospheric methane observations may solve ∇_xJ(x)=0 either analytically, or numerically using an adjoint method (Jacob et al., 2016). Unlike adjoint-based inversions, analytical solution provides direct error characterization of the optimal estimate $\widehat{\mathit{x}}$ through its error covariance matrix $\widehat{\mathbf{S}}$. The related averaging kernel matrix A describes the sensitivity of the optimal estimate $\widehat{\mathit{x}}$ to the true emissions x. The trace of A quantifies the degrees of freedom for signal (DOFS), i.e., the number of pieces of information in the observing system for constraining the methane emissions (DOFS ≤n).

The Jacobian matrix K for the inversion is constructed with the GEOS-Chem CTM (http://www.geos-chem.org, last access: 1 November 2017), which relates methane emissions to atmospheric concentrations through simulation of atmospheric transport. We use a nested version of GEOS-Chem as described by Kim et al. (2015) with $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ horizontal resolution over the North America window and adjacent oceans (9.75–60^∘ N, 130–60^∘ W), driven by GEOS-FP assimilated meteorological data from the NASA Global Modeling and Assimilation Office (GMAO). The same version of the GEOS-Chem has been applied to simulation of other chemical observations from the SEAC⁴RS campaign (Kim et al., 2015; Fisher et al., 2016; Marais et al., 2016; Travis et al., 2016; Zhu et al., 2016; Yu et al., 2016; Chan Miller et al., 2017). The boundary conditions for the nested-grid simulation are from a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ global simulation by Turner et al. (2015) using methane emissions optimized with three years of GOSAT satellite data. The GOSAT-optimized emissions have been independently evaluated with atmospheric methane observations from the NOAA surface network (Turner et al., 2015). The GEOS-Chem uses a 3-D archive of monthly average OH concentrations from Park et al. (2004) to compute the methane sink, with a lifetime of 8.9 years in the troposphere consistent with observational constraints (Prather et al., 2012; Turner et al., 2017). The sink is irrelevant for our North American simulation since ventilation of the domain is much faster (Wecht et al., 2014). Since we treat OH concentrations as decoupled from methane in the inversion, the relationship between emissions and concentrations is strictly linear, so that K fully describes the GEOS-Chem model for the purpose of the inversion.

The prior emission estimates for the inversion are taken from the $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ gridded version of the EPA anthropogenic greenhouse gas emission inventory for 2012 (Maasakkers et al., 2016) and the mean wetland emissions from the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ monthly WetCHARTs extended ensemble for 2013 (Bloom et al., 2017). Figure 1 (top panels) shows the distribution of these prior methane emissions over the inversion domain for August–September 2013. Emissions total 13.3 Tg a⁻¹ for anthropogenic sources and 13.0 Tg a⁻¹ for wetlands over these two months (expressed on an annual basis). Anthropogenic emissions in the EPA inventory have little seasonal or interannual variability (Turner et al., 2015; Maasakkers et al., 2016), while wetland emissions have a large seasonal variation.

The SEAC⁴RS DC-8 aircraft conducted 21 flights over the Southeast between 6 August and 21 September 2013. Methane was measured by gas chromatography from whole air flask samples and calibrated to the NOAA standard. Figure 2 (left panel) shows the SEAC⁴RS flight tracks and the spatial distribution of the methane flask measurements below 2 km altitude. The mean observed vertical profile is shown in the right panel of Fig. 2, and compared to the GEOS-Chem profile using the prior emissions. The model is unbiased in the free troposphere above 2 km, implying a successful representation of background methane by the boundary conditions. Model overestimation in the boundary layer below 2 km suggests that the prior US emissions are too high. For the inversion, we use the SEAC⁴RS observations below 2 km altitude averaged over the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ GEOS-Chem grid for individual flights. This represents 652 individual methane observations.

Figure 2 Boundary layer methane concentrations over the Southeast US measured during the SEAC⁴RS aircraft campaign (6 August–21 September 2013). Panel (a) shows the flight tracks in grey and the methane measurements at 0–2 km altitude. The three NOAA∕ESRL sites at SGP (Southern Great Plains, Oklahoma; 36.6^∘ N, 97.5^∘ W), WKT (Moody, Texas; 31.3^∘ N, 97.3^∘ W), and SCT (Beech Island, South Carolina; 33.4^∘, 81.8^∘ W) are indicated. SGP is co-located with the TCCON site at Lamont, Oklahoma. Panel (b) shows the mean methane vertical profiles over the Southeast US domain measured from the aircraft and simulated by GEOS-Chem using the prior and posterior emissions.

We use the residual error method (Heald et al., 2004) to estimate the diagonal elements of the observational error covariance matrix S_O. The method assumes that the mean model bias between the observations and the model is due to error in prior emissions to be corrected by the inversion. The residuals are the differences between observed and modeled values after removing the mean model bias. The residual error standard deviation (RSD) is taken to represent the observational error including contributions from the instrument and the transport model. Figure 3 shows the vertical profile of the RSD for the ensemble of the SEAC⁴RS data over the Southeast US. The RSD is about 60 ppb below 2 km and 20 ppb in the free troposphere above. Subsetting the data by latitudinal bands gives similar results. Thus we use 60 ppb for the standard deviation of the observational error (diagonal elements in S_O). The instrument precision is better than 2 ppb (Simpson et al., 2002), therefore most of that observational error is from the transport model. We take S_O to be diagonal since error correlations between boundary layer observations on the GEOS-Chem grid are not significant (Wecht et al., 2014).

Figure 3 Residual standard deviation (RSD) of the difference between SEAC⁴RS methane observations and the GEOS-Chem model with prior emissions, for 1 km altitude bins. The RSD defines the observational error standard deviations for the inversion as described in the text. Values are shown for two latitudinal bands.

The inversion can, in principle, optimize emissions at the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ grid resolution of the GEOS-Chem model, representing 3004 grid cells over the inversion domain. But the aircraft observations do not have sufficient information to constrain emissions at that resolution. In order to reduce the dimensionality of the state vector, we project the 3004 grid cells onto 216 elements of a Gaussian mixture model (GMM) with radial basis functions based on spatial proximity and source type patterns (Turner and Jacob, 2015). The use of the GMM allows us to retain high resolution of up to 25 km for major localized sources while degrading resolution in areas of weak or broadly distributed sources. Areas dominated by wetlands have resolution of 100–200 km in the GMM because they are broadly distributed.

The anthropogenic inventory of Maasakkers et al. (2016) and the wetlands inventory of Bloom et al. (2017) both include gridded error estimates that serve as the diagonal elements of the prior error covariance matrix S_A. Maasakkers et al. (2016) found no significant spatial error correlation at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ resolution in their inventory while a variogram analysis across the elements of the WetCHARTS ensemble indicates a spatial error correlation length scale of 130 km. Here we ignore this correlation and take S_A to be diagonal.

3 Results and discussion

Figure 1b shows the results of the inversion including the optimized posterior emissions, the corrections to the prior emissions, and the DOFS as measured by the diagonal elements of the averaging kernel matrix. Figure 4 (top panels) compare the observed boundary layer methane concentrations to the values simulated by GEOS-Chem with prior and posterior emissions (Fig. 1). The simulation with prior emissions has a positive bias that is effectively corrected when using posterior emissions. The coefficient of determination (R²) between model and observations increases from 0.30 to 0.50 when using posterior emissions. Figure 4 also evaluates the SEAC⁴RS inversion results with independent surface air observations from the three NOAA∕ESRL surface network sites in the region (Andrews et al., 2014) and with methane column observations from the TCCON site in Lamont, Oklahoma (Wunch et al., 2011; Wennberg et al., 2017). The posterior emissions improve the simulation of these independent data sets. GOSAT satellite observations are another source of independent data but the 2 month period is too sparse for useful evaluation (Wecht et al., 2014).

Figure 4 Evaluation of the SEAC⁴RS inversion of methane emissions in the Southeast US for the 6 August–21 September 2013 period. Panels (a, b) compare GEOS-Chem methane concentrations with the SEAC⁴RS observations, using prior emissions (left) and posterior emissions (right). Panels (c, d) compare GEOS-Chem methane concentrations with independent observations from the three NOAA∕ESRL surface sites in the inversion domain (see Fig. 2 and caption). Panels (e, f) compare GEOS-Chem methane columns with TCCON hourly column observations at Lamont, Oklahoma (Wennberg et al., 2017), after correction for stratospheric bias in the model (Turner et al., 2015). The 1:1 lines (dashed) and the reduced-major-axis(black solid line) linear regressions are also shown, along with the coefficients of determination (R²) and the slopes (±1σ) derived from the bootstrap method.

Total posterior emissions over the SEAC⁴RS domain are 15 % (4 Tg a⁻¹) lower than the prior estimate (Fig. 1). The inversion is able to constrain about 10 pieces of information in the spatial distribution of methane emissions as measured by the DOFS. It is strongly sensitive to the Gulf Coast and to large anthropogenic source areas such as the Floyd Shale in central Alabama. For the regions with averaging kernel sensitivity larger than 0.05, posterior emissions are 35 % lower than the prior estimate. The posterior errors are 18–30 % over these regions. The scaling factors show large downward corrections of prior emissions in Louisiana and Mississippi, and along the Gulf Coast, where wetlands are the dominant sources. There are also downward corrections in southern West Virginia, where coal mines are dominant, and in the Haynesville Shale gas production region of northern Louisiana and southern Arkansas. On the other hand, there are significant upward corrections for the coal mines of southern Illinois and for the Floyd Shale in central Alabama.

We can attribute the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ scaling factors from the inversion to specific methane source sectors by using the sector-resolved spatial patterns in the prior emission inventories, as described by Turner et al. (2015) but here with the improved anthropogenic source patterns from Maasakkers et al. (2016) and wetland emissions from (Bloom et al., 2017). Anthropogenic and wetland sources are well separated spatially in these inventories. Figure 5 compares our results with the prior emission totals for the different sectors in the Southeast US. We find a significant 27 % reduction in regional wetland emissions relative to the prior estimate (mean of the WetCHARTs extended ensemble). For the subdomain with averaging kernel sensitivity larger than 0.05 that reduction is 42 %. By contrast, we find no significant regional bias in the EPA anthropogenic inventory for any of the major source sectors for the period of August–September 2013. However, there are large local biases that tend to cancel each other on a regional scale (e.g., Haynesville Shale vs. Floyd Shale for natural gas). This suggests that methane emission factors for the oil and gas sector are more variable than assumed in the EPA (2016) inventory.

Figure 5 Prior and posterior methane emissions in the Southeast US (domain of Fig. 1) for August–September 2013. The prior anthropogenic emissions are from the EPA national inventory for 2012 (EPA, 2016; Maasakkers et al., 2016) and the prior wetland emissions are the means of the WetCHARTs extended ensemble (Bloom et al., 2017). Error bars (one standard deviation) on sectoral emissions are from the prior and posterior error variances of our inversion. Methane emissions in the subdomain with averaging kernel sensitivities larger than 0.05 (Fig. 1) are also indicated.

The WetCHARTs extended ensemble includes 18 wetland methane emission models intended to encompass the uncertainties in estimating wetland emissions (Bloom et al., 2017). The different models (ensemble members) use different datasets for wetland extent fraction A [m² wetlands m⁻² surface area], heterotrophic respiration rate R [mg C day⁻¹ m⁻² of wetland area], temperature-dependent factor $q_{\mathrm{10}}^{T/\mathrm{10}}$ of C respired as CH₄ [mg CH₄ mg⁻¹ C] where T is the surface skin temperature, and global scaling factors s. The wetland methane emission flux E [mg CH₄ m⁻² day⁻¹] at a time t and location x for each of these members is given by the following equation:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }E(t,x)=sA(t,x)R(t,x)q_{\mathrm{10}}^{T(t,x)/\mathrm{10}}.\end{array}$$

The 18-member ensemble consists of three temperature dependence factors ($q_{\mathrm{10}}=\mathrm{1},\mathrm{2},\mathrm{3}$), three global scale factors ($s=\mathrm{125},\mathrm{166},\mathrm{208}$), and two wetland extent maps (A) from the Global Lakes and Wetlands Database (GLWD; Lehner and Dölla, 2004) and GLOBCOVER (Bontemps et al., 2011). The heterotrophic respiration rate (R) is the median output from the carbon data model framework (CARDAMOM; Bloom et al., 2016), and is not varied across that ensemble.

Figure 6 shows the Southeast US wetland emissions for each WetCHARTs member, along with the root-mean-square error (RMSE) of its spatial distribution relative to our optimized posterior estimate on the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ grid. Consistency in spatial distribution with our optimized estimate is indicated by a low RMSE. We find that the specification of wetland extent is the most systematic source of error in wetland emission estimates; all GLOBCOVER-based models underestimate wetland emissions, while all GLWD-based models overestimate emissions. Estimates using q₁₀=1 (no temperature dependence in the CH₄ : C respiration ratio) exhibit the lowest RMSE values. The WetCHARTs ensemble mean used as prior for our inversion has the lowest RMSE, although this may reflect its influence on the posterior solution.

Figure 6 Range of wetland emission estimates in the Southeast US (domain of Fig. 1) for August–September 2013. The figure shows the spread of the WetCHARTs extended ensemble and compares with the posterior emission estimate from our inversion in terms of emission total and root-mean-square error (RMSE) on the $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ spatial grid. The RMSE of the posterior emission estimate with itself is zero by definition. WetCHARTs ensemble members use wetland areal extent data from either the GLOBCOVER (Bontemps et al., 2011) or GLWD (Lehner and Dölla, 2004) databases, as well as different estimates of temperature sensitivity q₁₀ and global scaling factors s (see Eq. 5 and text). The posterior wetland emission estimate from our inversion is shown as black solid diamond with its associated error standard deviation. The mean of the WetCHARTs ensemble used as prior for our inversion is shown as blue solid circle with its associated error standard deviation.

4 Conclusions

We used extensive boundary layer methane observations from the SEAC⁴RS aircraft campaign over the Southeast US in August–September 2013 to optimize methane emissions in that region with up to $\mathrm{0.25}{}^{\circ }\times \mathrm{0.3125}{}^{\circ }$ spatial resolution and with detailed error characterization. The inversion used new state-of-the-art inventories as prior information, including the gridded version of the EPA (2016) national anthropogenic inventory from Maasakkers et al. (2016) and the WetCHARTs wetlands extended ensemble from Bloom et al. (2017). The inversion domain over the Southeast US accounts for 45 % of national methane emissions in the EPA inventory, and for 56 % of wetland emissions over the contiguous US in the mean WetCHARTs estimate.

Our inversion results suggest that the EPA emission inventory has no significant bias on the regional scale for the major anthropogenic source sectors (livestock, oil and gas, waste, coal), while the mean of the WetCHARTs wetland ensemble needs to be reduced by 27 % over the inversion domain. These results are supported by independent methane observations from the NOAA∕ESRL surface network and from the TCCON site in Lamont, Oklahoma. The specification of wetland areal extent is the dominant source of error in the WetCHARTs ensemble. Results also indicate that a low temperature dependence for the CH₄ : C heterotrophic respiration ratio best explains the spatial variability of the posterior emissions. Despite regional agreement with the EPA anthropogenic inventory, we still find significant local discrepancies with the EPA inventory for the oil and gas sector, suggesting that methane emission factors are more variable than assumed in the inventory.