Accelerating methane growth rate from 2010 to 2017: leading contributions from the tropics and East Asia

After stagnating in the early 2000s, the atmospheric methane growth rate has been positive since 2007 with a significant acceleration starting in 2014. While the causes for previous growth rate variations are still not well determined, this recent increase can be studied with dense surface and satellite observations. Here, we use an ensemble of six multi-species atmospheric inversions that have the capacity to assimilate observations of the main species in the methane oxidation chain – namely, methane, formaldehyde, and carbon monoxide – to simultaneously optimize both the methane sources and sinks at each model grid. We show that the surge of the atmospheric growth rate between 2010–2013 and 2014–2017 is most likely explained by an increase of global CH₄ emissions by 17.5±1.5 Tg yr⁻¹ (mean ± 1σ), while variations in the hydroxyl radicals (OH) remained small. The inferred emission increase is consistently supported by both surface and satellite observations, with leading contributions from the tropical wetlands (∼ 35 %) and anthropogenic emissions in China (∼ 20 %). Such a high consecutive atmospheric growth rate has not been observed since the 1980s and corresponds to unprecedented global total CH₄ emissions.

1 Introduction

Methane (CH₄) is an important greenhouse gas highly relevant to climate mitigation, given its stronger warming potential and shorter lifetime than carbon dioxide (CO₂) (IPCC, 2013). Atmospheric levels of methane, usually measured as a dry air mole fraction, have nearly tripled since the industrial revolution according to ice core records (Etheridge et al., 1998; Rubino et al., 2019). This increase is mostly due to increases in anthropogenic emissions from agriculture (ruminant livestock and rice farming), fossil fuel use, and waste processing (Kirschke et al., 2013; Saunois et al., 2016; Schaefer, 2019). A significant portion of methane is also emitted from natural sources, including wetlands, inland freshwaters, geological sources, and biomass burning (although many of the wildfires may have anthropogenic origins) (Saunois et al., 2016). Methane has a lifetime of around 10 years in the atmosphere (Naik et al., 2013), with a dominant sink from oxidation by hydroxyl radicals (OH) in the troposphere (∼ 90 % of the total sink) (Saunois et al., 2019). In addition, its reactions with atomic chlorine (Cl), soil deposition, and stratospheric loss through reaction with a range of reactants (including O(¹D), Cl and OH) account for a minor portion of the total methane sink (Saunois et al., 2019).

Since the beginning of the direct measurement period in the early 1980s, the methane growth rate had been gradually declining until it reached a stagnation between the late 1990s and 2006, often referred to as the ”stabilization” period (Dlugokencky et al., 1998, 2003). However, methane has been increasing again since 2007 (Dlugokencky et al., 2009; Nisbet et al., 2014). A sharp increase of the growth rate was observed in 2014 from surface background stations (12.6 ± 0.5 ppb yr⁻¹; mean ± 1σ) (Nisbet et al., 2016; Fletcher and Schaefer, 2019; Nisbet et al., 2019), more than twice the average growth rate of 5.7 ± 1.1 ppb yr⁻¹ during the post stagnation period between 2007 and 2013. Since then, the CH₄ growth rate has remained high (8.6 ± 1.6 ppb yr⁻¹ for 2014–2017). Understanding the methane source and sink changes underlying the variations in methane growth rate can help us identify how methane sources respond to human activity, climate, or environmental changes, which are critical to climate mitigation efforts.

The attribution of the plateau and regrowth in atmospheric methane during the 2000s reached conflicting conclusions about the role of fossil fuel emissions (Hausmann et al., 2016; Simpson et al., 2012; Worden et al., 2017), agriculture or wetland emissions (Nisbet et al., 2016; Saunois et al., 2017; Schaefer et al., 2016), OH concentration (Rigby et al., 2017; Turner et al., 2017), and biospheric sinks (Thompson et al., 2018). The range of competing explanations exemplifies the complexity and uncertainty of interpolating limited observations of atmospheric methane and its ¹³C $/$ ¹²C isotopic ratio (expressed as δ¹³CH₄) to changes in different sectors of methane sources as well as its sinks (Turner et al., 2019; Schaefer, 2019). The situation now is more encouraging than the previous decade as we have continuous global satellite retrievals of the total column CH₄ dry air mole fraction (denoted as X${}_{{\mathrm{CH}}_{\mathrm{4}}}$) from the Greenhouse Gases Observing Satellite (GOSAT) with better precision and accuracy than previous instruments (Kuze et al., 2009; Parker et al., 2015; Jacob et al., 2016; Buchwitz et al., 2017; Houweling et al., 2017). The combined information from satellite and surface observations provides us a unique opportunity to understand the recent changes in methane growth rate with better spatial coverage.

