Attribution of the accelerating increase in atmospheric methane during 2010–2018 by inverse analysis of GOSAT observations

We conduct a global inverse analysis of 2010–2018 GOSAT observations to better understand the factors controlling atmospheric methane and its accelerating increase over the 2010–2018 period. The inversion optimizes anthropogenic methane emissions and their 2010–2018 trends on a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid, monthly regional wetland emissions, and annual hemispheric concentrations of tropospheric OH (the main sink of methane). We use an analytical solution to the Bayesian optimization problem that provides closed-form estimates of error covariances and information content for the solution. We verify our inversion results with independent methane observations from the TCCON and NOAA networks. Our inversion successfully reproduces the interannual variability of the methane growth rate inferred from NOAA background sites. We find that prior estimates of fuel-related emissions reported by individual countries to the United Nations are too high for China (coal) and Russia (oil and gas) and too low for Venezuela (oil and gas) and the US (oil and gas). We show large 2010–2018 increases in anthropogenic methane emissions over South Asia, tropical Africa, and Brazil, coincident with rapidly growing livestock populations in these regions. We do not find a significant trend in anthropogenic emissions over regions with high rates of production or use of fossil methane, including the US, Russia, and Europe. Our results indicate that the peak methane growth rates in 2014–2015 are driven by low OH concentrations (2014) and high fire emissions (2015), while strong emissions from tropical (Amazon and tropical Africa) and boreal (Eurasia) wetlands combined with increasing anthropogenic emissions drive high growth rates in 2016–2018. Our best estimate is that OH did not contribute significantly to the 2010–2018 methane trend other than the 2014 spike, though error correlation with global anthropogenic emissions limits confidence in this result.

1 Introduction

Methane is the second most important anthropogenic greenhouse gas after CO₂, with an emission-based radiative forcing of 0.97 W m⁻² since preindustrial times (Myhre et al., 2013). Methane is emitted to the atmosphere from a range of anthropogenic activities including fuel exploitation, agriculture, waste and wastewater treatment, and biomass burning. The main natural source is from wetlands, with minor contributions from geological seeps, forest fires, and termites. Atmospheric methane has a lifetime of 11.2±1.3 years against tropospheric oxidation by the hydroxyl radical (OH) (Prather et al., 2012). Minor sinks include stratospheric loss, oxidation by Cl atoms, and absorption by soils (Kirschke et al., 2013).