Atmospheric methane measurements can be linked quantitatively to regional sources and sinks by inverse modeling, where changes in the atmospheric transport are guided by meteorological reanalysis and fluxes are adjusted to match the temporal and spatial variations of the observations given their uncertainties in a Bayesian formalism (Chevallier et al., 2005). A number of inverse studies have explored the surface and GOSAT observations to improve methane emission estimates (Monteil et al., 2013; Cressot et al., 2014; Alexe et al., 2015; Miller et al., 2019; Ganesan et al., 2017; Maasakkers et al., 2019), but the recent acceleration of methane growth since 2014 has not been widely investigated (Nisbet et al., 2019; McNorton et al., 2018; Turner et al., 2019; Zhang et al., 2021). GOSAT satellite X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ retrievals agree with the surface methane observations on the acceleration of the increase in atmospheric methane burden over the period from mid-2009 to the end of 2017 (Fig. 1). However, the satellite column data show a smoother temporal variation in the global average growth rate. GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ growth rates over different regions show diverse temporal patterns with a higher variability than the global average (Fig. S1 in the Supplement), suggesting that satellite data sampling directly over the source regions could provide valuable information to track regional changes in CH₄ fluxes. Furthermore, species in the oxidation chain of methane – namely, methane–formaldehyde–carbon monoxide (CH₄–HCHO–CO) – with their reactions to OH as the common sink path, could provide additional constraints on the OH sink of methane. A recent study showed that HCHO levels can provide information about remote tropospheric OH concentrations (Wolfe et al., 2019), and the feedback of CO variations on OH is directly linked to the sink of CH₄ (Gaubert et al., 2017; Nguyen et al., 2020). Hence, satellite retrievals of X_HCHO from the Ozone Monitoring Instrument (OMI) (González Abad et al., 2016) and X_CO from the Measurements of Pollution in the Troposphere (MOPITT) instrument (Deeter et al., 2017) covering the study period could, in theory, provide additional constraints on regional variations of methane sinks.

Figure 1 Atmospheric methane mixing ratio changes. (a) Monthly time series of the global mean methane mixing ratio from mid-2009 to the end of 2017. The green curve represents methane mixing ratios in the marine boundary layer observed by the NOAA surface network (https://gml.noaa.gov/ccgg/mbl/, last access: 1 May 2020). The red curve represents the total column mixing ratio, X${}_{{\mathrm{CH}}_{\mathrm{4}}}$, seen by GOSAT satellite and averaged from all soundings over the land. The smooth curve fit shows a quadratic fit of the trend that accelerates in the latter part of the study period. (b) Smooth methane growth rate derived from the time series shown in (a) following the methods of Thoning et al. (1989).

Hence, we developed a multi-tracer variational inverse system, PYVAR-LMDZ, with the capacity to assimilate observations of the CH₄–HCHO–CO oxidation chain to better constrain the sources and sinks of these species at individual model grid cells (Chevallier et al., 2005; Pison et al., 2009; Fortems-Cheiney et al., 2012; Yin et al., 2015; Zheng et al., 2019). Given the observed changes in temporal and spatial variations of all three species, we optimize simultaneously (i) methane emissions, (ii) CO emissions, (iii) HCHO sources (surface emissions + chemical productions from VOC oxidation), and (iv) OH concentrations. These terms are optimized at a weekly temporal resolution and a 1.9^∘ by 3.75^∘ spatial resolution. In addition, we optimize the initial concentrations of all four species on the individual model grid cell. Here, we performed an ensemble of six inversions using different combinations of observational constraints (surface vs. satellite, single vs. multiple species) and alternative prior estimates of 3-D OH distributions. With the ensemble results we aim to (1) identify key regions that contribute to the methane growth rate acceleration from 2010 to 2017 and (2) evaluate the consistency of results inferred from surface and satellite observations. Inversion methods and observational datasets are documented in Sect. 2. We report estimates of the global methane budget change from 2010 to 2017 in Sect. 3 and discuss regional attributions and sources of uncertainties in Sect. 4. Section 5 summarizes this work and provides some perspectives for future studies.

2 Data and methods

2.1 Atmospheric observations

We assimilate surface and satellite methane observations in parallel to test the consistency of information brought by these two types of measurements. We also include versions assimilating HCHO and CO along with CH₄ to test the impacts of adding chemically related species. In total, there are three groups of observational constraints:

• S1

  Methane and CO measurements from surface stations;

• S2

  GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$;

• S3

  GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$, OMI X_HCHO, and MOPITT X_CO.

The assimilation is done from April 2009 to February 2018, and we analyze the results of the 8 full years of 2010–2017 with the starting and ending period being spin-up and spin-down phases to avoid the edge effect.

2.1.1 Surface observations

We include surface methane observations from a total of 103 stations (Fig. S2; Table S3), with leading contributions from the following networks: the US National Oceanic and Atmospheric Administration (NOAA, 58 stations), Australia's Commonwealth Scientific and Industrial Research Organisation (CSIRO, nine stations), Environment and Climate Change Canada (ECCC, eight stations), and the Advanced Global Atmospheric Gases Experiment (AGAGE, five stations) (Prinn et al., 2018). Measurements from different networks are calibrated to the WMO scale. Daily afternoon averages between 12:00 and 18:00 local time are used for the assimilation of the continuous in situ measurements to minimize uncertainties associated with boundary layer height modeling. CO observations from those stations are also assimilated in S1.

2.1.2 Satellite observations

The TANSO-FTS instrument onboard the Greenhouse Gases Observing Satellite (GOSAT) was launched by The Japan Aerospace Exploration Agency (JAXA) into a polar sun-synchronous orbit in early 2009. It observes column-averaged dry-air carbon dioxide and methane mixing ratios by solar backscatter in the shortwave infrared (SWIR) with near-unit sensitivity across the air column down to the surface (Butz et al., 2011; Kuze et al., 2016). Observations are made at around 13:00 local time with a circular pixel around 10 km in diameter. The distances between pixels both along and cross track are ∼ 250 km in the default observation mode, and the revisit time for the same observation location is 3 d. Denser observations over particular areas of interest are made in target mode. Here, we use GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ proxy retrievals (OCPR) version 7.2 from the University of Leicester, which has been well documented and evaluated against various observations. The retrieval has a single-observation precision of 14 ppb (∼ 0.7 %) and a regional bias of ∼ 4 ppb compared to TCCON stations (Parker et al., 2015, 2020). This product is also consistent with other GOSAT methane retrievals (Buchwitz et al., 2017). However, we note that there is limited spatial coverage of TCCON stations to fully evaluate GOSAT observations in the high latitudes and the tropics. We only assimilate GOSAT retrievals over land to minimize potential retrieval biases between nadir and glint viewing modes. The same GOSAT data are assimilated in both S2 and S3.

For the multi-tracer inversion, S3, we also include OMI X_HCHO retrievals version 3 from the Smithsonian Astrophysical Observatory (SAO) (González Abad et al., 2016) and MOPITT X_CO retrievals version 7 from NCAR (Deeter et al., 2017). All satellite retrievals are processed following the recommend quality flags and the application of corresponding prior profiles and retrieval averaging kernels when provided. We exclude data poleward of 60^∘. Individual retrievals that are located in the same model grid within 3 h intervals are averaged for further assimilation. The observation uncertainty contains the retrieval errors as reported by the data product plus model errors whose standard deviations are empirically set as 1 % for CH₄, 30 % for CO, and 30 % for HCHO based on previous experiments (Fortems-Cheiney et al., 2012; Cressot et al., 2014; Yin et al., 2015).

Satellite retrievals of the three species (CH₄, HCHO, and CO) we use here are generally sensitive to the entire vertical column with some differences toward the lower troposphere. GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ retrievals using shortwave infrared (SWIR) radiances have approximately uniform sensitivity to methane at all pressure levels (Parker et al., 2015). OMI HCHO retrievals using ultraviolet (UV) radiance are sensitive to the entire column with some decline in the lowest atmospheric layers (González Abad et al., 2016). For MOPITT, we use the multispectral total column CO retrieval products that combine near-infrared (NIR) and thermal infrared (TIR) radiances and hence have an enhanced sensitivity to the lower troposphere (Deeter et al., 2014). Such subtle differences in the vertical sensitivities of the three retrievals, as well as their different vertical profiles and lifetimes, may influence the ways the observations of the three species provide information about OH, which is another source of uncertainty in addition to the model and observation errors.

2.1.3 Ground-based total column measurements

Ground-based X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ retrievals from the Total Carbon Column Observing Network (TCCON) from 27 stations are used for an independent evaluation of the posterior model states. TCCON is a network of Fourier transform spectrometers (FTSs) from near-infrared (NIR) solar absorption spectra, designed to retrieve precise total column abundances of CO₂, CH₄, N₂O, and CO to validate satellite observations (Wunch et al., 2011).

2.2 Inverse modeling

2.2.1 Variational inverse system

We use a Bayesian variational inversion system, PYVAR-LMDz, which uses LMDz-INCA as the chemistry transport model (CTM) (Hourdin et al., 2013; Hauglustaine, 2004). This inversion system has been documented and evaluated by a series of studies focusing on tracers including CH₄ (Pison et al., 2009; Locatelli et al., 2015; Cressot et al., 2014), HCHO (Fortems-Cheiney et al., 2012), CO (Yin et al., 2015; Zheng et al., 2019), and CO₂ (Chevallier et al., 2005, 2010). We use a recently developed version that has the capacity to assimilate observations of the major tracers in the CH₄ oxidation chain – namely, CH₄–HCHO–CO – with OH being their common sink path, to optimize the sources and sinks for all these species simultaneously. Here, we use a simplified chemistry scheme that assumes methane being oxidized into formaldehyde in a single step. We expect this simplification to have a relatively small impact on the inverse results of methane, given that all pathways of methane oxidation result in formaldehyde as an intermediate product. In addition, HCHO production from non-methane VOC oxidation is simulated upstream with a full-chemistry model, so that the correction on OH from the inversion will not directly feedback to the VOC oxidation. This should not be an issue as we optimize the production of HCHO instead of VOC emissions, but the impact of VOC on OH recycling is not accounted for. Future studies using a full chemistry scheme to optimize methane and OH simultaneously would be helpful to diagnose potential impacts of this simplification on the derived methane lifetime.

The CTM version we use here has a horizontal resolution of 1.875^∘ × 3.75^∘ (latitude, longitude) and a vertical resolution of 39 eta levels. Atmospheric transport is guided by the ERA-Interim meteorological reanalyses (Dee et al., 2011) to represent changes in the dynamics. Given observational information of the three species, the system optimizes the following quantities at each grid cell at a weekly resolution: (i) surface emissions of CH₄, (ii) surface emissions of CO, (iii) scaling factors for the sum of HCHO emissions and its chemical production from hydrocarbon oxidation, (iv) scaling factors of the OH concentration, and (v) the initial state of all the four species CH₄, HCHO, CO, and OH. The assimilation is performed continuously for the entire study period to avoid errors in temporal segmentation. The minimization of the cost function is solved iteratively until it reaches a reduction of 99 % in the gradient of the cost function or a minimum of 45 iterations. The reduced chi-squared (J divided by the number of observations) is about 0.5, which is much lower than 1 because of observation error inflation to compensate for the fact that we do not account for observation error correlations, following the findings of Chevallier (2007).

2.2.2 Prior estimates of surface methane fluxes and OH fields

We use prior estimates of climatological methane emissions from various sectors except for biomass burning. This choice is made to avoid prior assumptions about the interannual variations (IAVs) or trends in the surface emissions so that IAVs in the posterior fluxes are primarily driven by assimilated observations. The exception made for fire emissions is due to their non-Gaussian distribution and large variations across different seasons and years where the bottom-up estimates based on satellite-derived burned areas bring valuable prior information to guide the solution. The emission datasets from different sectors are listed in Table S2, and their spatial distributions are shown in Fig. S3. Note that soil deposition is treated as negative fluxes from the land to the atmosphere, and the emissions reported in this study are hence the net methane fluxes from the land to the atmosphere. The Gaussian uncertainty is set as 70 % and 100 %, respectively, for gridded CH₄ and CO emissions, whereas 200 % for chemical HCHO productions and 20 % for OH. Those errors are chosen empirically given the spreads across different bottom-up estimates. The a priori spatial error correlations are defined by an e-folding length of 500 km over the land and 1000 km over the ocean. Temporal error correlations are defined by an e-folding length of 2 weeks. We do not account for error correlations across species.

We include two alternative prior estimates for the OH concentration: one based on a full chemistry simulation by the model LMDZ-INCA (Hauglustaine, 2004), noted as INCA-OH hereafter, and one from the TransCom model intercomparison experiment for methane and related species (Patra et al., 2011), noted as TransCom-OH. The two OH fields have contrasting 3-D distributions that could help to evaluate the impact of OH distributions on the resultant methane fluxes (Yin et al., 2015). In particular, the two OH fields have different Northern to Southern hemisphere ratios: ∼ 1.2 for INCA and ∼ 1 for TransCom. Similar to the prior estimates of the emissions, there are no interannual variations in the prior estimates of OH fields. Note that for the case of assimilating surface observations (S1), the spatial error correlation of OH is set to 1 within six latitudinal zones (90–60^∘ S, 60–30^∘ S, 30–0^∘ S, 0–30^∘ N, 30–60^∘ N, and 60–90^∘ N) and 0 across them, i.e., the zonal mean OH is optimized instead of per grid cell given limited observational constraints. In summary, we include six inversions here with three different observational constraints and each pairing with two different prior estimates of global OH distributions (Table S1).

2.2.3 Information content analysis

While the variational inverse system has the advantage of optimizing large state vectors of fluxes for multiple species at high spatial and temporal resolutions, it is computationally too expensive to calculate the error covariances of posterior fluxes. Hence, we perform additional analytical inversions for aggregated source regions to estimate the information content of available methane observations on regional emissions and posterior error covariances. In this configuration, the state vector x becomes monthly regional emissions from 18 regions across the globe (regional mask shown in Fig. 8) plus a background term and the impacts of changes in OH are not accounted for. The transport model and observation operator K, relating each element of x to observable quantities y can be numerically simulated. Using x_a to represent the prior and S_a and S_ε to represent the error covariance matrices of the state vector x and of the observation vector y, the posterior solution is expressed as

$$\begin{array}{}\text{(1)}& \widehat{\mathit{x}}={\mathit{x}}_a+\mathbf{G}(\mathit{y}-\mathbf{K}{\mathit{x}}_a),\end{array}$$

where

$$\begin{array}{}\text{(2)}& \mathbf{G}={\mathbf{S}}_a{\mathbf{K}}^T({\mathbf{KS}}_a{\mathbf{K}}^T+{\mathbf{S}}_{\mathit{\epsilon }}{)}^{-\mathrm{1}}.\end{array}$$

Here, G represents the gain matrix that describes the sensitivity of the fluxes to observations, i.e., $\mathbf{G}=\partial \widehat{\mathit{x}}/\partial \mathit{y}$. The error covariance matrix $\widehat{\mathbf{S}}$ of $\widehat{\mathit{x}}$ can be derived as

$$\begin{array}{}\text{(3)}& \widehat{\mathbf{S}}=({\mathbf{K}}^T{\mathbf{S}}_{\mathit{\epsilon }}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_a^{-\mathrm{1}}{)}^{-\mathrm{1}}.\end{array}$$

The ability of an observational system to constrain the true value of the state vector can be represented by the sensitivity of the posterior solution $\widehat{\mathit{x}}$ to the true state x, commonly termed as the averaging kernel matrix $\mathbf{A}=\partial \widehat{\mathit{x}}/\partial \mathit{x}$, as the product of the gain matrix G and the Jacobian matrix $\mathbf{K}=\partial \mathit{y}/\partial \mathit{x}$, so that A=GK (Rodgers, 2000). This complementary analysis provides us important estimates of how much information content the surface and satellite methane observations can provide on regional methane emission changes.

3 Changes in the global CH₄ budget from 2010 to 2017

3.1 Changes in atmospheric methane growth rate

In general, the observed global average methane growth rate is well captured by the posterior model states both at the surface sites and through the total column, irrespective of which data are being assimilated (Fig. 2b and c and Table S4). Sampled from the same ensemble of posterior model states, the surface growth rates show a sharp increase in 2014 (Fig. 2b), whereas more gradual increase is found in the column average (Fig. 2c). The agreement across different inversions demonstrates that differences in the temporal variations of the growth rates seen by surface and GOSAT observations are primarily due to 3-D sampling differences rather than by some inconsistency between those two types of observations. This contrast suggests that the sharp increase in the surface methane growth rate in 2014 could have been amplified by a sampling effect of the sparse surface network as also shown by a longer record (Pandey et al., 2019). Surface in situ observations with high precision and accuracy provide critical anchor points for monitoring the background methane concentrations in the boundary layer, while satellite retrievals are sensitive to the entire atmospheric column filling in continental gaps that are not effectively covered by surface stations. The consistency between the two observation approaches demonstrates a robust constraint on the acceleration of the atmospheric growth rate on the global scale.

Figure 2 Global CH₄ emissions and atmospheric growth rates from 2010 to 2017. (a) Surface CH₄ fluxes of the prior (black triangles) and posterior estimates (color coded). The circles represent versions using INCA-OH (denoted with IN as suffix), while the squares represent versions using TransCom-OH (denoted with TR as suffix). The horizontal lines mark the average emissions of the two periods, 2010–2013 and 2014–2017. (b) Deseasonalized methane growth rates smoothed for variations shorter than 90 d in the posterior model states sampled at the 103 surface stations included in inversion S1. (c) $X_{{\mathrm{CH}}_{\mathrm{4}}}$ growth rates in the posterior model states sampled at the measurement time and location of GOSAT retrievals included in inversions S2 and S3.

3.2 Changes in global CH₄ emissions

Posterior global CH₄ emissions derived from all six inversions show similar interannual variations (IAVs) regardless of which observations are assimilated or which prior OH fields are used (Fig. 2a). As stated in the methods, the prior CH₄ emission IAVs only account for fire emissions, while the other emission sectors are represented by climatological means; hence, the IAVs of the posterior emissions are primarily driven by methane observations. Surface and satellite observations derive generally consistent IAV results. The choice of the prior OH fields has a notable effect on the magnitude of the optimized global emissions but not on the inferred temporal changes. Inversions using INCA-OH derives on average 20 ± 1.5 Tg yr⁻¹ higher emissions due to a larger OH sink (higher Northern Hemisphere OH concentrations). Therefore, in this study, we focus primarily on the IAV of methane fluxes that are directly relevant to changes in the methane growth rate while avoiding systematic differences across different inversions.

Global CH₄ emissions increased by 17.5 ± 1.5 Tg yr⁻¹ between 2010–2013 and 2014–2017 (the uncertainty range represents the standard deviation of the six inversions throughout this study). On average, the increase amounts to a linear trend of 4.1 ± 1.2 Tg yr⁻² over the 8 years, corresponding to a nearly 1% increase per year. The lowest annual total emission occurred in 2012 and the highest in 2017. Current global CH₄ emissions are thus at a maximum level within the past million years with high growth rates similar to the 1980s, during which the total methane sink was, however, not as high as today due to a lower CH₄ burden.

Figure 3 Global average posterior scaling factors on OH. Note that OH is only optimized by the system if other tracers in addition to CH₄ are assimilated (S1 and S3).

3.3 Variations attributed to OH

The inversion adjusted surface emission levels given the two different prior OH fields (Fig. 2a), indicating that there is not enough information to constrain the magnitudes of the sources and sinks of the three species separately with their atmospheric observations. Changes in the inferred OH concentrations are less than 1 % on the global scale, with a small increase during 2010–2014 followed by a small decline thereafter (Fig. 3). The resulting decrease in OH since 2014, albeit small in magnitude, occurs in both the surface-driven (S1) and satellite-driven (S3) inversions, most notably in the Southern Hemisphere (Fig. S4). Inflating the prior OH uncertainty up to ±50 % at each model grid only results in larger scaling factors on the OH distribution but not higher temporal variations. The resultant small interannual variations in the posterior OH field are in line with a modeling study that showed a high OH recycling probability and hence a weak sensitivity to emission perturbations (Lelieveld et al., 2016). Some atmospheric chemistry models simulate a slightly larger year-to-year variability (1 %–4 %) (Holmes et al., 2013; Turner et al., 2018), while recent data-constrained estimates using observed ozone columns, water vapor, methane, model-simulated NO_x, and Hadley cell width suggest a relatively stable OH level over the past several decades (Nicely et al., 2018). In addition, compared to earlier box model studies that infer around 5 % OH IAV from methyl chloroform (MCF) and δ¹³CH₄ observations (Turner et al., 2017; Rigby et al., 2017), a recent box model study that accounts for model biases related to tracer specific dynamics suggest a smaller IAV in OH (Naus et al., 2019).

Recent GOSAT inverse studies explored optimizing gridded annual anthropogenic methane emissions and associated trends, regional monthly wetland emissions, and global (or hemispheric) annual OH concentration with an analytical inversion scheme, where it is possible to compute the full posterior error covariance matrix (Maasakkers et al., 2019; Zhang et al., 2021). The results suggest a strong negative error correlation between global anthropogenic emissions and methane lifetime ($r=-\mathrm{0.8}$), moderate correlations between wetland emissions and methane lifetime ($r=-\mathrm{0.4}$), and moderate correlations between the OH trend and wetland or anthropogenic emission trends ($r=-\mathrm{0.6}$) (Zhang et al., 2021). Hence, assimilating GOSAT data alone, the inversion has limited information to separate the sources and the sinks. With our multi-species variational inverse system, it is computationally too costly to estimate the posterior error covariances using a Monte Carlo approach. Given the strong error correlations between the source and sink terms identified by Zhang et al. (2021), we cannot rule out the possibility that numerically it might be easier for the optimization system to adjust surface emissions of the three species to fit the observations rather than modifying OH to adjust the sink terms in the absence of mechanistic chemical feedback in the chemical transport model. The feedback effects are mostly tested using box models at the current stage (Prather, 1994; Nguyen et al., 2020); future studies accounting for these effects in a 3-D inversion would be helpful to diagnose its impacts on estimated changes in methane lifetimes.

4 Regional contributions

4.1 Changes in zonal CH₄ emissions

Similar zonal emission increases between 2010–2013 and 2014–2017 are found across the six inversions (Fig. 4a), even though they produce different latitudinal distributions of CH₄ fluxes (Fig. 4b). Both satellite and surface data suggest that the largest increase occurred in the southern tropics (0–30^∘ S, 7.5 ± 2.1 Tg yr⁻¹) and the northern mid-latitudes (30–60^∘ N, 6.5 ± 0.8 Tg yr⁻¹), while a moderate increase is found in the northern high latitudes (60–90^∘ N, 1.3 ± 0.5 Tg yr⁻¹). For the northern tropics (0–30^∘ N), most versions suggest a small increase, but one version assimilating surface data suggests a small decline. Different versions agree on the overall spatial distribution of the inferred emission trends, with the most significant increase seen in East China, tropical South America, tropical Africa, and Russia (Fig. 5). Opposing trends are noted in former Indochina and Southeast Asia that result in more divergent estimates across the different inversions in the 0–30^∘ N zone.

Figure 4 (a) Emission change between 2010–2013 and 2014–2017 in the five latitudinal zones. The error bars represent the standard deviation of changes in CH₄ fluxes between the two periods. (b) Zonal fluxes estimated by different versions for the period 2010–2013 and 2014–2017. The mean values for each 4-year period are shown and errors bars represent their 1σ standard deviations.

Differences of zonal flux distributions are noted across versions, most notably between surface and satellite data constraints. For the same observational constraints, inversions using INCA OH fields result in higher Northern Hemisphere emissions compared to the cases using TransCom OH fields due to a higher north-to-south hemispheric OH ratio of the former. Compared to the results assimilating surface observations (S1), assimilating GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ retrievals (S2 and S3) allocates smaller emissions in the northern mid- and high latitudes (30–60 and 60–90^∘ N) but higher emissions in the tropics and subtropics (0–30^∘ N and 0–30^∘ S) (Fig. 4b). Such difference is, to a large extent, related to a latitudinal-dependent difference between model states that fit surface data and that fit GOSAT data. Specifically, the posterior model states of S1 that fit surface observations show positive biases against GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ in the northern mid- and high latitudes but negative ones in the tropics (Fig. S5). Symmetrically, the posterior model states of S2 and S3, which fit GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ well, show negative biases in the northern mid- and high latitudes against surface observations (Fig. S7), while the biases turn positive gradually toward the tropics and the Southern Hemisphere (Fig. S16). However, no latitude-dependent biases are found between GOSAT-assimilated posterior model states (S2 and S3) against TCCON total column measurements, and the magnitude of remaining biases are in line with GOSAT data validation (Parker et al., 2015). Yet S1 shows similar model bias structure against TCCON as compared to GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ (Fig. S7), suggesting discrepancies in the vertical distribution of methane concentrations between the model and the satellite retrievals. Such a bias pattern between model and surface or GOSAT data has been identified by previous inverse studies (Alexe et al., 2015; Turner et al., 2015; Miller et al., 2019; Maasakkers et al., 2019), which is likely related to biases in the model representation of the stratosphere. An empirical bias correction on the GOSAT data so that the assimilated model states also agree with surface observations are typically applied by some studies. Here, since we focus on the IAVs of the posterior fluxes where systematic biases do not impact such results, we did not apply an empirical bias correction to the GOSAT data. Future studies to correct those biases with mechanistic understandings will be very valuable.

Figure 5 Spatial distribution of trends in the posterior CH₄ emissions from 2010 to 2017. Panels (a), (c), and (e) show results using INCA-OH and panels ((b), (d), (f) use TransCom-OH. Each row represents one type of observational constraint. The black crosses denote trends that are statistically significant at a 95 % confidence level.

4.2 Information content of observations on regional fluxes

To assess the extent to which the surface and satellite observations can inform us about changes of methane fluxes in distinct regions, we conducted an information content analysis for a total of 18 regions (see Sect. 2.2.3). The regional mask following the convention of the Global Carbon Project (Saunois et al., 2019) is shown in Fig. 8. Note that this analysis assumes all atmospheric methane changes result from surface flux changes and hence does not account for potential contributions from changes in OH or other sink processes. The results suggest that in most cases GOSAT data provide more constraints on regional emissions than the surface observations (Fig. 6). This is particularly obvious in the tropics and subtropics, including the Amazon, eastern Brazil, southern South America, northern Africa, tropical Africa, southern Africa, the Mideast, India, and Southeast Asia. This is because fewer surface sites exist in those regions but satellite data have a better coverage. Consequently, the posterior errors in the optimized emissions constrained by satellite data are less correlated across different regions compared to the case with surface data constraints only (Fig. S8). The error covariances suggest that the surface observations alone, mostly located in the background boundary layer, is insufficient to separate tropical emissions from the three continents – South America, Africa, and Asia. In contrast, the cross-error terms in the GOSAT inversion are much smaller, suggesting that to a large extent emissions from different regions can be individually constrained by these X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ observations.

Figure 6 Averaging kernels (AKs) of regional emissions to observations over that region during each month. Diagonal terms of the AK matrix are shown using data from 2010 as an example.

4.3 Regional emission changes

Looking into regional emission changes, the differences in the posterior CH₄ emissions between the last and the first 4 years of our study period (2014–2017 vs. 2010–2013) are shown in Fig. 7, while regional masks of the 18 sub-regions and regional annual emission anomalies are shown in Fig. 8. The most substantial increases between the two periods occurred in the Amazon, China, and tropical Africa, by 4.2 ± 1.2, 3.7 ± 1.0, and 2.1 ± 0.8 Tg yr⁻¹, respectively (Fig. 7). Changes in the three regions amounts to nearly 60 % of the global emission increase. This increase does not necessarily imply linear trends in emissions as there are considerable interannual variations in the derived emissions (Fig. 8); in particular, the first period includes a La Niña year (2011) during which high wetland emissions are reported (Pandey et al., 2017) and the latter period includes a strong El Niño year (2015) during which large fire emissions are reported (Yin et al., 2016; Worden et al., 2017). While all the six inversions agree on such a regional pattern, the multi-tracer versions (S3) that optimize OH concentrations simultaneously with the surface methane fluxes infer smaller CH₄ emission increases compared to the version assimilating GOSAT alone (S2). This difference could stem from adapting the regional mean OH level that converts the same concentration change to different emission changes. In addition, differences between S2 and S3 could result from the variational inversion reaching different approximations of the cost function minimum.

Figure 7 Regional emission changes between 2010–2013 and 2014–2017 ranked from the highest to the smallest changes. The color-coded markers represent individual inversions, the gray stars represent the ensemble mean, and the horizontal error bars denote the standard deviation of all versions. The regional mask is shown in Fig. 8.

To gain further understanding of observed changes in regional CH₄ emissions, we attribute our inversion emission anomaly estimates into the following categories based on our prior bottom-up emission inventory: fossil fuel (oil, gas, coal mining, industry, residential, transport, and geological), waste (landfills and wastewater), agriculture (enteric fermentation, manure management, and rice cultivation), wetlands (including inland water), and fire (including biofuel). We acknowledge the fact that this prior information has significant uncertainties as evidenced by the large spread across different bottom-up inventories (Saunois et al., 2016). The proportion of the different sectors remains unchanged in each grid cell throughout all years, except for fire, because we use a climatological estimates for prior emissions. Our emission attribution thus reflects a likelihood of contributing processes at a given location and season, which is larger and most useful in regions where emissions are predominately contributed by a specific sector (Figs. S9 and S10).

For the Amazon, wetlands are the major contributor to CH₄ emissions according to the bottom-up emission inventories and hence our identified source for the increase, showing an average trend of 0.8 ± 0.1 Tg yr⁻² over the 8 study years with shorter-term interannual variations (Fig. 8). Fire emissions from this region were high during the 2010 drought but did not rise significantly in the recent 2015 El Niño, which is in agreement with previous fire emission estimates based on CO and CO₂ (Gatti et al., 2014; Liu et al., 2017). No significant trend in the anthropogenic emissions are noted up to 2014 according to the most recent updates from the Community Emissions Data System (Hoesly et al., 2018) (Fig. S11). Our inferred wetland emissions in the 2011 La Niña show the highest positive anomaly in the 2010–2013 period, consistent with previous estimates covering this particular period (Pandey et al., 2017). Wetland methane emissions come from anaerobic degradation of organic matter and hence depend on organic carbon inputs and inundation areas and exponentially on temperature (Whalen, 2005). Consistent behaviors between the time and locations of anomalies in the GOSAT X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ and changes in wetland extent have been documented with the focus on seasonally flooded wetlands (Parker et al., 2018). An intensification of Amazon flooding extremes has been documented based on water levels in the Amazon river, with anomalously high flood levels and long flood durations since 2012 (Barichivich et al., 2018), which could result in higher wetland CH₄ emissions.

Figure 8 Regional emission anomalies relative to the 2010–2013 mean for the sectoral attribution based on prior information. Note the scales on the y axis are different for each subplot.

For the other tropical regions, significant increases are also attributed to wetland emissions, in particular to tropical Africa (1.5 ± 0.7 Tg yr⁻¹; Fig. 8). The increasing tropical Africa wetland emissions are consistent with a recent regional inversion using GOSAT data at a high spatial resolution of 0.5^∘ × 0.625^∘, which finds a positive trend of 1.5–2.1 Tg yr⁻² in the region from 2010 to 2016, mainly from wetlands in the Sudd in South Sudan (Lunt et al., 2019). Smaller wetland emission increases are found in the other tropical regions including eastern Brazil (0.3 ± 0.1 Tg yr⁻¹), northern Africa (0.2 ± 0.1 Tg yr⁻¹), and southern South America (0.1 ± 0.1 Tg yr⁻¹). However, other emission sources also play a significant role in these regions, in particular agricultural emissions (Chang et al., 2019). Thus future studies with additional constraints on wetland emissions are needed to better quantify wetland-related changes. Only in Southeast Asia, is the major contribution to different CH₄ emissions between the two periods from fire associated with the strong El Niño in 2015 (Yin et al., 2016; Liu et al., 2017). No significant increases are noted for India, consistent with a previous regional study focusing on the 2010–2015 period (Ganesan et al., 2017).

The sectoral breakdown of emissions from China suggests a substantial increase in anthropogenic sources from fossil fuels, agriculture, and waste, adding up to an overall trend of 1.0 ± 0.2 Tg yr⁻² between 2010 and 2017 (Fig. 8). As stated above, this attribution relies on the relative contribution of different sectors from the prior information and does not account for structural changes in time. A recent inverse study focusing on Asian emissions from 2010 to 2015 derived nearly the same magnitude of emission trend for China (Miller et al., 2019); a continued increase is confirmed here beyond 2015 until the end of the record in 2017. In contrast, a global inversion that uses a different prior emissions estimate and separates the mean anthropogenic emissions and trends in the state vector found a smaller trend in anthropogenic emissions over China (0.39 ± 0.27 Tg yr⁻²) for the period 2010–2018 and a trend of 0.72 ± 0.39 Tg yr⁻² focusing on the period 2010–2016 (Zhang et al., 2021). The numbers are comparable given the differences in the inverse setups and the chemical transport models being used.

Russia also contributed significant increase in CH₄ emissions by 1.7 ± 0.7 Tg yr⁻¹ between 2010–2013 and 2014–2017 (Fig. 7), possibly from both fossil fuel extraction in northern Russia and extensive peatland areas (Fig. 8). The surface-driven and satellite-driven inversions identify slightly different source regions for the rise (Fig. 5). The surface-driven inversions attribute most of the increases to the European part of Russia where anthropogenic emissions dominate, whereas the satellite-driven inversions attribute more changes to the West Siberian Plain where more wetlands are located (Terentieva et al., 2016). As there are both fossil fuel and wetland sources in the West Siberian Plain (Fig. S10), further information is needed to disentangle relative contributions between anthropogenic and natural wetland sources. For the other extratropical regions showing significant CH₄ emission increases, the increase in Canada (1.1 ± 0.4 Tg yr⁻¹) was mostly attributed to wetlands (Fig. 8), with interannual variations consistent with previous regional inversions (Sheng et al., 2018). A relatively small increase in the US is found after 2014 (0.7 ± 0.2 Tg yr⁻¹) with nearly flat emissions before, which is consistent with previous studies finding no trend over the US before 2012 (Saunois et al., 2017; Bruhwiler et al., 2017).

Relying on the prior distribution to approximate possible contributions from wetlands in the mid- and high latitudes, the increase between 2010–2013 and 2014–2017 amounts to 0.9 ± 0.5, 0.6 ± 0.4, and 0.1 ± 0.06 Tg yr⁻¹ for Russia, Canada, and the US. Up to 2012, high-latitude wetland emissions are not identified as significant contributors to increasing atmospheric methane (Saunois et al., 2017). The positive trend in high-latitude wetland emissions found here could be the first sign of the impact of the fast warming observed at these latitudes. Adding up all wetland contributions across the globe, changes in wetland emissions dominate the interannual variations in the emission anomaly (Fig. S12a). The general increase in wetland CH₄ fluxes is in line with observed atmospheric δ¹³CH₄ that shows a general negative trend at all latitudes (Fig. S12b), as biogenic sources like wetlands are more δ¹³CH₄ depleted than the other ones (Sherwood et al., 2017). In addition, anthropogenic emissions from agriculture and waste management are also associated with a biogenic δ¹³CH₄ signature.

5 Conclusions

Our ensemble of inversions assimilating surface or satellite CH₄ observations, as well as chemically related species to partly constrain the OH sink, suggests that the recent acceleration in CH₄ growth rate from 2010 to 2017 is most likely induced by increases in methane emissions. The derived global emissions point to an unprecedented new maximum in global total methane emissions. The most substantial increases during the 8 study years come from the tropics and East Asia. Given our prior knowledge on the distribution of different CH₄ sources, natural wetland emissions show the largest increase with dominant contributions from the tropics. Such an increase would result in potential positive feedback to climate warming (Zhang et al., 2017). The second-largest increase comes from anthropogenic emissions in China. The continuation of existing surface CH₄ and δ¹³CH₄ observations and GOSAT/GOSAT-2 X${}_{{\mathrm{CH}}_{\mathrm{4}}}$ retrievals, the newly available TROPOspheric Monitoring Instrument (TROPOMI) observations with frequent global mapping capability (Hu et al., 2018), and the new methane space missions such as the MEthane Remote sensing Lidar missioN (MERLIN) (Bousquet et al., 2018) will bring further insight into regional methane budget changes and their climate sensitivity. Here, we tested the consistency of using different observational constraints and different prior OH distributions. The sensitivity of prior emission estimates, associated error characteristics, and transport model errors could be further explored by model intercomparison studies through community effort such as Patra et al. (2011). Future studies using spatial–temporal variations in the observed atmospheric δ¹³CH₄ and spatially resolved isotopic source signatures (Ganesan et al., 2018) will provide further constraints on the source attribution. At the same time, a process-based understanding of the wetland CH₄ emissions and effective anthropogenic emission regulation measures are urgently needed to meet climate mitigation goals.