Unlike the steady rise in atmospheric CO₂, the rise of methane has taken place in fits and starts. Observations from the NOAA network (Dlugokencky, 2020) (https://www.esrl.noaa.gov/gmd/ccgg/trends_ch4/, last access: 22 June 2020) show a period of stabilization in the early 2000s, followed by renewed growth after 2007 that has accelerated since 2014. Annual growth rates averaged 0.50 % a⁻¹ for 2014–2018 compared to 0.32 % a⁻¹ for 2007–2013. The growth of atmospheric methane concentrations, if continued at current rates in coming decades, may significantly negate the climate benefit of CO₂ emission reduction (Nisbet et al., 2019).

However, our understanding of the drivers behind the methane growth rate is still limited, preventing reliable projections for future changes. Explanations have differed for the renewed growth of atmospheric methane since 2007. A concurrent increase in atmospheric ethane has been interpreted as evidence of an increase in oil and gas emissions (Hausmann et al., 2016; Franco et al., 2016). However, the assumption that the ethane-to-methane emission ratio should be stable is questionable (Lan et al., 2019). Meanwhile, a concurrent shift towards isotopically lighter methane has been attributed to an increase in microbial sources either from livestock or wetlands (Schaefer et al., 2016; Nisbet et al., 2016). Worden et al. (2017) pointed out that the trend towards isotopically lighter methane could be explained by decreases in fire emissions that are isotopically heavy. Based on methyl chloroform observations, Turner et al. (2017) and Rigby et al. (2017) suggested that a decrease in the OH sink may be the cause of the methane regrowth.

To better interpret the methane budget and its recent trends, we present an inverse analysis of global 2010–2018 methane observations from the GOSAT instrument. GOSAT provides a long record (starting in 2009) of global high-quality observations of column methane mixing ratios (Kuze et al., 2016; Buchwitz et al., 2015). A number of inverse analyses previously used GOSAT observations to constrain methane emission estimates (Fraser et al., 2013; Monteil et al., 2013; Cressot et al., 2014; Alexe et al., 2015; Turner et al., 2015; Pandey et al., 2016, 2017a; Miller et al., 2019; F. Wang et al., 2019a; Lunt et al., 2019; Maasakkers et al., 2019; Janardanan et al., 2020; Tunnicliffe et al., 2020; Yin et al., 2020). Maasakkers et al. (2019) used 2010–2015 GOSAT observations to optimize gridded methane emissions, global OH concentrations, and their 2010–2015 trends. They concluded that increasing methane emissions were driven mainly by India, China, and tropical wetlands. Our analysis is based on that of Maasakkers et al. (2019) but extends it to 2018 in order to interpret the post-2014 acceleration. We implement for that purpose a number of major improvements to the Maasakkers et al. (2019) methodology including in particular (1) separate optimization of subcontinental wetland emissions to resolve their seasonal and interannual variability, (2) correction of stratospheric methane forward model biases based on ACE-FTS solar occultation satellite data (Waymark et al., 2014), (3) prior estimates of global fuel exploitation emissions using national reports submitted to the United Nations Framework Convention on Climate Change (UNFCCC) (Scarpelli et al., 2020), and (4) optimization of annual hemispheric OH concentrations.

2 Methods

2.1 GOSAT observations

The observation vector for the inversion (y) consists of column-averaged dry-air methane mole fractions during 2010–2018 observed by the TANSO-FTS instrument on board the Greenhouse Gases Observing Satellite (GOSAT) (Kuze et al., 2009). The satellite is in polar sun-synchronous low-Earth orbit and observes methane by nadir solar backscatter in the 1.65 µm shortwave infrared absorption band. Observations are made at around 13:00 local solar time. We use the University of Leicester version 9 CO₂ proxy retrieval (Parker et al., 2020a). The retrieval has been extensively validated against ground-based column observations from the Total Carbon Column Observing Network (Wunch et al., 2011). Validation has also been performed for the model XCO₂ used in the CO₂ proxy retrieval (Parker et al., 2015) and for a specific region (i.e., the Amazon) against aircraft profile observations (Webb et al., 2016). Overall, the retrieval has a single-observation precision of 13.7 ppb and a regional bias of 4 ppbv (Parker et al., 2020a), which is sufficient for a successful methane inversion (Buchwitz et al., 2015). The inversion ingests a total of 1.5 million successful GOSAT retrievals. Previous inversions of GOSAT data often excluded high-latitude observations because of seasonal bias, large retrieval errors at low solar elevations, and forward model errors for the stratosphere (Bergamaschi et al., 2013; Turner et al., 2015; Z. Wang et al., 2017; Maasakkers et al., 2019). The exclusion of high-latitude observations limited the capability of the inversions to resolve emissions at high latitudes such as from boreal wetlands and oil and gas activity in Russia (Maasakkers et al., 2019). Here we use an improved model bias correction scheme (Sect. 2.5) and include these high-latitude observations in the inversion.

2.2 State vector

The state vector (x) is the ensemble of variables that we seek to optimize in the inversion. In this work, the state vector includes (1) mean 2010–2018 methane emissions from non-wetland sources (all anthropogenic and natural emissions excluding wetlands) on a global $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid (1009 elements), (2) linear trends of non-wetland emissions on that same grid (1009 elements), (3) wetland emissions from 14 subcontinental regions for individual months (1512 elements) (Fig. 1), and (4) annual mean tropospheric OH concentrations in the Northern and Southern Hemisphere (18 elements). The reason to treat wetland and non-wetland emissions separately is that wetland emissions have large seasonal and interannual uncertainties (compared to anthropogenic emissions); coarsening the spatial resolution when optimizing wetland emissions allows us to estimate monthly values for individual years (Bloom et al., 2017). This is a significant improvement over the inverse analysis of Maasakkers et al. (2019), wherein interannual and seasonal errors in prior wetland emissions were not addressed by the inversion.

Another improvement in the state vector definition relative to Maasakkers et al. (2019) is to optimize annual mean OH concentrations in each hemisphere rather than just globally. Y. Zhang et al. (2018) previously found with an observing system simulation experiment that it should be possible to constrain annual mean hemispheric OH concentrations from satellite methane observations. Patra et al. (2014) suggested that global chemical transport models (CTMs) are often biased in their inter-hemispheric OH gradient relative to methyl chloroform observations, and such bias, if not corrected, would propagate to the solution for methane emissions.

Figure 1 Spatial distribution of mean 2010–2018 methane emissions used as prior estimates in the inversion of GOSAT data. The top panel shows wetland emissions, and the bottom panel shows non-wetland emissions. Blue boxes indicate the 14 subcontinental regions for which wetland emissions are optimized for individual months (Sect. 2.2): (1) Alaska + western Canada, (2) eastern Canada, (3) northern Europe, (4) Siberia, (5) temperate North America, (6) Latin America, (7) the Mediterranean, (8) East Asia, (9) the Amazon, (10) sub-Saharan Africa, (11) tropical South Asia, (12) Argentina, (13) southern Africa, and (14) Indonesia + Australia.

2.3 Prior estimates

Prior estimates for methane sources and sinks (x_a) are compiled from an ensemble of bottom-up studies. Figure 1 shows the spatial distribution of prior emission estimates. For gridded $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ anthropogenic emissions, we use as default the EDGAR v4.3.2 global emission inventory for 2012 (https://edgar.jrc.ec.europa.eu/, last access: 1 December 2017) (Janssens-Maenhout et al., 2017). We supersede it for the US with the gridded version of the Environmental Protection Agency greenhouse gas emission inventory for 2012 (Maasakkers et al., 2016). We further supersede it globally for fuel (oil, gas, and coal) exploitation with the inventory of Scarpelli et al. (2020) for 2012, which spatially disaggregates the national emissions reported to the United Nations Framework Convention on Climate Change (UNFCCC) (https://di.unfccc.int/, last access: 22 June 2020). All anthropogenic emissions are assumed to be aseasonal, except manure management for which we apply local temperature-dependent corrections (Maasakkers et al., 2016) and rice cultivation for which we apply gridded seasonal scaling factors from B. Zhang et al. (2016).

For the prior estimates of natural emissions, we take monthly wetland emissions during 2010–2018 from the WetCHARTS v1.0 extended ensemble mean (Bloom et al., 2017) for each subcontinental domain in Fig. 1. To test the impact of wetland spatial distribution within the subcontinental domains on inversion results, we performed a sensitivity inversion in which prior WetCHARTS emissions in Africa (regions 10 and 13 in Fig. 1) are increased by a factor of 3 in the Sudd wetland of South Sudan and decreased by a factor of 2.5 in the Congo Basin, following Lunt et al. (2019) and as shown in Fig. S1. Daily global emissions from open fires are taken from GFEDv4s (van der Werf et al., 2017), which accounts for high methane emissions from peatland fires (Liu et al., 2020). For geological sources, we scale the spatial distribution from Etiope et al. (2019) to a global total of 2 Tg a⁻¹ inferred from preindustrial-era ice core ¹⁴CH₄ data (Hmiel et al., 2020). Termite emissions are from Fung et al. (1991), totalling 12 Tg a⁻¹.

The prior estimates for 2010–2018 trends in non-wetland emissions are specified as zero on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid, except for interannual variability in fire emissions given by GFEDv4s. In this manner, all information on the posterior estimate of anthropogenic emission trends is from the GOSAT observations.

The prior estimates for the hemispheric tropospheric OH concentrations are based on a GEOS-Chem full chemistry simulation (Wecht et al., 2014). The monthly 3-D OH concentration fields from this full chemistry simulation are also used in the forward model. We optimize hemispheric OH concentrations as the methane loss frequency [s⁻¹] due to oxidation by tropospheric OH (kⁱ) in the Northern and Southern Hemisphere (i= north or south):

$$\begin{array}{}\text{(1)}& \begin{array}{l}{k}^i=\frac{{\int }_{\mathrm{troposphere},i}{k}^{\prime }\left(T\right)\left[\mathrm{OH}\right]n_{{\mathrm{CH}}_{\mathrm{4}}}\mathrm{d}v}{{\int }_{\mathrm{atmosphere}}{n}_{{\mathrm{CH}}_{\mathrm{4}}}\mathrm{d}v},\end{array}\end{array}$$

where $n_{{\mathrm{CH}}_{\mathrm{4}}}$ is methane number density [molec. cm⁻³], v is volume, and $k^{\prime }\left(T\right)=\mathrm{2.45}\times {\mathrm{10}}^{-\mathrm{12}}{e}^{-\mathrm{1775}/T}$ cm³ molec.⁻¹ s⁻¹ is the temperature-dependent oxidation rate constant (Burkholder et al., 2015). In this definition, the denominator of Eq. (1) integrates over the entire atmosphere, and the numerator integrates over the hemispheric troposphere. Hence, global methane loss frequency (or inverse lifetime; k) due to oxidation by tropospheric OH can be expressed as the sum of hemispheric values ($k=\mathrm{1}/\mathit{\tau }=k^{\mathrm{north}}+k^{\mathrm{south}}$, where τ is the global lifetime due to oxidation by tropospheric OH). Our prior estimates from Wecht et al. (2014) are 0.050 a⁻¹ for k^north and 0.043 a⁻¹ for k^south, which translates to a τ of 10.7 years and a north-to-south inter-hemispheric OH ratio of 1.16. In comparison, the methyl chloroform proxy infers τ of 11.2±1.3 years (Prather et al., 2012) and an inter-hemispheric ratio in the range 0.85–0.98 (Montzka et al., 2000; Prinn et al., 2001; Krol and Lelieveld, 2003; Bousquet et al., 2005; Patra et al., 2014), while the ACCMIP model ensemble yields a τ of 9.7±1.5 years and an inter-hemispheric ratio of 1.28±0.10 (Naik et al., 2013).

The Bayesian inversion requires error statistics for the prior estimates. The prior error covariance matrix (S_a) is constructed as follows. For mean non-wetland emissions, we assume 50 % error standard deviation for individual grid cells and 20 % for each source category when aggregated globally. For linear trends in non-wetland emissions, we specify an absolute error standard deviation of 5 % a⁻¹ for individual grid cells. For wetland emissions, we take the full error covariance structure (including spatial and temporal error covariance) derived from the WetCHARTs ensemble members for the 14 subcontinental regions (Bloom et al., 2017). For annual hemispheric OH concentrations, we assign 5 % independent errors for individual years on top of a 10 % error for the 2010–2018 mean.

2.4 Forward model

We use the GEOS-Chem CTM v11.02 as a forward model (F) for the inversion (Wecht et al., 2014; Turner et al., 2015; Maasakkers et al., 2019) that relates atmospheric methane observations (y) to the state vector to be optimized (x) as y=F(x). The simulation is conducted at $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ horizontal resolution with 47 vertical layers (∼30 layers in the troposphere) and is driven by 2009–2018 MERRA-2 meteorological fields (Gelaro et al., 2017) from the NASA Global Modeling and Assimilation Office (GMAO). The prior simulation is conducted from 2010 to 2018. The initial conditions are from Turner et al. (2015) with an additional 1-year spin-up starting from January 2009 to establish methane gradients driven by synoptic-scale transport (Turner et al., 2015). We further set the initial conditions on 1 January 2010 to be unbiased by removing the zonal mean biases relative to GOSAT observations. Thus, we attribute any model departures from observations over the 2010–2018 period in the inversion to errors in sources and sinks over that period.

We use archived 3-D monthly fields of OH concentrations from a GEOS-Chem full chemistry simulation (Wecht et al., 2014) to compute the removal of methane from oxidation by tropospheric OH. Other minor loss terms include stratospheric oxidation computed with archived monthly loss frequencies from the NASA Global Modeling Initiative model (Murray et al., 2012), tropospheric oxidation by Cl atoms computed with archived Cl concentration fields from X. Wang et al. (2019b), and monthly soil uptake fields from Murguia-Flores et al. (2018). The inversion does not optimize these minor sinks. The loss from oxidation by Cl is 5.5 Tg a⁻¹, accounting for ∼1 % of methane loss. It is lower than the previous estimate of 9 Tg a⁻¹ (Sherwen et al., 2016) used by Maasakkers et al. (2019) but is consistent with a recent analysis of methane and CO isotopic signatures (Gromov et al., 2018). Use of monthly soil uptake fields from the Murguia-Flores et al. (2018) climatology of 2000–2009 data is another update relative to Maasakkers et al. (2019) and yields a global soil sink of 34 Tg a⁻¹.

2.5 Forward model bias correction

The GEOS-Chem-simulated methane columns have a latitude-dependent background bias relative to the GOSAT data (Turner et al., 2015). This is thought to result from excessive meridional transport in the stratosphere, a common problem in global models (Patra et al., 2011). In particular, coarse-resolution global models have difficulty resolving polar vortex dynamics that control the distribution of stratospheric methane in the winter–spring hemisphere (Stanevich et al., 2020). The GEOS-Chem model evaluation with stratospheric sub-columns derived from ground-based TCCON measurements shows that the stratospheric bias varies seasonally (Saad et al., 2016). Previous GEOS-Chem-based inversions of GOSAT data (Turner et al., 2015; Maasakkers et al., 2019) developed correction schemes by fitting differences between the prior model simulation and background GOSAT observations as a second-order polynomial function of latitude. However, these correction schemes did not consider the seasonal variation of the stratospheric biases. Moreover, they could falsely attribute high-latitude model–GOSAT differences to stratospheric model bias rather than to errors in prior emissions. Therefore, previous studies excluded high-latitude observations from their analyses (Turner et al., 2015; Maasakkers et al., 2019).

Here we improve the stratospheric bias correction by using satellite observations from ACE-FTS v3.6 (Waymark et al., 2014; Koo et al., 2017). ACE-FTS is a solar occultation instrument launched in 2003 that measures vertical profiles of stratospheric methane (Bernath et al., 2005). We compute correction factors to GEOS-Chem stratospheric methane sub-columns as a function of season and equivalent latitude based on the ratios of stratospheric methane sub-columns between the ACE-FTS and GEOS-Chem prior simulation for 2010–2015 (Fig. 2). A global scaling factor (1.06) is applied to these correction factors to enforce mass conservation for methane in the stratosphere so that the correction does not introduce a spurious stratospheric source and sink in the model simulation. We use equivalent latitude, computed on the 450 K isentropic surface from MERRA-2 reanalysis fields, as one of the predictors for parameterization. The equivalent latitude is a coordinate based on potential vorticity (PV) that maps PV to latitude based on areas enclosed by PV isopleths (Butchart and Remsberg, 1986), and it is often used to represent the influence of high-altitude dynamics (e.g., polar vortex) on chemical tracers (Engel et al., 2006; Hegglin et al., 2006; Strahan et al., 2007). We use the same stratospheric bias correction for all years because the correction does not vary significantly for individual years (Fig. S2). Figure 2 shows that GEOS-Chem model biases are largely confined to high latitudes of the winter–spring hemisphere. By having our correction factors be dependent on equivalent latitude and season, we specifically account for the overly weak polar vortex dynamical barrier in the model as the cause of the stratospheric bias (Stanevich et al., 2020).

Figure 2 GEOS-Chem stratospheric bias correction based on ACE-FTS observations. The figure shows the ACE-FTS to GEOS-Chem ratio of stratospheric methane sub-columns as a function of equivalent latitude and season, averaged over the 2010–2015 period. Grey shading represents the fitting uncertainty.

2.6 Inversion procedure

We perform the inversion by minimizing the Bayesian cost function (Brasseur and Jacob, 2017):

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

Here, the Jacobian matrix $\mathbf{K}(=\frac{\partial \mathit{y}}{\partial \mathit{x}})$ is a linearized description of the forward model (F) that relates y (observations) to x (state vector). We explicitly compute the Jacobian matrix by perturbing each individual element of x independently in GEOS-Chem simulations and calculating the sensitivity of y to that perturbation. x_a is the prior estimate for x and S_a is the prior error covariance matrix (Sect. 2.3). S_O is the observation error covariance matrix including contributions from the instrument error and the forward model error. We take S_O to be diagonal and compute the variance terms from the statistics of the residual error (${\mathit{\epsilon }}_O=\mathit{y}-\mathit{F}\left({\mathit{x}}_a\right)-\stackrel{\mathrm{‾}}{\mathit{y}-\mathit{F}\left({\mathit{x}}_a\right)}$, where the overbar denotes annual averages in a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell) that represents the random component of model–observation differences (Heald et al., 2004). This method of constructing S_O has been previously applied to GOSAT observations by Turner et al. (2015) and Maasakkers et al. (2019). The observational error standard deviation derived in this manner averages 13 ppbv. γ is the regularization parameter taken to be 0.05 following Y. Zhang (2018) and Maasakkers et al. (2019) to account for missing error covariance structure in S_O.

Minimizing J(x) (Eq. 2) by solving $\mathrm{d}J/\mathrm{d}\mathit{x}=\mathrm{0}$ analytically (Rodgers, 2000; Brasseur and Jacob, 2017) yields a best posterior estimate of the state vector ($\widehat{\mathit{x}}$) and the associated posterior error covariance matrix ($\widehat{\mathbf{S}}$) characterizing the error statistics of $\widehat{\mathit{x}}$:

$$\begin{array}{}\text{(3)}& {\displaystyle }\widehat{\mathit{x}}={\mathit{x}}_a+{\left(\mathit{\gamma }{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}}\mathit{\gamma }{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{K}{\mathit{x}}_a\right),\text{(4)}& {\displaystyle }\widehat{\mathbf{S}}={\left(\mathit{\gamma }{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$

From there we derive the averaging kernel matrix $\mathbf{A}=\partial \widehat{\mathit{x}}/\partial \mathit{x}$ describing the sensitivity of the solution to the true state:

$$\begin{array}{}\text{(5)}& \mathbf{A}=\mathbf{I}-\widehat{\mathbf{S}}{\mathbf{S}}_a^{-\mathrm{1}}.\end{array}$$

The trace of the averaging kernel matrix is referred to as the degrees of freedom for signal (DOFS) (Rodgers, 2000) and represents the number of independent pieces of information on the state vector that are constrained by the inversion. We refer to the diagonal terms of A as averaging kernel sensitivities, which measure the ability of the observations to quantify the individual elements of the state vector (Sheng et al., 2018c; Maasakkers et al., 2019).

The posterior solution is often presented in reduced dimensionality. For example, spatially resolved emission and trend estimates on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid can be aggregated to countries, regions, or global and/or regional emissions from individual source sectors (oil and gas, livestock, etc.). Let W be a matrix that represents the linear transformation from the full state vector to a reduced state vector. The posterior estimation of the reduced state vector (${\widehat{\mathit{x}}}_{\mathrm{red}}$) is computed as

$$\begin{array}{}\text{(6)}& {\widehat{\mathit{x}}}_{\mathrm{red}}=\mathbf{W}\widehat{\mathit{x}},\end{array}$$

with posterior error covariance matrix

$$\begin{array}{}\text{(7)}& {\widehat{\mathbf{S}}}_{\mathrm{red}}=\mathbf{W}\widehat{\mathbf{S}}{\mathbf{W}}^{\mathrm{T}}\end{array}$$

and averaging kernel matrix

$$\begin{array}{}\text{(8)}& {\mathbf{A}}_{\mathrm{red}}={\mathbf{WAW}}^{*},\end{array}$$

where ${\mathbf{W}}^{*}={\mathbf{W}}^{\mathbf{T}}(\mathbf{W}{\mathbf{W}}^{\mathbf{T}}{)}^{-\mathrm{1}}$ is the pseudo-inverse of W. The regional and global budget terms and their error covariance structures obtained by using this approach are consistent with the full inversion. In the case of aggregation by sectors, we construct W on the basis of the relative contribution of the sector to the prior emissions in each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell. This does not assume that the prior distribution of sectoral emissions is correct, only that the relative allocation within a given $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell is correct.

3 Evaluation of the inversion results

We conduct a posterior simulation driven by the optimized (posterior) distributions of methane emissions, emission trends, and OH concentrations to evaluate the inversion. The posterior simulation results are compared with the training data (GOSAT) and independent evaluation data including TCCON total column measurements (https://tccondata.org/, last access: 22 June 2020) (Wunch et al., 2011) and NOAA surface measurements (https://www.esrl.noaa.gov/gmd/ccgg/flask.php, last access: 22 June 2020) (Dlugokencky et al., 2020). Figure 3 shows the GEOS-Chem comparison to the GOSAT data. As expected for a successful inversion, the posterior simulation achieves a better fit to GOSAT observations than the prior simulation both spatially and temporally, with root mean square errors reduced by 70 % (prior: 44 ppbv; posterior: 13 ppbv). The prior simulation shows a negative bias that grows with time and has a large seasonal structure presumably associated with errors in wetland emissions. The prior biases also have prominent spatial patterns, particularly in the extratropical Northern Hemisphere and the tropics. All these error features largely vanish in the posterior simulation through the optimized adjustment of the state vector (Fig. 3).

Figure 3 Difference in methane columns between GEOS-Chem simulations and GOSAT observations. Results are shown for GEOS-Chem using prior (a, c) and posterior (b, d) state vector estimates as well as for the spatial distribution averaged during 2010–2018 (a, b) and monthly time series of zonal means in different latitude bands (c, d). Note the different color scales in (a) and (b). Tick marks on the x axes in (c) and (d) represent January in each year. Figure S6 plots panel (d) in an expanded ordinate scale.

Figure 4 presents evaluations against independent 2010–2018 observations from TCCON and NOAA sites arranged by latitude. Values are shown as the root mean square error (RMSE) for individual sites. Figure 4 shows that the inversion substantially improves the agreement between simulations and observations for all TCCON sites and almost all NOAA surface sites. Average root mean square errors are reduced by 60 % for TCCON sites (prior: 38 ppbv; posterior: 15 ppbv) and by 42 % for NOAA surface sites (prior: 43 ppbv; posterior: 25 ppbv). The seasonal component of the errors (root mean square errors computed from monthly mean model–observation differences after annual mean biases are removed; not shown in the figure) is reduced on average by 42 % for TCCON sites (prior: 6.5 ppbv; posterior: 3.8 ppbv) and 30 % for surface sites (prior: 10 ppbv; posterior: 7 ppbv), primarily owing to optimized seasonal variations in wetland emissions. In addition, we do not find large interannual variability of posterior biases that could be associated with climate oscillations such as ENSO (Fig. S3).

Figure 4 Root mean square errors of prior and posterior GEOS-Chem simulations relative to TCCON observations of dry column methane mixing ratios (a) and NOAA observations of surface air mixing ratios (b, c). Observation sites are arranged by latitude. Data are for 2010–2018. Site names are shown along with their latitude and longitude (more information about these sites can be found at https://tccon-wiki.caltech.edu/ – last access: 22 June 2020 and https://www.esrl.noaa.gov/gmd/dv/site/index.php?program=ccgg – last access: 22 June 2020). A mountaintop TCCON site located at Zugspitze, Germany (zs; ∼3000 m a.s.l.), is excluded because the terrain effect on the total column is not resolved by the coarse-resolution model.

The posterior error covariance matrix $\widehat{\mathbf{S}}$ resulting from analytically solving the Bayesian problem allows us to analyze the error structure of posterior estimates. Figure 5 plots the posterior joint probability density functions (PDFs) for pairs of global budget terms and their trends (computed following Eqs. 6–7). A strong negative error correlation in the inversion results is found between global anthropogenic emissions and methane lifetime ($r=-\mathrm{0.8}$), reflecting the limited ability of the inversion to separate the two. In contrast, error correlations between wetland emissions and methane lifetime ($r=-\mathrm{0.4}$) as well as between wetland and anthropogenic emissions ($r=-\mathrm{0.2}$) are much smaller. We find moderate error correlations between the OH trend and either wetland or anthropogenic emission trends ($r=-\mathrm{0.6}$). Improved separation of global budget terms and their trends may be achieved by including additional information from surface observations (Lu et al., 2020) and from thermal infrared satellite observations (Y. Zhang et al., 2018).

Figure 5 Error correlations between global anthropogenic emissions, wetland emissions, and tropospheric OH concentrations (methane lifetime against oxidation by tropospheric OH; τ) in the inverse solution. Results are shown for both 2010–2018 mean values and 2010–2018 trends. The error correlations are presented as joint probability density functions for pairs of reduced global state vector elements. Confidence ellipses represent a probability of 0.1 (innermost) to 0.9 (outermost) at intervals of 0.1. The error correlation coefficients are inset.

Figure 6 further examines the error aliasing of estimates for anthropogenic and wetland emissions within or between regions. For this analysis, anthropogenic emissions optimized on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid are aggregated to 14 subcontinental regions used for estimating wetland emissions. We find only moderate negative error correlations ($r=-\mathrm{0.1}$ to −0.5) between estimates for anthropogenic and wetland emissions within the same region (diagonal of top left quadrant), suggesting that the inversion is able to separate the two. Cross-region error correlations are generally small for anthropogenic emissions (bottom left quadrant of Fig. 6) but have a complex structure for wetland emissions (top right quadrant of Fig. 6). For example, errors are positively correlated between sub-Saharan Africa and southern Africa wetlands (r=0.6) but are negatively correlated between eastern Canada and northern Europe wetlands ($r=-\mathrm{0.4}$).

Figure 6 Posterior error correlations between regional anthropogenic and wetland emissions. To examine error aliasing at a regional scale, anthropogenic emissions resolved on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid are aggregated to the 14 subcontinental regions in Fig. 1 used for optimizing wetland emissions. Numbers (1–14) indicate the region index as in Fig. 1. “A” and “W” stand for anthropogenic and wetland emissions, respectively. For example, 5A stands for anthropogenic emissions from temperate North America.

4 Results and discussion

4.1 Anthropogenic emissions

Figure 7 shows the correction factors from the inversion to 2010–2018 mean non-wetland emissions (posterior-to-prior ratios) along with the associated averaging kernel sensitivities (corresponding diagonal terms of the averaging kernel matrix). We achieve 179 independent pieces of information (DOFS) for constraining the emissions on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid. The spatial distribution of averaging kernel sensitivities largely follows the pattern of prior emissions (right panel of Fig. 5). The inversion provides strong constraints in major anthropogenic source regions such as East Asia, South Asia, and South America. The constraints are generally weaker over North America, Europe, and Africa, indicating that the inversion provides more diffuse spatial information in these regions.

We find that the prior emission inventory significantly overestimates anthropogenic emissions in eastern China (Fig. 7). The posterior estimate of Chinese anthropogenic emissions (47±1 Tg a⁻¹) is 30 % lower than the prior estimate (67 Tg a⁻¹) and is also lower than the latest national report from China to the UNFCCC of 55 Tg a⁻¹ for 2014 (Fig. 8). Based on the relative contribution of each sector in the prior inventory, we attribute ∼60 % of this downward correction to coal mining. The overestimation of anthropogenic emissions from China has been reported by previous global and regional GOSAT inversion studies using different EDGAR inventory versions as prior estimates (Monteil et al., 2013; Thompson et al., 2015; Turner et al., 2015; Maasakkers et al., 2019; Miller et al., 2019).

Figure 7 Corrections to prior estimates of 2010–2018 mean non-wetland methane emissions. (a) Posterior-to-prior emission ratios. Figure S7 shows the same corrections as posterior–prior emission differences. (b) Averaging kernel sensitivities (diagonal elements of the averaging kernel matrix). The averaging kernel sensitivities measure the ability of the observations to constrain the posterior solution (0: not at all, 1: fully). The sum of averaging kernel sensitivities defines the degrees of freedom for signal (DOFS), which is inset.

Another major downward correction is for the oil- and gas-producing regions in Russia. We estimate Russia's anthropogenic emissions to be 20±1 Tg a⁻¹ as opposed to the prior estimate of 34 Tg a⁻¹ (Fig. 8), and the difference is mainly attributable to the oil and gas sector (posterior: 15 Tg a⁻¹; prior: 27 Tg a⁻¹). This finding is consistent with Maasakkers et al. (2019), though they used a different oil and gas emission inventory. Russia has been revising downwards its national emission estimates submitted to the UNFCCC, and our posterior estimate of total anthropogenic emissions is closer to the country's latest submission to the UNFCCC for 2010–2018 (16 Tg a⁻¹; Fig. 8). The new submission revises oil and gas methane emissions downward by a factor of 3 relative to the previous submission used as a prior estimate in our inversion (Scarpelli et al., 2020).

We find large upward corrections to the prior inventory over East Africa (Mozambique, Zambia, Tanzania, Ethiopia, Uganda, Kenya, and Madagascar) and over South America (Brazil). A previous inversion suggested that corrections for these regions may be due to an underestimation of prior wetland emissions (Maasakkers et al., 2019). Our inversion, which optimizes wetland and anthropogenic emissions separately, attributes the underestimation to anthropogenic emissions (see also Sect. 4.3 for wetland results), though there is some error aliasing between them ($r=-\mathrm{0.5}$ for sub-Saharan Africa, −0.4 for southern Africa; Fig. 6). We find that the upward correction over eastern Africa generally remains robust in a sensitivity inversion whereby prior wetland emissions in a neighboring region (Sudd in South Sudan) are increased by a factor of 3 (Figs. S4 and 8). Based on prior sectoral information, these underestimations are most likely due to livestock emissions, whose bottom-up estimates have large uncertainties in these developing regions (Herrero et al., 2013).

Figure 8 National and regional estimates of 2010–2018 mean methane emissions from anthropogenic sources. Included are the top five individual countries in our posterior estimates, the European Union (including the United Kingdom), and East Africa (including Mozambique, Zambia, Tanzania, Ethiopia, Uganda, Kenya, and Madagascar). The UNFCCC record is from https://di.unfccc.int (last access: 10 July 2020). Non-Annex I countries do not report yearly emissions to the UNFCCC, and for those we use the latest UNFCCC submissions (Brazil, 2015; China, 2014; Ethiopia, 2013; India, Madagascar, Kenya, 2010; Uganda, Zambia, 2000; Mozambique, Tanzania, 1994). Inset are the averaging kernel sensitivities for these national and regional aggregated results, which measure the ability of the observations to constrain the posterior solution (0: not at all, 1: fully). The dot for East Africa represents the estimate inferred from a sensitivity inversion with the prior spatial distribution of African wetlands perturbed.

Another upward correction pattern in South America is located near Venezuela, a major oil-producing country in the region. The large correction in Venezuela likely reflects underestimation of fossil fuel emissions by the national estimates for 2010 reported to the UNFCCC. Upward corrections are also seen in central Asia (Iran, Turkmenistan), where the regional posterior estimates (10.1±0.9 Tg a⁻¹) are 49 % higher than the prior (6.8 Tg a⁻¹), with adjustments mainly attributed to the oil and gas sector. This region has previously been identified by satellite observations as a methane emission hot spot (Buchwitz et al., 2017; Varon et al., 2019; Schneising et al., 2020).

The inversion finds small upward corrections over the US (prior: 28 Tg a⁻¹; posterior: 33±1 Tg a⁻¹) (Fig. 8), resulting mainly from underestimation of emissions from the oil and gas sector in the western and south-central US (Fig. 7). This result is consistent with a high-resolution inversion over the US using the 2010–2015 GOSAT data, which spatially allocates the correction more specifically to oil and gas basins (Maasakkers et al., 2020). A number of previous studies have found that emissions from oil and gas production are underestimated in the national US inventory (e.g., Kort et al., 2014; Smith et al., 2017; Peischl et al., 2018; Alvarez et al., 2018; Y. Zhang et al., 2020; Gorchov Negron et al., 2020).

Small downward corrections with a diffuse pattern are found over Europe. The posterior estimate of anthropogenic emissions for the European Union (including the UK) is 16±0.7 Tg a⁻¹, slightly lower than the prior estimate (21 Tg a⁻¹) and the UNFCCC national reports (19 Tg a⁻¹ for 2014) (Fig. 8). Source sector attribution is difficult in this case because of spatial overlap between emission sectors. The inversion finds only small adjustments to prior emissions for India (prior: 32 Tg a⁻¹; posterior: 33±0.6 Tg a⁻¹) even though the information content is relatively large, confirming the prior inventory used in the inversion. Our estimate, however, is much higher than a previous inversion study for India (Ganesan et al., 2017) (22 Tg a⁻¹), the results of which are in agreement with India's UNFCCC report (20 Tg a⁻¹ for 2010) (Fig. 8). Our inversion attributes the discrepancy with the UNFCCC submission mainly to the livestock sector.

4.2 Anthropogenic emission trends

Figure 9 shows the spatial distribution of 2010–2018 trends for anthropogenic emissions inferred from GOSAT observations, along with the associated averaging kernel matrix sensitivities. The GOSAT observations provide 42 pieces of information to constrain the spatial distribution of anthropogenic emission trends, suggesting that, compared to mean emissions, the inversion is only able to constrain more diffuse spatial patterns for emission trends. These constraints are strongest in China and India, but there is also fairly strong information aggregated over other continental regions. The prior estimate assumed zero anthropogenic trends anywhere; therefore, the posterior trends are driven solely by the GOSAT information.

Figure 9 Anthropogenic methane emission trends for 2010–2018, as informed by GOSAT observations. (a) Relative emission trends on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid. Absolute emission trends are shown in Fig. S8. (b) Averaging kernel sensitivities that measure the ability of the observations to constrain the posterior solution (0: not at all, 1: fully). The degrees of freedom for signal (DOFS) are inset.

Significant positive trends of anthropogenic emissions are found in South Asia (0.58±0.16 Tg a⁻¹ a⁻¹ or 1.4±0.4 % a⁻¹; Pakistan and India), East Africa (0.22±0.10 Tg a⁻¹ a⁻¹ or 1.4±0.6 % a⁻¹; Ethiopia, Tanzania, Uganda, Kenya, and Sudan), West Africa (0.28±0.10 Tg a⁻¹ a⁻¹ or 4.4±1.4 % a⁻¹; Nigeria, Niger, Ghana, Côte d'Ivoire, Mali, Benin, Burkina Faso),and Brazil (0.19±0.15 Tg a⁻¹ a⁻¹ or 0.8±0.6 % a⁻¹). Based on prior sectoral fractions in each $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell, we attribute these positive trends mostly to the livestock sector (0.31 Tg a⁻¹ a⁻¹ in South Asia, 0.13 Tg a⁻¹ a⁻¹ in East Africa, 0.09 Tg a⁻¹ a⁻¹ in West Africa, and 0.17 Tg a⁻¹ a⁻¹ in Brazil). Indeed, according to the United Nations Food and Agriculture Office (UNFAO) statistics (http://www.fao.org/faostat, last access: 22 June 2020), all these regions have had rapid increases in livestock population. The fastest-growing cattle populations in the world reported by the UNFAO over the 2010–2016 period were in Pakistan (1.4×10⁶ heads a⁻¹), Ethiopia (1.2×10⁶ heads a⁻¹), Tanzania (1.1×10⁶ heads a⁻¹), and Brazil (0.9×10⁶ heads a⁻¹). Moreover, our inversion results for these regional trends in livestock emissions are generally consistent in magnitude with the trends from bottom-up livestock emission inventories (FAOSTAT, IPCC tier I method; EDGAR v4.3.2 and v5, hybrid tier I method; Chang et al., 2019, IPCC tier II method) (Fig. 10). Because our inversion assumes zero prior trends in anthropogenic emissions, the inferred trends are solely informed by satellite observations and are independent of the bottom-up trends in Fig. 10.

A positive trend in anthropogenic emissions (0.39±0.27 Tg a⁻¹ a⁻¹ or 0.8±0.6 % a⁻¹) is found over China driven by coal mining (northern China) and rice (southern China), but the magnitude of the trend is smaller than previous inverse analyses of satellite and surface observations for time horizons before 2015 (Bergamaschi et al., 2013; Thompson et al., 2015; Patra et al., 2016; Saunois et al., 2017; Miller et al., 2019; Maasakkers et al., 2019). We infer a much stronger trend (0.72±0.39 Tg a⁻¹ a⁻¹ or 1.6±0.9 % a⁻¹) for China if we restrict the GOSAT record to 2010–2016. Our results thus suggest that anthropogenic emission trends in China peaked midway within the 2010–2018 record. Indeed, coal production in China peaked in 2013 (Sheng et al., 2019).

Figure 10 Regional trends in anthropogenic methane emissions from livestock. Our GOSAT inversion results for 2010–2018 (with error standard deviations) are compared to estimates from different bottom-up inventories over the 2005–2017 period: Chang et al. (2019), FAOSTAT (2020), EDGAR v5 (Crippa et al., 2019), and EDGAR v4.3.2 (Janssens-Maenhout et al., 2017). Results are shown for South Asia (India and Pakistan), West Africa (Nigeria, Côte d'Ivoire, Mali, Niger, Burkina Faso, Cameroon, Ghana, and Benin), East Africa (Ethiopia, Kenya, Uganda, and Tanzania), and Brazil. By assuming zero prior trends in anthropogenic emissions, our inversion does not use trend information in any of these bottom-up inventories; the trends inferred by the inversion are solely from the GOSAT observations.

The inversion does not find significant 2010–2018 trends in anthropogenic emissions over the US ($-\mathrm{0.06}\pm \mathrm{0.21}$ Tg a⁻¹ a⁻¹, $-\mathrm{0.2}\pm \mathrm{0.6}$ % a⁻¹). This is generally consistent with the lack of a trend in US emissions over 2000–2014 in inversions collected by the Global Carbon Project (Bruhwiler et al., 2017) and over 2010–2015 in a North America regional inversion using GOSAT data (Maasakkers et al., 2020). It contradicts the 2 % a⁻¹–3 % a⁻¹ positive trend inferred from direct analysis of GOSAT enhancements (Turner et al., 2016; Sheng et al., 2018a) and the inference of large positive trends based on increasing concentrations of ethane and propane (Franco et al., 2016; Hausmann et al., 2016; Helmig et al., 2016). Bruhwiler et al. (2017) pointed out that the inference of methane emission trends from a simple analysis of GOSAT data could be subject to various biases including variability in background and seasonal sampling, which would be accounted for in an inversion. Increasing ethane-to-methane and propane-to-methane emission ratios over years may confound the inference of methane emission trends from ethane and propane records (Lan et al., 2019).

Small negative trends are found in central Asia (Uzbekistan, Kazakhstan, Turkmenistan, Afghanistan; $-\mathrm{0.17}\pm \mathrm{0.16}$ Tg a⁻¹ a⁻¹), Europe ($-\mathrm{0.19}\pm \mathrm{0.15}$ Tg a⁻¹ a⁻¹), and Australia ($-\mathrm{0.12}\pm \mathrm{0.07}$ Tg a⁻¹ a⁻¹). The decrease in central Asia is attributed mainly to oil and gas, and the decrease in Australia is attributed to coal mining and livestock. Trends over Europe and Russia are too diffuse to be separated by sectors. No significant trend is found in Russia ($-\mathrm{0.01}\pm \mathrm{0.25}$ Tg a⁻¹ a⁻¹).

4.3 Wetland emissions

From the inversion we infer wetland emissions for 14 subcontinental regions (Fig. 1) and for individual months from 2010 to 2018, thus allowing seasonal and interannual variability to be optimized. This achieves 167 independent pieces of information (DOFS) for wetland emissions. The results are presented as mean seasonal cycles (Fig. 11) and time series of annual means (Fig. 12). Recent studies have found that the mean WetCHARTs inventory used here as a prior estimate is too high in the Congo Basin and too low in the Sudd region (Lunt et al., 2019; Parker et al., 2020b; Pandey et al., 2021). Our inversion is unable to resolve this subregional spatial correction pattern because of coarse resolution in the wetland state vector (Fig. 1). We therefore perform a sensitivity inversion with modified prior estimates following Lunt et al (2019) (Fig. S1), which finds a 20 % (3 Tg a⁻¹) increase in estimates of the 2010–2018 average for sub-Saharan Africa and a 7 % (0.6 Tg a⁻¹) increase for southern Africa relative to the base inversion (Fig. S5). Interannual trends and seasonal cycles are almost unchanged between the two inversions (Fig. S5).

As shown in Figs. 11 and 12, our results find lower wetland emissions than the mean of the WetCHARTs ensemble (prior estimate) over the Amazon, temperate North America (the US), and eastern Canada. The previous inversion of GOSAT data by Maasakkers et al. (2019) also found overestimation of emissions by WetCHARTs in the coastal US and Canadian wetlands but did not have significant corrections over the Amazon. The overestimation of wetland emissions in the US and eastern Canada is also reported in analyses of aircraft measurements in the southeastern US (Sheng et al., 2018b) and surface observations in Canada (Baray et al., 2021), both of which used WetCHARTs v1.0 as prior information. The downward correction of North American emissions is consistent with recent WetCHARTs updates (v1.2.1) that substantially reduce methane emissions across regions categorized as partial wetland complexes (Lehner and Döll, 2004; Bloom et al., 2017).

The seasonal cycles inferred from the inversion are in general consistent with prior information (Fig. 11), although averaging kernel sensitivities indicate that we only have limited constraints on the seasonality, particularly for high-latitude regions in Northern Hemisphere winter. This was generally expected given the limited GOSAT observational coverage at high latitudes during winter months. The inversion infers a sharp and late (May–June) onset of methane emissions across boreal wetlands, in contrast to an early and gradual increase starting from March predicted by the prior inventory. This feature in posterior estimates appears to be consistent with aircraft and surface observations over Canada's Hudson Bay Lowlands (Pickett-Heaps et al., 2011) and eddy flux measurements over Alaskan Arctic tundra (Zona et al., 2016), while the gradual onset in the prior inventory agrees with aircraft observations over Alaska by Miller et al. (2016). The negative fluxes right before the onset may reflect strong soil sinks during spring thaw over these high-latitude regions (Jørgensen et al., 2015). The inversion also indicates stronger seasonal cycles than the prior inventory in sub-Saharan Africa and tropical South Asia, which suggests that WetCHARTs may underestimate the sensitivity of wetland emissions to precipitation but overestimate their sensitivity to temperature.

Our posterior estimates of 2010–2018 trends in wetland emissions vary greatly by region and are generally larger than the prior estimates from WetCHARTs (Fig. 12). Large positive trends are found in the tropics (Amazon: +1.0 Tg a⁻¹ a⁻¹; sub-Saharan Africa: +0.6 Tg a⁻¹ a⁻¹; southern Africa: +0.4 Tg a⁻¹ a⁻¹) and high latitudes (Siberia: +0.4 Tg a⁻¹ a⁻¹). Increasing Amazonian wetland emissions may have been driven by intensification of flooding (Barichivich et al., 2018) or increasing temperature (Tunnicliffe et al., 2020) in the region over the past decades. Our result of increasing tropical Africa wetland emissions is consistent with a previous regional analysis of GOSAT data, which found a positive trend of 1.5–2.1 Tg a⁻¹ a⁻¹ in the region for 2010–2016 attributed mainly to wetlands, particularly the Sudd in South Sudan (Lunt et al., 2019). Compared to steady and linear increases in the tropics, boreal Siberia and northern Europe show abrupt increases in 2016–2018 for reasons that are unclear (Fig. 12). Decreasing but weaker trends are found in tropical Southeast Asia (−0.2 Tg a⁻¹ a⁻¹) and Australia (−0.1 Tg a⁻¹ a⁻¹).

Figure 11 Seasonal variation in wetland emissions for 14 subcontinental regions (Fig. 1). Values are means for 2010–2018. The prior estimate is the mean of the WetCHARTs inventory ensemble (Bloom et al., 2017). The posterior estimate is from our inversion of GOSAT data. Vertical bars are error standard deviations. The reduction of error in the posterior estimate measures the information content provided by the GOSAT data. Scales are different between panels.

Figure 12 Wetland emission trends for 2010–2018. The figure shows annual mean emissions for the prior estimate (mean of WetCHARTs inventory ensemble) and the posterior estimate after inversion of GOSAT data. Values are for the 14 subcontinental regions in Fig. 1 (panels with a white background) and also aggregated for the extratropics and tropics (panels with a grey background). The trends are from ordinary linear regression. Inset are prior and posterior 2010–2018 average annual emissions (Tg a⁻¹) with 2010–2018 trends (Tg a⁻¹ a⁻¹) in parentheses. Significant trends at the 95 % confidence level are denoted with *. Note the differences in scales between panels.

4.3.1 OH concentration

Our inversion infers a global methane lifetime against oxidation by tropospheric OH of 12.4±0.3 a, which is 15 % longer than the prior estimate (10.7±1.1 a) and is at the higher end of the inference from the methyl chloroform proxy (11.2±1.3 years) (Prather et al., 2012). We find that the downward correction for OH concentrations is mainly in the Northern Hemisphere. The north-to-south inter-hemispheric OH ratio is 1.02±0.05 in the posterior estimate compared to 1.16 in the prior estimate and 1.28±0.10 in the ACCMIP model ensemble (Naik et al., 2013). It is more consistent with the observation-based estimate of 0.97±0.12 (Patra et al., 2014). No significant 2010–2018 trend is seen in the methane lifetime (Fig. 13). The OH concentration in 2014 is 5 % lower than average, which may relate to particularly large peatland fires over Indonesia in 2014 that would be very large sources of carbon monoxide (CO) as a sink for OH (Pandey et al., 2017b).

Figure 13 Methane loss frequency and lifetime against oxidation by tropospheric OH for 2010–2018. Values are annual means with error standard deviations. The loss frequency (k) is as calculated by Eq. (1) and the lifetime (τ) is the inverse. The prior estimate from Wecht et al. (2014) assumes no 2010–2018 trend in OH concentrations; the slight variability seen in the figure is due to temperature.

4.4 Attribution of the 2010–2018 methane trend

Figure 14 shows the 2010–2018 annual methane growth rates inferred from NOAA background surface measurements (Dlugokencky, 2020) and from our inversion of GOSAT data. There is general consistency between the two, with both showing the highest growth rates in 2014–2015 and a general acceleration after 2014. Our inversion achieves a much better match to the interannual variability of the NOAA record than the previous work of Maasakkers et al. (2019), in large part because of our optimization of the interannual variability of wetland emissions.

The bottom panel of Fig. 14 attributes the interannual variability in the methane growth rate to individual processes as determined by the inversion. The growth rate $G_j=[\mathrm{d}m/\mathrm{d}t{]}_j$ in year j (where m is the global methane mass) is determined by the balance between sources and sinks:

$$\begin{array}{}\text{(9)}& G_j=E_j-k_j{m}_j-L_j.\end{array}$$

Here, E_j denotes the global emissions in year j, for which the inversion provides further sectoral breakdown, k_j is the loss frequency against oxidation by tropospheric OH (Eq. 1), m_j is the total methane mass, and L_j represents the minor sinks not optimized by the inversion. Let $\mathrm{\Delta }{E}_j=E_j-E_o$, $\mathrm{\Delta }{k}_j=k_j-k_o$, and $\mathrm{\Delta }{m}_j=m_j-m_o$ represent the changes relative to 2010 conditions (E_o, k_o, m_o) taken as a reference. We can then write

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& G_j=\left(E_o+\mathrm{\Delta }{E}_j\right)-\left(k_o+\mathrm{\Delta }{k}_j\right)\left(m_o+\mathrm{\Delta }{m}_j\right)-(L_o+\mathrm{\Delta }{L}_j)\\ & \approx \left(E_o-k_o{m}_o-L_o-k_o\mathrm{\Delta }{m}_j\right)+\mathrm{\Delta }{E}_j-m_o\mathrm{\Delta }{k}_j,\end{array}\end{array}$$

where we have neglected the minor terms Δk_jΔm_j and ΔL_j. Here, the growth rate G_j in year j is decomposed into three terms: (1) a relaxation to steady state based on 2010 conditions ($E_o-k_o{m}_o-L_o-k_o\mathrm{\Delta }{m}_j$), (2) a perturbation to emissions (ΔE_j) that can be further disaggregated by sectors, and (3) a perturbation to OH concentrations (m_oΔk_j).

We see from the bottom panel of Fig. 14 that the legacy of the 2010 imbalance sustains a positive growth rate throughout the 2010–2018 period, but this influence diminishes towards the end of the record. The 2010–2018 trend in anthropogenic emissions is linear by design of the inversion and sustains a steady growth rate over the 2010–2018 period as the legacy of the 2010 imbalance declines. Figure 15 shows the sectoral breakdown of the anthropogenic trend. The increase in global anthropogenic emissions totalling 1.9±0.8 Tg a⁻¹ a⁻¹ is driven by livestock (0.70±0.36 Tg a⁻¹ a⁻¹; South Asia, tropical Africa, Brazil), rice (0.44±0.18 Tg a⁻¹ a⁻¹; East Asia), and wastewater treatment (0.33±0.13 Tg a⁻¹ a⁻¹; Asia). We find an insignificant positive global trend in emissions from fuel exploitation (oil, gas, and coal) (0.18±0.4 Tg a⁻¹ a⁻¹) (Fig. 15).

The bottom panel of Fig. 14 also shows that the spike in the methane growth rate in 2014–2015 is attributed to an anomalously low OH concentration in 2014 (5 % lower than 2010–2018 average; Fig. 13), partly offset by low wetland emissions and anomalously high fire emissions in 2015, mostly from Indonesia peatlands (Worden et al., 2017). Smoldering peatland fires are particularly large sources of methane (Liu et al., 2020). The large fire emissions are informed by the GFED inventory because the interannual variability of fire emissions is not optimized in our inversion. Despite their small magnitude relative to wetland and anthropogenic emissions globally, anomalous fire emissions can be an important contributor to methane interannual variability (Worden et al., 2017; Pandey et al., 2017b) both directly by releasing methane and indirectly by decreasing OH concentrations through CO emissions.

In addition to the 2014–2015 extremum, the NOAA surface observations show an acceleration of methane growth during the latter part of the 2010–2018 record (Nisbet et al., 2019), and this is reproduced in our inversion (Fig. 14). This accelerating growth is driven by strong wetland emissions, particularly in 2016–2018, on top of increasing anthropogenic emissions (Fig. 14). Our inversion infers a relatively steady 2010–2018 increase from tropical wetlands (in particular the Amazon and tropical Africa) and a 2016–2018 surge from extratropical wetlands (in particular boreal Eurasia) (Fig. 12). More work is needed to understand the drivers of changes in wetland emissions.

Figure 14 The 2010–2018 annual growth rates in global atmospheric methane. (a) Comparison of annual growth rates inferred from our inversion of GOSAT data and from the NOAA surface network (Dlugokencky, 2020). Average methane growth rates for the period are inset. (b) Attribution of annual growth rates in the GOSAT inversion to perturbations of emissions (anthropogenic, wetlands, fires) and OH concentrations relative to 2010 conditions. The purple bar shows the relaxation of the 2010 budget imbalance to steady state. See the text for details explaining the breakdown.

We estimate from the inversion global mean methane emissions for 2010–2018 of 512±4 Tg a⁻¹ (wetlands: 145 Tg a⁻¹; anthropogenic: 336 Tg a⁻¹) and a sink of 489±4 Tg a⁻¹. This posterior global emission is lower than the prior estimate (538 Tg a⁻¹) and the 538–593 Tg a⁻¹ range reported by the Global Carbon Project for 2008–2017 (Saunois et al., 2020). Compared to prior emissions, we estimate lower emissions for wetlands and fossil fuel and higher emissions for livestock and rice (Figs. 12 and 15). Meanwhile, we estimate a methane lifetime against tropospheric OH oxidation of 12.4±0.3 years, which is at the high end of 11.2±1.3 years based on the methyl chloroform proxy (Prather et al., 2012), though strong error correlations between wetland emissions and methane lifetime suggest that our inversion has limited ability to constrain both independently (Fig. 5).

Figure 15 The 2010–2018 global methane anthropogenic emissions and emission trends partitioned by individual sectors. Posterior estimates are from our inversion of GOSAT data. Prior estimates for anthropogenic emission trends are zero. Error bars in (b) show posterior error standard deviations for emission trends. Posterior error standard deviations for mean emissions are small and are thus not shown in (a).

5 Conclusions

We quantified the regional and sectoral contributions to global atmospheric methane and its 2010–2018 trend through the inversion of GOSAT observations. The inversion jointly optimizes (1) 2010–2018 anthropogenic emissions and their linear trends on a $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid, (2) wetland emissions in 14 subcontinental regions for individual months, and (3) annual mean hemispheric OH concentrations for individual years. An analytical solution to the optimization problem provides closed-form estimates of posterior error covariances and information content, allowing us in particular to diagnose error correlations in our solution. The separate optimization of wetland and anthropogenic emissions allows us to resolve interannual and seasonal variations in posterior wetland emissions. Our inversion introduces additional innovations, including the correction of stratospheric model biases using ACE-FTS satellite data, and a new bottom-up inventory for emissions from fossil fuel exploitation based on national reports to the UNFCCC (Scarpelli et al., 2020).

Our optimization of 2010–2018 mean anthropogenic emissions on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid provides strong information in source regions as measured by averaging kernel sensitivities. We find that estimates of anthropogenic emissions reported by individual countries to the UNFCCC are too high for China (coal emissions) and Russia (oil and gas emissions) and too low for Venezuela (oil and gas) and the US (oil and gas). We also find that tropical livestock emissions are larger than previous estimates, particularly in South Asia, Africa, and Brazil. Our posterior estimate of anthropogenic emissions in India (33 Tg a⁻¹) is much higher than its most recent (2010) report to the UNFCCC (20 Tg a⁻¹), mostly because of livestock emissions.

The 2010–2018 trends in methane emissions on the $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid are successfully quantified in source regions. We find that large growth in anthropogenic emissions occurs in tropical regions including South Asia, tropical Africa, and Brazil that can be attributed to the livestock sector. This finding is consistent with trends in livestock populations. There has been little discussion in the literature about increasing agricultural methane emissions in these developing countries (Jackson et al., 2020). Our results also show a 2010–2018 increase in Chinese emissions, but the inferred rate of the increase is smaller than previously reported in inversions focused on earlier periods, likely caused by leveling of coal emissions in China. The 2010–2018 emission trend in the US is insignificant on the national scale.

We find that global wetland emissions are lower than the mean WetCHARTs emissions used as a prior estimate, mostly because of the Amazon. Wetland emissions over North America are also lower, consistent with previous studies. In both cases, posterior estimates are all well within the full WetCHARTs uncertainty range (Bloom et al., 2017). The seasonality of wetland emissions inferred by the inversion is in general consistent with WetCHARTs. An exception is in boreal wetlands where we find negative fluxes in April–May, possibly reflecting methane uptake as the soil thaws. The inversion infers increasing wetland emissions over the 2010–2018 period, superimposed on large interannual variability, in both the tropics (the Amazon, tropical Africa) and extratropics (Siberia).

Our optimization of annual hemispheric OH concentrations yields a global methane lifetime of 12.4±0.3 years against oxidation by tropospheric OH, with an inter-hemispheric OH ratio of 1.02. Our best estimate is that the global OH concentration has no significant trend over 2010–2018 except for a 5 % dip in 2014.

Taking all these methane budget terms together, our inversion of GOSAT data estimates global mean methane emissions for 2010–2018 of 512 Tg a⁻¹, with 336 Tg a⁻¹ from anthropogenic sources, 145 Tg a⁻¹ from wetland sources, and 31 Tg a⁻¹ from other natural sources. Our inferred growth rate of methane over that period matches that observed at NOAA background sites, including peak growth rates in 2014–2015 and an overall acceleration over the 2010–2018 period. We attribute the 2014–2015 peaks in methane growth rates to low OH concentrations (2014) and high fire emissions (2015), and we attribute the overall acceleration to a sustained increase in anthropogenic emissions over the period and strong wetland emissions in the latter part of the period. Most of the increase in anthropogenic emissions is attributed to livestock (in tropics), with contributions from increases in rice and wastewater emissions (Asia). Our best estimate indicates a positive trend from fuel exploitation, but this trend is statistically insignificant given the uncertainty of the inversion. Our finding is in general consistent with a previous 2010–2015 inversion of GOSAT data (Maasakkers et al., 2019), although here we use a longer record and capture the interannual variability better. Our results also agree with isotopic data, indicating that the rise in methane is driven by biogenic sources (Schaefer et al., 2016; Nisbet et al., 2016). The increase in tropical livestock emissions is quantitatively consistent with bottom-up estimates. More work is needed to understand interannual variations in wetland emissions.