Investigation of the renewed methane growth post-2007 with high-resolution 3-D variational inverse modeling and isotopic constraints

We investigate the causes of the renewed growth of atmospheric methane (CH₄) amount fractions after 2007 by using variational inverse modeling with a three-dimensional chemistry-transport model. Together with CH₄ amount fraction data, we use the additional information provided by observations of CH₄ isotopic compositions (¹³C : ¹²C and D : H) to better differentiate between the emission categories compared to the differentiation achieved by assimilating CH₄ amount fractions alone. Our system allows us to optimize either the CH₄ emissions only or both the emissions and the source isotopic signatures (δ_source(¹³C,CH₄) and δ_source(D,CH₄)) of five emission categories. Consequently, we also assess, for the first time, the influence of applying random errors to both emissions and source signatures in an inversion framework. As the computational cost of a single inversion is high at present, the methodology applied to prescribe source signature uncertainties is simple, so it can serve as a basis for future work. Here, we investigate the post-2007 increase in atmospheric CH₄ using the differences between 2002–2007 and 2007–2014. When random uncertainties in source isotopic signatures are accounted for, our results suggest that the post-2007 increase (here defined using the two periods 2002–2007 and 2007–2014) in atmospheric CH₄ was caused by increases in emissions from (1) fossil sources (51 % of the net increase in emissions) and (2) agriculture and waste sources (49 %), which were slightly compensated for by a small decrease in biofuel- and biomass-burning emissions. The conclusions are very similar when assimilating CH₄ amount fractions alone, suggesting either that random uncertainties in source signatures are too large at present to impose any additional constraint on the inversion problem or that we overestimate these uncertainties in our setups. On the other hand, if the source isotopic signatures are considered to be perfectly known (i.e., ignoring their uncertainties), the relative contributions of the different emission categories are significantly changed. Compared to the inversion where random uncertainties are accounted for, fossil emissions and biofuel- and biomass-burning emissions are increased by 24 % and 41 %, respectively, on average over 2002–2014. Wetland emissions and agricultural and waste emissions are decreased by 14 % and 7 %, respectively. Also, in this case, our results suggest that the increase in CH₄ amount fractions after 2007 (despite a large decrease in biofuel- and biomass-burning emissions) was caused by increases in emissions from (1) fossil fuels (46 %), (2) agriculture and waste (37 %), and (3) wetlands (17 %). Additionally, some other sensitivity tests have been performed. While the prescribed interannual variability in OH can have a large impact on the results, assimilating δ(D,CH₄) observations in addition to the other constraints has only a minor influence. Using all the information derived from these tests, the net increase in emissions is still primarily attributed to fossil sources (50 ± 3 %) and agriculture and waste sources (47 ± 5 %). Although our methods have room for improvement, these results illustrate the full capacity of our inversion framework, which can be used to consistently account for random uncertainties in both emissions and source signatures.

1 Introduction

Atmospheric methane (CH₄) has a large influence on both climate and atmospheric chemistry. The globally averaged tropospheric CH₄ amount fraction has increased by a factor of 2.6 since pre-industrial times (Gulev et al., 2021), and it reached a new high of 1912 nmol mol⁻¹ in 2022 (global average from marine surface sites;  Lan et al., 2023). Neglecting indirect effects related to ozone, water vapor and nitrogen oxide production, a consequence of this large increase in CH₄ amount fractions since the pre-industrial era is that CH₄ now contributes 16 % of the current radiative forcing from well-mixed greenhouse gases (carbon dioxide, methane, nitrous oxide, halogens) (Forster et al., 2021). CH₄ therefore makes the second largest contribution to the additional greenhouse effect, behind carbon dioxide (CO₂). CH₄ amount fractions increased quasi-continuously since the pre-industrial era until 1999, but then stabilized between 1999 and 2006. The growth resumed after 2007, at a rate exceeding 10 $\mathrm{nmol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ for some years. Nisbet et al. (2019) pointed out that the dramatic increase in the CH₄ burden is contrary to pathways compatible with the goals of the Paris Agreement of the 2015 United Nations Framework Convention on Climate Change and that urgent action is required to bring CH₄ back to a pathway more in line with the goals of the Paris Agreement. A proper understanding of the CH₄ budget could highly facilitate such actions by increasing the effectiveness of mitigation policies.

CH₄ is emitted into the atmosphere by multiple sources (wetlands, livestock, rice cultivation, waste, fossil fuel exploitation, biomass burning, …), with distinct processes involved (microbial, thermogenic, pyrogenic). This species is mainly removed from the atmosphere through oxidation by the hydroxyl radical (OH), which represents about 92 % of the total sink (Saunois et al., 2020; Thanwerdas et al., 2022b). Other sinks include oxidation by atomic oxygen (O¹D), chlorine (Cl) and methanotrophs in the soil, which contribute about 1.5 %, 1.5 % and 5 %, respectively, to the total removal of CH₄ (Saunois et al., 2020; Thanwerdas et al., 2022b). Note that these numbers come with non-negligible uncertainties and vary from one study to another.

Estimating these sources and sinks is challenging, especially at the global scale, yet it is necessary to better understand the CH₄ budget and to anticipate its evolution. The scientific community has developed two approaches to estimate CH₄ emissions at different scales. On the one hand, bottom-up approaches aim to estimate these emissions using both inventories that combine statistical activity data with emission factors for anthropogenic emissions (e.g., Höglund-Isaksson, 2012, 2017; Janssens-Maenhout et al., 2019) and process-based models for natural and fire emissions (e.g., van der Werf et al., 2017; Poulter et al., 2017). Bottom-up estimates provide valuable sectorial and regional information, although their global emissions are not constrained by atmospheric observations. On the other hand, top-down approaches use inversion methods (Newsam and Enting, 1988; Enting and Newsam, 1990) and chemistry-transport models (CTMs) to statistically optimize model parameters (e.g., emissions) and minimize model–observation differences (e.g., Houweling et al., 2017, and references therein). These approaches provide posterior estimates that are consistent with both atmospheric observations (e.g., CH₄ amount fractions) and prior estimates (typically derived from bottom-up estimates). The inversion problem is considered to be “ill posed” because a wide range of surface flux configurations can equally well explain the observational data in the atmosphere. Since the 1980s, surface monitoring networks have nevertheless significantly increased the spatial coverage and the precision of their observations, narrowing the range of possible flux configurations and improving the relevance of inversion methods.

Although Saunois et al. (2020) recently showed that the consistency between top-down and bottom-up estimates has improved over time, the 1999–2006 plateau and the subsequent renewed growth still generate considerable attention and controversy (Rice et al., 2016; Schwietzke et al., 2016; Schaefer et al., 2016; Nisbet et al., 2016; Patra et al., 2016; Bader et al., 2017; Turner et al., 2017; Rigby et al., 2017; Worden et al., 2017; Saunois et al., 2017; Morimoto et al., 2017; McNorton et al., 2018; Thompson et al., 2018; Schaefer, 2019; Nisbet et al., 2019; Fujita et al., 2020; Zimmermann et al., 2020; Jackson et al., 2020; Chandra et al., 2021). Most of these studies suggest that the renewed growth is partially explained by an increase in microbial emissions (wetlands, livestock and/or rice cultivation), while some of them localized the increase to the tropics (Nisbet et al., 2016; Patra et al., 2016; Schaefer et al., 2016; Schwietzke et al., 2016). Multiple studies have also concluded that the renewed growth was driven by an increase in both microbial and fossil fuel emissions (Rice et al., 2016; Patra et al., 2016; Bader et al., 2017; Worden et al., 2017; Saunois et al., 2017; McNorton et al., 2018; Thompson et al., 2018; Jackson et al., 2020; Chandra et al., 2021; Lan et al., 2021; Basu et al., 2022), although they provide a very wide range of individual contributions. An increase in fossil fuel emissions was also supported by an independent work that used ethane-based approaches (Hausmann et al., 2016). However, other studies found that these emissions decreased or stabilized (Schwietzke et al., 2016; Schaefer et al., 2016; Fujita et al., 2020) over the period of renewed growth. Some studies also found that an increase in emissions was not the main driver and that a large decrease in OH concentrations could explain the recent variations (Turner et al., 2017; Rigby et al., 2017).

Such controversy partly arises from the difficulty in separating contributions from individual CH₄ sources. Despite the high number of observations over some regions, many of the sources are co-located, and isolating the contribution from each source to the local increase in CH₄ amount fractions is challenging. Atmospheric carbon and hydrogen isotope compositions, δ(¹³C,CH₄) and δ(D,CH₄), can help us to differentiate co-emitted emission categories because each CH₄ production process (microbial, thermogenic, pyrogenic) has its own characteristic isotopic signature (Sherwood et al., 2021, 2017). δ(¹³C,CH₄) and δ(D,CH₄) are generally defined as the deviation of the sample's atomic isotopic ratio ($R_{\mathrm{13}}=X{(}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}})/X{(}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}})$ or $R_{\mathrm{D}}=X\left({\mathrm{CH}}_{\mathrm{3}}\mathrm{D}\right)/X\left({\mathrm{CH}}_{\mathrm{4}}\right))$ relative to a specific standard ratio:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\delta }\left({}^{\mathrm{13}}\mathrm{C},{\mathrm{CH}}_{\mathrm{4}}\right)={\displaystyle \frac{R_{\mathrm{13}}}{R_{\mathrm{PDB}}}}-\mathrm{1}={\displaystyle \frac{X\left({}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\right)/X\left({}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}\right)}{R_{\mathrm{PDB}}}}-\mathrm{1}\text{(2)}& {\displaystyle }\mathit{\delta }\left(\mathrm{D},{\mathrm{CH}}_{\mathrm{4}}\right)={\displaystyle \frac{R_D}{R_{\mathrm{VSMOW}}}}-\mathrm{1}={\displaystyle \frac{X\left({\mathrm{CH}}_{\mathrm{3}}\mathrm{D}\right)/X\left({\mathrm{CH}}_{\mathrm{4}}\right)}{R_{\mathrm{VSMOW}}}}-\mathrm{1}.\end{array}$$

Here, X(¹²CH₄), X(¹³CH₄), X(CH₃D) and X(CH₄) denote the ¹²CH₄, ¹³CH₄, CH₃D and CH₄ amount fractions, respectively. $R_{\mathrm{PDB}}=\mathrm{1.12372}\times {\mathrm{10}}^{-\mathrm{2}}$ is the standard ratio of Pee Dee Belemnite (PDB) (Craig, 1957) and $R_{\mathrm{VSMOW}}=\mathrm{1.5595}\times {\mathrm{10}}^{-\mathrm{4}}$ is the Vienna Standard Mean Ocean Water (VSMOW) ratio (Hagemann et al., 1970; Wit et al., 1980). δ(¹³C,CH₄) and δ(D,CH₄) are expressed in ‰. Broadly summarized, CH₄ sources have a δ(¹³C,CH₄) isotopic signature, hereinafter denoted by δ_source(¹³C,CH₄), of between −65 ‰ and −55 ‰ for microbial sources, between −45 ‰ and −35 ‰ for thermogenic sources, and between −25 ‰ and −15 ‰ for pyrogenic sources (Sherwood et al., 2021, 2017). However, the full distributions of isotopic signatures are wider than these ranges, with overlaps between the distributions for the different production processes. Similarly, there are different δ_source(D,CH₄) distributions for different production processes. Their approximate ranges are between −350 ‰ and −100 ‰ for thermogenic sources, between −400 ‰ and −250 ‰ for microbial sources, and between −250 ‰ and −175 ‰ for pyrogenic sources (Sherwood et al., 2021, 2017). Notably, the microbial and thermogenic δ_source(D,CH₄) distributions show a smaller overlap than that of the corresponding distributions for δ_source(¹³C,CH₄), and thermogenic sources have signatures that are less distinguishable from others.

Variations in atmospheric isotopic composition are not caused by sources only. Reactions between sink species (OH, O¹D and Cl) and CH₄ have rates that depend on the isotopologue. This effect is called fractionation and is represented, for a specific reaction, using the ratio of the reaction rates achieved with the lightest and the heaviest members of a couple of isotopologues (e.g., ¹²CH₄ and ¹³CH₄). The fractionation effect explains why the atmospheric isotopic composition is not equal to the flux-weighted mean source signature for all the CH₄ sources. It acts to shift this mean source composition towards less negative values when CH₄ enters the atmosphere and gets removed by the sinks. This effect is particularly important for δ_source(D,CH₄) because the flux-weighted mean source signature and the observed isotopic composition are approximately −330 ‰ and −95 ‰, respectively (Sherwood et al., 2017). For δ(¹³C,CH₄), this effect is smaller, shifting the source signature from approximately −53.6 ‰ to −47.3 ‰ in the atmosphere (Sherwood et al., 2017).

The post-2007 CH₄ increase is notably associated with a decrease of 0.2 ‰–0.3 ‰ in δ(¹³C,CH₄) since 2007 (Rice et al., 2016; Schaefer et al., 2016; Nisbet et al., 2016, 2019). Such significant isotopic variations provide an additional atmospheric constraint to better estimate the relative contribution of CH₄ sources to this renewed atmospheric CH₄ growth. Some of the aforementioned studies that focused on the drivers of both the plateau and the renewed growth were conducted using inversion methods that included isotopic constraints. Those studies implemented either three-dimensional (3-D) CTMs coupled with analytical inversion methods that estimated emissions for aggregated large regions only (Rice et al., 2016; McNorton et al., 2018), box models with analytical inversion methods (Schwietzke et al., 2016; Turner et al., 2017; Rigby et al., 2017), or 2-D CTMs with variational inversion methods (Thompson et al., 2018). Analytical methods are not fit to use for large-dimension problems, i.e., those with both a large number of optimized variables and observations, and these methods generally need to aggregate emissions across large regions. By contrast, variational inversion methods can easily both optimize the emissions at the grid-cell scale (the model's horizontal resolution) and assimilate large observational datasets. Furthermore, 3-D CTMs can better capture the spatial variability of sources, sinks and observations than box models and 2-D CTMs.

This paper utilizes the system designed by Thanwerdas et al. (2022a) to investigate changes in CH₄ emissions from 1998–2018 by running 3-D variational inversions at the grid-cell scale. The original system has been improved and can assimilate both ¹³CH₄ : ¹²CH₄ and CH₃D : CH₄ observations. The optimization of source isotopic signatures is also tested here because, at present, they remain a large source of uncertainty (Sherwood et al., 2017; Turner et al., 2017; Feinberg et al., 2018) that should be considered.

In Sect. 2, we provide a detailed methodology describing the inversions performed in this study. In Sect. 3, the results are presented. First, we evaluate the agreement between model outputs and assimilated data and compare our simulations to independent data. As a second step, we provide an analysis of posterior emissions and isotopic signatures estimated by the reference inversion and the sensitivity tests. To the best of our knowledge, the methodology developed by Basu et al. (2022) is the only one that presents high similarities to ours. They investigated the same problem with a variational inversion framework and a 3-D CTM. However, substantial differences exist between our techniques. In their paper, they included a comparison between their work and Thanwerdas et al. (2022a). Based on our new results, we propose an updated comparison in Sect. 3.9. Conclusions are drawn and a discussion is provided in Sect. 4.

2 Methods

2.1 The chemistry-transport model

The general circulation model (GCM) LMDz is the atmospheric component of the coupled model from the Institut Pierre-Simon Laplace (IPSL-CM), developed at the Laboratoire de Météorologie Dynamique (LMD) (Hourdin et al., 2006). The version of LMDz used here is an offline version dedicated to the inversion framework created by Chevallier et al. (2005): the precalculated meteorological fields provided by the online version of LMDz are given as input to the model, considerably reducing the computation time. The model is built at a horizontal resolution of 3.8^∘ × 1.9^∘ (96 grid cells in longitude and latitude), with 39 hybrid sigma-pressure levels reaching an altitude of about 75 km. The time step of the model is 30 min and the output values have a resolution of 3 h. Horizontal winds are nudged towards the ECMWF meteorological analyses (ERA-Interim) in the online version of the model. Vertical diffusion is parameterized by the local approach of Louis (1979), and deep convection processes are parameterized by the scheme of Tiedtke (1989). The offline model LMDz is coupled with the Simplified Atmospheric Chemistry System (SACS) to represent CH₄ oxidation by radicals (Pison et al., 2009; Thanwerdas et al., 2022a).

We simulate atmospheric ¹²CH₄ and ¹³CH₄ amount fractions to retrieve both the CH₄ amount fractions and the δ(¹³C,CH₄) signal. Four clumped isotopologues (¹²CH₄, ¹²CH₃D, ¹³CH₄ and ¹³CH₃D) are simulated in one sensitivity simulation (see Sect. 2.7) to retrieve both the δ(¹³C,CH₄) and δ(D,CH₄) compositions.

Oxidations by OH, O(¹D) and Cl are included in the chemical scheme of LMDz-SACS. Time-varying 3-D fields of OH and O(¹D) with daily resolution, simulated beforehand with the LMDz-INCA chemistry model (Hauglustaine et al., 2004), are prescribed for each oxidant species to simulate the associated chemical loss. The same meteorological data were used for generating these fields and running the simulations presented in this study.

The resulting OH field, named OH-INCA, exhibits a global mean tropospheric mass-weighted concentration of 11.1 × 10⁵ cm⁻³ over 1998–2018, consistent with the previous estimates from Zhao et al. (2019) (11.7 × 10⁵ cm⁻³) and Prather et al. (2012) ((11.2 ± 1.3) × 10⁵ cm⁻³) and well within the range derived from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) (10.3–13.4 × 10⁵ cm⁻³; Voulgarakis et al., 2013). It is, however, slightly larger than the very recent estimate from Zhao et al. (2023), obtained by constraining OH with observations of its precursors. The interhemispheric ratio is 1.14, lower than the mean value of 1.3 inferred by Zhao et al. (2019) but more consistent with recent estimates from Zhao et al. (2023) and an interhemispheric parity obtained from methyl-chloroform-based inversions (Bousquet et al., 2005; Patra et al., 2014). Global concentrations of OH-INCA increased by 4 % between 2002 and 2014.

As suggested by Thanwerdas et al. (2022b), the Cl concentrations derived by Wang et al. (2021) are prescribed here for all simulations. Their work suggests that the Cl sink accounts for only 0.8 % of the total CH₄ oxidation, which is lower than other estimates used in the literature (1.8 %–5 %; Allan et al., 2007; Sherwen et al., 2016; Hossaini et al., 2016).

The fractionation effect must also be represented in the modeling framework. Table 1 provides the fractionation coefficients applied for each loss reaction. For the OH sink, we adopted the estimate derived by Saueressig et al. (2001). Burkholder (2019) recommends using the Saueressig et al. (2001) rates but suggests increasing the uncertainty in the OH fractionation to account for the estimate (1.0054) from Cantrell et al. (1990). As shown by Basu et al. (2022), switching estimates from that of Saueressig et al. (2001) to that of Cantrell et al. (1990) has a large influence on the results, despite the fact that those authors did not optimize source signatures in their setup. As Saueressig et al. (2001) indicate that their data is of considerably higher experimental precision and reproducibility than that from previous studies, in particular Cantrell et al. (1990), we prefer to allocate computational time to a sensitivity inversion testing a different OH field rather than testing a different OH fractionation coefficient. In addition, these estimates of fractionation coefficients come with uncertainty ranges that we could also consider in our inversions (e.g., with a Monte Carlo approach). In our case, the main limitation remains the large computational cost of one inversion (see Sect. 2.9). In the future, we hope to be able to increase the number of sensitivity tests and account for this uncertainty. For this work, the values we adopt are the best estimates for each fractionation coefficient.

Saueressig et al. (2001)Saueressig et al. (2001)Saueressig et al. (2001)Saueressig et al. (2001)Saueressig et al. (1995)Saueressig et al. (1996)Snover and Quay (2000)Snover and Quay (2000)Reeburgh et al. (1997)Tyler et al. (1994)King et al. (1989)

Table 1 Fractionation coefficients for loss reactions with OH, O(¹D) and Cl and soil uptake. T denotes the temperature.

2.2 Inverse modeling with a variational approach

Inversions were performed using the Community Inversion Framework (CIF; Berchet et al., 2021). This framework was designed to rationalize and bridge development efforts made by the scientific community within the same flexible, transparent and open-source system. This system was recently enhanced by Thanwerdas et al. (2022a) to allow it to assimilate δ(¹³C,CH₄) together with CH₄ observations and to optimize both CH₄ emissions and the source signatures δ_source(¹³C,CH₄) at the same time. For the purpose of this study, δ(¹³C,CH₄) and δ(D,CH₄) observations are assimilated together in the same inversion, and the system optimizes both the source signatures (δ_source(¹³C,CH₄) and δ_source(D,CH₄)).

The notations introduced here to describe the variational inversion method follow the convention defined by Ide et al. (1997) and Rayner et al. (2019). x is the control vector and includes all the variables optimized by the inversion system. Prior information about the control variables is included in the vector x^b. Its associated errors are assumed to be unbiased and Gaussian and are described within the error covariance matrix B.

Here, the observation vector y^o includes all available observations, namely the atmospheric CH₄ amount fraction, δ(¹³C,CH₄) and δ(D,CH₄) data for one sensitivity test. The associated errors are also assumed to be unbiased and Gaussian and are described within the error covariance matrix R. This matrix accounts for all errors contributing to mismatches between simulated and observed values.

ℋ is the observation operator that projects the control vector x into the observation space. This operator mainly consists of the CTM but is also followed by spatial, time and isotope-conversion operators. Following Thanwerdas et al. (2022a), prescribed source signatures and CH₄ fluxes are first combined to generate isotope fluxes. These fluxes are then fed to the model to simulate the mixing ratios of the different isotopes over the time period considered. After the forward run, the simulated fields are interpolated to produce simulated equivalents of the observed amount fractions and isotopic compositions at specific locations and times, ensuring that a comparison between simulations and observations is possible. Adjoint versions of these forward operations are also implemented in order to perform the complementary adjoint run.

In a variational formulation of the inversion problem that allows ℋ to be nonlinear, the cost function J is defined as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}J\left(\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}^b\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}^b\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathcal{H}\left(\mathit{x}\right)-{\mathit{y}}^{\mathrm{o}}\right)}^T{\mathbf{R}}^{-\mathrm{1}}\left(\mathcal{H}\left(\mathit{x}\right)-{\mathit{y}}^{\mathrm{o}}\right).\end{array}\end{array}$$

Here, the minimum of J is reached iteratively with the descent algorithm M1QN3 (Gilbert and Lemaréchal, 1989), which requires several computations (40–50) of the gradient of J with respect to the control vector x:

$$\begin{array}{}\text{(4)}& \mathrm{\nabla }{J}_{\mathit{x}}={\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}^b\right)+{\mathcal{H}}^{*}\left({\mathbf{R}}^{-\mathrm{1}}\left(\mathcal{H}\left(\mathit{x}\right)-{\mathit{y}}^{\mathrm{o}}\right)\right).\end{array}$$

ℋ^* denotes the adjoint operator of ℋ.

The reference inversion (INV_REF) assimilates CH₄ and δ(¹³C,CH₄) observations over 1998–2018. CH₄ emissions and δ(¹³C,CH₄) source signatures for five categories of emissions – biofuel and biomass burning (BB), wetlands (WET), fossil fuels and geological sources (FFG), agriculture and waste (AGW), and other natural sources (NAT) – are optimized. The initial conditions for CH₄ and δ(¹³C,CH₄) are also optimized (see Sect. 2.6).

2.3 Prior emissions and uncertainties

For prior CH₄ emissions, we adopt the bottom-up estimates compiled for the inversions performed as part of the Global Methane Budget and described in detail in Saunois et al. (2020). In short, anthropogenic (including biofuel) and fire emissions are based on EDGARv4.3.2 (Emissions Database for Global Atmospheric Research version 4.3.2; Janssens-Maenhout et al., 2019) and GFED4s (Global Fire Emissions Database version 4s; van der Werf et al., 2017), respectively. Statistics from British Petroleum (BP) and the Food and Agriculture Organization of the United Nations (FAO) have been used to extend EDGARv4.3.2 (which ends in 2012) until 2017. The natural source emissions are based on averaged literature values: Poulter et al. (2017) for wetlands, Kirschke et al. (2013) for termites, and Lambert and Schmidt (1993) and Etiope (2015) for geological (onshore) sources and oceanic sources that include geological (offshore) and hydrate sources. Prior emissions for 2018 are set equal to the emissions in 2017. Globally averaged emissions over 1998–2018 are listed in Table 2.

BB emissions are the combination of biomass-burning emissions from GFED4s and biofuel-burning emissions from EDGARv4.3.2. FFG emissions are the combination of oil, gas, coal, industrial and transport emissions from EDGARv4.3.2 and geological (onshore) source emissions from Etiope (2015), whose global emissions were scaled down to 15.0 Tg  a⁻¹ in the protocol of Saunois et al. (2020). AGW emissions are the combination of enteric fermentation, rice agriculture, manure management and waste emissions from EDGARv4.3.2. NAT emissions are the combination of termite and oceanic emissions, i.e., emissions from all natural sources apart from wetlands and geological sources.

Table 2 Information about emissions and flux-weighted isotopic signatures for the different categories. Emissions and source signatures are averaged over 1998–2018. The uncertainty (unc.) indicates the prior uncertainty as a percentage of the square of the maximum of the prior emissions over the cell and its eight neighbors during each month (or over a continental region for the signatures). This uncertainty is used to fill the matrix B. The number of optimized scaling factors (optim.) can be either (1) 3PMPG: three scaling factors per month and per grid cell, (2) PYR: one scaling factor per year and per continental region, or (3) PR: one scaling factor per continental region for the full assimilation window.

Emissions are optimized at the grid-cell scale (one scaling factor per grid cell). For each category, diagonal elements of the matrix B are filled with the variances set to 100 % of the square of the maximum of the prior emissions over the cell and its eight neighbors during each month. Spatial error correlations (off-diagonal elements) are prescribed using an e-folding correlation length of 500 km on land and 1000 km over the oceans, without any correlation between land and ocean grid points. No temporal error correlations are prescribed.

2.4 Prior source signatures and uncertainties

δ_source(¹³C,CH₄) values for each emission category are also optimized and therefore included in the control vector. Prior information is built using the references given in Table 3. When regional information could be found, regional source signature values were prescribed for 11 continental regions (see Fig. 1, lower-right panel) for each subcategory. Isotopic signatures for subcategories are flux-weighted averaged to create signatures for categories, which results in signatures that are grid-cell dependent, although signatures for subcategories are set constant within a continental region. Signatures for wetlands are the only ones prescribed at the grid-cell scale, following Ganesan et al. (2018). We optimize the source signatures at the regional scale rather than at the grid-cell scale to avoid substantial posterior differences between two adjacent grid cells. At present, there would not be enough data to corroborate, explain, or reject such differences. We therefore apply only one scaling factor per continental region and per category.

Livestock source signatures have likely been decreasing over time since the 1990s due to changes in the C₃ $/$ C₄ diet within the major livestock producing countries (Chang et al., 2019). Also, FFG regional source signatures can vary over time due to variations in the contributions from different sectors (coal, oil and gas) to the emissions of a specific region (Schwietzke et al., 2016; Feinberg et al., 2018). For AGW and FFG source signatures, we therefore optimize one scaling factor per year for each continental region. As for the other emission categories, only one scaling factor for the entire period and for each continental region is optimized. Error covariances are prescribed following the same methods as those applied to CH₄ emissions.

As for δ_source(D,CH₄), we adopted global values suggested by Warwick et al. (2016) and in agreement with the intervals given by Röckmann et al. (2016). One exception is for WET sources, for which the boreal (−360 ‰) and tropical (−320 ‰) regions are differentiated. All values are summarized in Table 2. For each category and each continental region, only one scaling factor is optimized for the entire period.

Table 3 Global flux-weighted values and references for δ_source(¹³C,CH₄) source signatures associated with the different emission categories and subcategories. Values for subcategories are taken from the literature and prescribed either globally (G), regionally (R, see Fig. 1, lower-right panel), or at the pixel scale (P). E19: Etiope et al. (2019); CH19: Chang et al. (2019); GA18: Ganesan et al. (2018); TH18: Thompson et al. (2018); SH17: Sherwood et al. (2017); SH16: Schwietzke et al. (2016); WA16: Warwick et al. (2016); ZA16: Zazzeri et al. (2016); TO12: Townsend-Small et al. (2012); KL10: Klevenhusen et al. (2010); BO06: Bousquet et al. (2006); BR01: Bréas et al. (2001); SA01: Sansone et al. (2001); CH00: Chanton et al. (2000); HO00: Holmes et al. (2000); CH99: Chanton et al. (1999); BE98: Bergamaschi et al. (1998); LE93: Levin et al. (1993).

Figure 1 Prior estimates of δ_source(¹³C,CH₄) isotopic signatures for each of the five emission categories averaged over the 1998–2018 period. The regions over which the values are optimized are shown in the lower-right panel. WET source signatures are dependent on the latitude, with more depleted values occurring in boreal regions than in tropical regions. BB source signatures are dependent on the vegetation (C₃ $/$ C₄). Burning C₄ vegetation in tropical regions releases CH₄ that is more ¹³C enriched than the CH₄ released when burning C₃ vegetation. The AGW source signature is dependent on the country/region and the C₃ versus C₄ livestock diet. FFG source signatures mainly depend on both the location and the contributions from coal, oil, gas and geological sources to the total FFG emissions of a specific country/region. For example, China's large ¹³C-enriched coal emissions greatly contribute to the FFG source signature in this region, which is notably ¹³C enriched compared to other regions.

δ_source(¹³C,CH₄) uncertainty values that are used to fill the diagonal elements of the matrix B are summarized in Table 2. These values have been chosen by compiling data from several studies (Sherwood et al., 2017; Ganesan et al., 2018; Feinberg et al., 2018; Zazzeri et al., 2016). The observed variability over the whole globe (standard deviation σ or minimum–maximum range) presented in these studies for each category is compiled and applied here as an uncertainty (1σ), thus adopting the same value for all regions. Note that for BB sources, Sherwood et al. (2017) indicates a global standard deviation of about 20 %. However, this value is not weighted by the proportion of C₃ versus C₄ vegetation. Therefore, we inflated this uncertainty to 30 % to account for the uncertainty in the type of vegetation. δ_source(D,CH₄) uncertainty values were derived from the minimum–maximum ranges suggested by Röckmann et al. (2016). We could have also used the standard deviation provided by Sherwood et al. (2017). However, as the amount of data for AGW, BB, WET and NAT source signatures is very low compared to the δ_source(¹³C,CH₄) values, we prefer to use larger uncertainties and examine whether the assimilation of δ(D,CH₄) observations modifies the results. For future studies, additional δ(D,CH₄) data to derive realistic regional estimates, especially for non-fossil sources, would be invaluable.

We acknowledge the fact that our methods are not perfect and that the prescribed regional uncertainties might be too large compared to the regional observed uncertainties that are currently estimated, especially for δ_source(¹³C,CH₄). As our inversion system is used for the first time over a time period exceeding 10 years, it is difficult to predict the influence of the setup on the results. As the time required to run an inversion is very high at the moment (see Sect. 2.9), we prefer to assess the behavior of this system in response to a simple (and probably slightly loose) setup and to estimate whether such uncertainties are small enough to help better constrain the CH₄ budget. Additionally, source signature data representativeness is generally poor, owing in part to the small number of samples and the lack of data for several regions, particularly for non-fossil sources. It might therefore be challenging to derive a robust uncertainty using a data-driven approach for each region of the world. However, there is definitely room for improvement, and future work will build on the present work to improve this methodology and assess and prescribe better uncertainties.

Note that random uncertainties are only one side of the coin when it comes to uncertainties. Systematic uncertainties must also be investigated when using isotopic constraints. In particular, Oh et al. (2022) derived source signature maps for wetland sources that carry systematic uncertainties. Typically, inverse modelers address such uncertainties by conducting numerous inversions using parameters designed to account for these systematic errors. The high computational cost associated with our system prevents us from running a large number of inversions. Nevertheless, only one scaling factor is applied for each region and each category. Consequently, there is a strong regional correlation between random errors. To some degree, this approach enables the detection and correction of regional systematic errors by our system. It is important to note, however, that this correction does not rely on any existing sensitivity analysis (e.g., Lan et al., 2021; Basu et al., 2022), but rather uses solely the information provided by atmospheric isotopic observations.

2.5 Observations

2.5.1 Assimilated data

Our study uses observations from the NOAA Global Monitoring Laboratory (NOAA GML) Global Greenhouse Gas Reference Network. CH₄ measurements are made by NOAA GML (Lan et al., 2022), and isotopic measurements are made at the Stable Isotope Laboratory at the Institute of Arctic and Alpine Research (INSTAAR) (White et al., 2021, 2016). This ensemble was selected to provide the largest number of consistent CH₄ and isotopic data since 1998. Seventy-nine stations (four of which were mobile stations) provided CH₄ measurements between 1998 and 2018 (not necessarily over the full period), 22 stations provided δ(¹³C,CH₄) measurements between 1998 and 2018, and 15 stations provided δ(D,CH₄) measurements between 2005 and 2010 (see Fig. 2). Missing CH₄ instrumental errors are filled with the maximum value of this error at the station over the monitoring period. For δ(¹³C,CH₄) and δ(D,CH₄) measurements, missing instrumental errors are filled with a value of 0.1 ‰ and 3 ‰, respectively (Quay et al., 1999). Variances (diagonal elements) in the covariance matrix R are defined as the sum of the instrumental and model errors (variances). For each station and each year, we used the residual standard deviation (RSD) between the measurements and a fitting curve function as a proxy for the model error (Thanwerdas et al., 2022a; Locatelli et al., 2015; Bousquet et al., 2006). The fitting function includes three polynomial parameters (quadratic) and eight harmonic parameters, sine and cosine, as in Masarie and Tans (1995). We also remove outliers that are outside three times the residual standard deviation, as such extreme values cannot be reasonably reproduced at the horizontal grid resolution of LMDz. Typical values for observation errors are 20 nmol mol⁻¹ for CH₄, 0.3 ‰ for δ(¹³C,CH₄) and 7 ‰ for δ(D,CH₄).

2.5.2 Satellite data used for comparison

The Greenhouse Gases Observing Satellite (GOSAT), which carries a Fourier-transform spectrometer within its Thermal And Near-infrared Sensor for carbon Observation (TANSO-FTS) (Kuze et al., 2016), provides radiance measurements in a spectral band centered on a value close to 1.6 µm, in which CH₄ has a high absorption capacity. The University of Leicester's retrieval algorithm is able to produce column-average dry air amount fractions of CH₄ from these radiances. Although this quantity is commonly referred to with the symbol XCH₄ in the existing literature, it is denoted here by $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ because this is considered a more valid notation. We use version 9.0 of the GOSAT Proxy $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ dataset provided by the University of Leicester (Parker et al., 2020) to evaluate the $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ after the inversion process. To this end, vertical profiles simulated by LMDz-SACS are sampled at the observation location and time and convolved with the retrieval of the prior vertical profiles and column averaging kernels provided by the University of Leicester. Finally, within each grid cell, all the individual $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ differences between the satellite observations and model outputs are averaged.

Satellite data are not assimilated here because inversions assimilating both satellite and surface data have not yet been performed with LMDz-SACS. Before using satellite data and isotope data together in an inversion, we need to rigorously assess the added value of satellite data without isotope constraints.

Figure 2 Locations of surface CH₄, δ(¹³C,CH₄) and δ(D,CH₄) stations in the NOAA GML network. Samples from several stations are retrieved and analyzed by INSTAAR to provide δ(¹³C,CH₄) and δ(D,CH₄) observations. More information about the stations can be found in Appendix A. Note that mobile stations (AOC, PAO, POC and WPC) are each indicated by a single point for clarity.

2.6 Initial conditions

To infer the initial CH₄ and δ(¹³C,CH₄) conditions for 1998, we run an inversion between 1988 and 1998 using the same prior emissions and isotopic signatures as that of INV_REF. We assimilate CH₄ measurements from the NOAA GML network (56 stations) and δ(¹³C,CH₄) measurements retrieved at five stations across the globe by the University of Washington (UW) between 1988 and 1996 (Quay et al., 1999; Bousquet et al., 2006), which we offset by 0.1 ‰ to account for measurement differences between INSTAAR and UW (Umezawa et al., 2018).

We also run a forward simulation from 1998–2010 to obtain a good spatial distribution of the δ(D,CH₄) field and then apply a global offset to match the observed mean δ(D,CH₄) value between 2005 and 2010. As we acknowledge that both methods are not perfect, considering the equilibration times of these isotopic compositions (Tans, 1997), we also prescribe large uncertainties in these initial conditions: 10 % for CH₄, 3 % for δ(¹³C,CH₄) and 20 % for δ(D,CH₄). To optimize the initial conditions, the globe is regularly discretized using latitudinal and longitudinal bands. A step of 30^∘ is applied to generate the bands, resulting in 6 × 12 = 72 regions. One scaling factor is optimized for each of these regions.

2.7 Description of the sensitivity tests

The reference inversion was first introduced in Sect. 2.2, and its setup is detailed in previous sections. Three sensitivity tests were conducted to investigate the influence of the setup of our system on posterior estimates when assimilating isotopic observations:

• INV_CH4 is an inversion that only assimilates CH₄ observations, not δ(¹³C,CH₄) nor δ(D,CH₄) observations, and thus only optimizes CH₄ emissions for the five categories.

• INV_DD assimilates CH₄, δ(¹³C,CH₄) and δ(D,CH₄) observations and optimizes both δ_source(¹³C,CH₄) and δ_source(D,CH₄) source signatures. Note that δ(D,CH₄) observations only span the period from 2005–2010, and therefore the full run cannot be fully constrained by this data.

• INV_LOCKED is an inversion that assimilates CH₄ and δ(¹³C,CH₄) observations but does not optimize δ(¹³C,CH₄) source signatures (fixed to prior values). This run considers source signatures fixed to prior values and thus investigates the influence of overconstrained isotopic signatures on posterior estimates.

We also investigate the influence of the OH interannual variability (IAV) on our results. The OH IAV in the troposphere is usually derived from inversions using CH₃CCl₃ constraints (Montzka et al., 2011; Rigby et al., 2017; Turner et al., 2017; Naus et al., 2019) or using global atmospheric chemistry–climate models (He et al., 2020; Dalsøren et al., 2016). CH₃CCl₃ inversion-based studies suggest a decrease in post-2005 OH after a peak in 2000–2002. By contrast, the chemistry modeling studies derive a post-2005 stabilization after a quasi-continuous increase between 1990 and 2005, consistent with the OH IAV estimated by LMDz-INCA (see Fig. 3). We therefore perform two more sensitivity tests:

• INV_TURNER is designed to investigate the influence of the IAV on our results. We apply the IAV suggested by Turner et al. (2017) to the OH-INCA field. The associated OH field is named OH-TURNER, and its global concentrations decrease by 7 % between 2002 and 2014.

• INV_FLATOH removes the IAV from our OH-INCA field by prescribing the concentrations for the year 2000 over the full period. The associated field is named OH-FLAT.

All the sensitivity tests are summarized in Table 4.

Figure 3 Time series of the global volume-weighted tropospheric OH annual concentrations for the 1998–2018 period. OH-INCA was simulated by the LMDz-INCA chemistry model. OH-TURNER was obtained by applying the IAV from Turner et al. (2017) to the OH-INCA field. OH-FLAT has no interannual variability and concentrations are set equal to those of OH-INCA in 2000. In 1980, the OH-INCA mean concentration is very close to 10.0 × 10⁵ cm⁻³, which was also taken as a reference value by Turner et al. (2017). This year is therefore taken as a reference to derive the anomalies.

Table 4 Description of the sensitivity tests.

2.8 Analysis period

Thanwerdas et al. (2022a) suggested that the results of the inversion should be discarded for the period up to 2–3 years after the beginning and the period 2–3 years before the end of the assimilation window. Although the term “spin-up” is not quite appropriate for the beginning of the window because the reason for discarding those results is slightly different, the outcome is still similar. The spin-up time, for an inversion, typically refers to a period at the beginning of the inversion when the errors in the assumed initial concentration field might influence the posterior fluxes. If the adopted spin-up period is too short, the inversion may fit the data by compensating for errors in the initial conditions with artificial emission adjustments. If the initial concentrations are also optimized, which is the case here, this effect can be reduced (Houweling et al., 2017). However, source signatures are optimized, and, for a certain period after the start of the inversion, it is easier for the system to optimize the initial δ(¹³C,CH₄) fields than the source signatures to fit the δ(¹³C,CH₄) data. Thanwerdas et al. (2022a) found that the optimized source signatures slowly move away from the prior value over time. After 2–3 years, the posterior value finally reaches a new and rather stable state. In other words, as the influence of the initial conditions on the isotopic composition decreases, the system prefers to optimize the source signatures, hence slowly reaching the posterior value. The equilibration time for δ(¹³C,CH₄) and δ(D,CH₄) is larger than the CH₄ equilibration time (Tans, 1997), and therefore the signatures are more affected than the fluxes. For the end of the inversion, the reason for discarding the results is a lack of constraints, resulting in a slow return to the prior value. In this case, the use of the term “spin-down” is correct.

In addition, the strong 1997–1998 El Niño event leads to fire emission anomalies of about 20 Tg  a⁻¹ according to GFED4s data and studies (Bousquet et al., 2006; Langenfelds et al., 2002). Similarly, the following 1999–2000 La Niña event produced a wetland emission anomaly that persisted until 2002 (Zhang et al., 2018). Finally, OH concentrations were also likely affected by this El Niño–Southern Oscillation (ENSO) phase (Zhao et al., 2020a). We therefore only analyze the results of the 2002–2014 period to limit the consequences of these effects. This period of time is large enough to explain the variations that caused the post-2007 CH₄ renewed growth and the associated δ(¹³C,CH₄) shift to more negative values.

2.9 Computational aspects

The same convergence criterion was used for all inversions in order to ensure consistency between the results. The minimization process was stopped when at least 35 iterations (forward + adjoint runs) had been performed and the gradient norm ratio had fallen below 1 % of its initial value for four successive iterations.

A similar number of iterations (approx. 40) were necessary for all sensitivity tests. About 260 CPU hours were necessary to run a single iteration on LSCE (Laboratoire des Sciences du Climat et de l'Environnement) computational clusters consisting of Intel^® Xeon^® Gold 5317 central processing units (CPUs) with a frequency of 3.00 GHz. For this work, eight CPUs were run in parallel, resulting in a runtime of 32.5 h for a single iteration. More CPUs could not increase the overall performance because of some I/O (input/output) limitations of our offline model. With only one tracer to simulate, INV_CH4 therefore required about 2 months to reach the convergence criterion. Because the runtime is proportional to the number of simulated tracers, double this runtime was needed for the other inversions (two tracers), except for INV_DD, which required four times this runtime (four tracers).

While the number of CPU hours needed for these complex inversions remains reasonable, the overall runtime is excessive. It is therefore an important limitation of our system. Further developments of parallelization methods are being implemented to enable a significant reduction of the computational cost (e.g., Chevallier, 2013). This method consists of breaking down the full assimilation window into multiple sub-windows and running smaller inversions in parallel for each sub-window. If source signatures remain constant, we expect the results to closely resemble those of a longer-term window inversion. Conversely, if source signatures are optimized, the influence of the initial conditions on the atmospheric isotopic composition might persist over a time that is larger than the length of the sub-window. In this case, source signatures might remain unchanged and the results could be impacted. Therefore, it is crucial to rigorously validate this parallelization method before interpreting its outcomes.

3 Results

In this section, we first verify the quality of the model's fit to the constraining observations and evaluate it against independent data (Sects. 3.1 and 3.2). After this, we examine the posterior estimates of our reference inversion for emissions and source signatures and compare them to prior estimates (Sects. 3.3 and 3.4). Subsequently, we attribute the CH₄ increase and downward shift in δ(¹³C,CH₄) post-2007 to changes in CH₄ emissions and source signatures (Sects. 3.5 and 3.6). Finally, we analyze the sensitivity of our results to setup modifications (Sects. 3.7 and 3.8).

3.1 Model–observation agreement

Before analyzing the optimized emissions and source signatures, we verify the quality of the model's fit to the constraining observations and evaluate it against independent data. A good fit shows that the system is operational over long time periods and that posterior emissions and source signatures are consistent with the observed state of the atmosphere.

The observed globally averaged CH₄ amount fraction as well as the observed globally averaged δ(¹³C,CH₄) isotopic composition at the surface are well captured by all the posterior simulations (see Fig. 4b, d). INV_REF shows a Pearson's moment correlation coefficient r of 0.994 for CH₄ (RMSE is 2.8 nmol mol⁻¹) and 0.936 for δ(¹³C,CH₄) (RMSE is 0.04 ‰). The posterior simulation therefore captures the observations much better than the prior simulation does (RMSEs of 71.2 nmol mol⁻¹ and 1.44 ‰, respectively). The inversion that best captures the δ(¹³C,CH₄) isotopic composition is INV_DD, with a RMSE of 0.02 ‰. This shows that assimilating δ(D,CH₄) observations slightly increases the agreement with isotopic observations without performing additional iterations.

The 2002–2007 δ(¹³C,CH₄) stabilization is well reproduced by the model in INV_REF, with a mean RMSE of 0.02 ‰ over the period. However, the post-2007 trend is not as consistent (0.05 ‰), mainly due to an overestimation of the decreasing rate (0.03 ‰  a⁻¹ against 0.02 ‰  a⁻¹). The simulated δ(¹³C,CH₄) seasonal cycle amplitude is also slightly smaller than the observed one (0.12 ‰ against 0.14 ‰, respectively), although the two signals are well phased. Note that our results are, however, still within the prescribed observation uncertainty range.

For the sake of completeness, we also provide a comparison between δ(D,CH₄) observations and prior and posterior simulations from INV_DD in Appendix B (Fig. B1). After the inversion, simulations capture the observed δ(D,CH₄) data much better, reducing the RMSE from 9.3 ‰ to 1.2 ‰. Although the δ(D,CH₄) data are much more limited than the δ(¹³C,CH₄) data, linear regressions indicate a small negative trend (−0.23 ± 0.12 ‰  a⁻¹) between 2005 and 2009. Additionally, the trend is positive between 2005 and 2007 (+0.86 ±  0.42 ‰  a⁻¹) and negative between 2007 and 2009 (−0.47 ± 0.25 ‰  a⁻¹). This shows that δ(D,CH₄) observations might also carry some information about the renewed growth of CH₄ post-2007.

Figure 4 Posterior agreement between INV_REF and assimilated observations. Panels (a) and (c) show the posterior CH₄ and δ(¹³C,CH₄) biases at the surface stations over the 2002–2014 period. Marine boundary layer stations are indicated by squares rather than circles. Panels (b) and (d) show the observed (solid black line) and simulated (solid red line) globally averaged trends in CH₄ and δ(¹³C,CH₄). The red-shaded area shows the minimum and maximum values over the sensitivity tests. The gray-shaded area shows the standard error of the globally averaged observed trend. This error is based on the error prescribed in the matrix R, i.e., the sum of the measurement and model errors. The same figure with prior data is provided in Appendix B (Fig. B2).

The model–observation agreement varies across the stations for both CH₄ and δ(¹³C,CH₄) (see Fig. 4a and c). CH₄ at marine boundary layer (MBL) stations (i.e., where site samples consist mainly of well-mixed MBL air) is very well reproduced by the model (mean RMSE of 17.7 nmol mol⁻¹ and mean bias of −1.87 nmol mol⁻¹). The model has more difficulties in simulating amount fractions at several polluted stations, such as Lac La Biche, Canada (54.95^∘ N, 112.45^∘ W), Shangdianzi, People's Republic of China (40.65^∘ N, 117.12^∘ E), Anmyeon-do, Republic of Korea (36.54^∘ N, 126.33^∘ E), or the Southern Great Plains, United States (36.62^∘ N, 97.48^∘ W), presumably owing to transport errors, representation errors and/or inaccurate estimates of CH₄ prior fluxes around these stations. The δ(¹³C,CH₄) at MBL stations is generally correctly simulated, with RMSEs of 0.2 ‰–0.3 ‰, comparable to the prescribed uncertainties. However, posterior simulations slightly overestimate δ(¹³C,CH₄) in Northern America and underestimate it in Central America and temperate North America. This suggests an over- and underestimation of flux-weighted source signatures in these regions, respectively. This is further investigated in Sect. 3.4.

Figure 5 Mean posterior model–observation differences in $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ averaged over 2010 and gridded at the model resolution. Model outputs were obtained using the posterior estimates of INV_REF and by applying the averaging kernels provided by the University of Leicester. Model–observation differences at assimilated surface stations are also displayed for comparison. We acknowledge that the temporal sampling by surface stations is not identical to that of satellite data, which might affect the comparison. To limit this effect, only stations providing at least one observation for each month of 2010 are displayed to reduce the seasonal influence. MBL stations are indicated by squares rather than circles.

3.2 Comparison of model-optimized with satellite-derived column-average CH₄ amount fractions

For comparison, we performed one forward simulation with posterior fluxes obtained with INV_REF to compare our simulated $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ to independent (i.e., not assimilated) satellite observations in 2010 and evaluate the optimized atmosphere (see Fig. 5). The posterior mean bias is −13.0 nmol mol⁻¹, indicating that the GOSAT observations are higher overall than our optimized $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$, even after the inversion. Further analysis reveals that the mean bias (−17.5 nmol mol⁻¹) in the tropics (30^∘ S–30^∘ N) is larger in absolute value than the bias in the northern mid-latitudes (30–60^∘ N; −5.7 nmol mol⁻¹) or in the northern high latitudes (60–90^∘ N; +3.0 nmol mol⁻¹). Parker et al. (2020) also reported observing a negative $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ mean bias (−6.55 nmol mol⁻¹) when using the TM5 model and posterior estimates deduced from a surface-based inversion. Similar to ours, the TM5 biases were mostly located in the tropics and northern mid-latitudes. Ostler et al. (2016) reported that model errors when simulating stratospheric CH₄ amount fractions could contribute to the $\stackrel{\mathrm{‾}}{X}\left({\mathrm{CH}}_{\mathrm{4}}\right)$ bias. However, they did not find a strong improvement for LMDz and TM5 when replacing model simulations with MIPAS (Michelson Interferometer for Passive Atmospheric Sounding) stratospheric CH₄. While this suggests that the biases in the GOSAT simulation presented here are probably not caused by stratospheric discrepancies, it is important to note that the same authors conclude that current satellite measurements of stratospheric CH₄ may lack the precision necessary to eliminate these biases.

If we assume that the biases are solely the result of emission discrepancies, our findings indicate that the estimated tropical posterior emissions from our inversions might still be underestimated. Although posterior biases are lower at the surface stations in the tropics (the mean value is −5.3 nmol mol⁻¹), the number of stations is limited in this area, especially in South America and Central Africa. Saunois et al. (2020), using configurations similar to ours but no isotopic constraints, found that differences between emissions from GOSAT-based and surface-based inversions mainly occurred in the tropical regions.

The in situ-only and the GOSAT-only inversions performed by Lu et al. (2021) provided 113 and 212 independent pieces of information, respectively, highlighting that additional constraints can be gained using satellite data. As tropical fluxes likely had a significant influence on the renewed increase of CH₄ around 2007, it would be interesting to assimilate satellite observations to increase the constraints in the tropics. Furthermore, jointly assimilating satellite observations and isotopic observations might be valuable because tropical emissions largely dominate the total global release of CH₄. Consequently, a change in tropical emissions might influence the global flux-weighted source signature and impact the results of an inversion performed with isotopic constraints. Our system is capable of performing such a joint assimilation; however, since we are analyzing the influence of adding isotope constraints here, we prefer to assimilate only surface data as a first step.

3.3 Posterior–prior emission differences

Global emissions are estimated by INV_REF at 590 Tg  a⁻¹ when averaged over the 2002–2014 period, which is larger than prior estimates by 28 Tg  a⁻¹ (see Fig. 7). This change mainly arises from increases in Asia (+15 Tg  a⁻¹), Central and South America (+6 Tg  a⁻¹), and Africa (+3 Tg  a⁻¹). About 50 % and 25 % of the increase in Asia is due to the AGW and FFG categories, respectively, suggesting that the prior estimates in EDGARv4.3.2 for these regions are underestimates. However, global emissions estimated by inverse modeling are strongly dependent on the chemical loss prescribed in the CTM. OH is responsible for most of this loss and therefore the prescribed OH field greatly influences the results of the inversion. In particular, Zhao et al. (2020b) showed that a 1 × 10⁵ cm⁻³ increase in prescribed OH concentrations leads to an increase in CH₄ global posterior emissions of 40 Tg  a⁻¹. Here, we are more interested in the trends in various CH₄ emission categories and their contributions to the total emissions. Therefore, we have only used a single OH field with several trends. The influences of the trends on the dedicated inversions are further discussed in Sect. 3.8.

The posterior global distribution of emissions across the individual categories is only slightly different from the prior one (see Fig. 6). Relative posterior–prior emission differences averaged over 2002–2014 are larger for WET (+7 %) than for AGW (+4 %), FFG (+6 %) or BB (+5 %). The increase in WET emissions is mainly located in the Amazon basin (43 %) and is responsible for a small shift in the WET contribution to the total emissions (from 32.1 % to 32.6 %), which is offset by a similar reduction in the AGW contribution (from 37.9 % to 37.4 %). The tropics (90^∘ S–30^∘ N), the northern mid-latitudes (30–60^∘ N) and the high latitudes contribute about 60 %, 35 % and 5 % to the global emissions, respectively. Apart from a small reduction in the contribution from high latitudes, the latitudinal distributions of emissions are not modified.

Figure 6 Prior (left) and posterior (right) contributions from three latitudinal regions to global CH₄ emissions for each emission category. Posterior emissions are taken from INV_REF. The latitudinal bands are (1) the tropics (90^∘ S–30^∘ N), (2) the northern mid-latitudes (30–60^∘ N) and (3) the high latitudes (60–90^∘ N). Emissions and source signatures for each region and category are given in the associated bars. The total emissions and global source signatures for each category are given on the top of each bar.

3.4 Posterior–prior source signature differences

Global and regional source isotopic signatures are calculated using a flux-weighted average to account for the global and regional source mixture, respectively. Therefore, they may be modified by the system due to a source mixture change and/or a source signature change in a specific region.

The inversion system shifts the global source signature considerably upward from −54.5 ‰ to −52.5 ‰ (see Fig. 6). The global signature is highly constrained by the fractionation coefficients and the concentrations of radicals (OH, Cl and O¹D) prescribed in the CTM. Our posterior global source signature is higher (less negative) compared to other estimates (Sherwood et al., 2017; Rigby et al., 2017; Schaefer et al., 2016), mainly because we chose to prescribe Cl concentrations and an OH fractionation that are at the low ends of the existing ranges. Each additional percent of oxidation caused by the prescribed Cl sink would lead to a drop in the global source signature of about 0.5 ‰ (Thanwerdas et al., 2022b; Strode et al., 2020). Using other recent estimates of Cl tropospheric concentrations (see Sect. 2.1), our posterior global source signature would range between −54.5 ‰ and −52.5 ‰. In addition, the application of the fractionation value derived by Cantrell et al. (1990) instead of that derived by Saueressig et al. (2001) would likely shift the global signature downward by another 1.5 ‰.

AGW, BB and FFG source signatures are shifted upward by 0.7 ‰, 0.7 ‰ and 1.3 ‰, respectively. Most notably, the posterior global flux-weighted WET signature is considerably higher (−56.6 ‰) than the prior estimate (−60.8 ‰) (see Fig. 6) with regard to recent estimates (Sherwood et al., 2017; Feinberg et al., 2018; Ganesan et al., 2018; Oh et al., 2022). This global shift mainly arises from upward regional source signature shifts in the tropics (+3.2 ‰) and in the northern mid-latitudes (+7.1 ‰), which together contribute 97 % of the posterior–prior global WET isotopic signature difference. The remaining contribution is due to an increase in the contribution from tropical WET emissions. Our posterior global estimate of the WET source signature strongly disagrees with the recent estimates. In Appendix B, Fig. B3 compares our prior and posterior signatures to observations from Supplement Data 1 provided by Oh et al. (2022). Overall, prior estimates show a better agreement with observations than posterior estimates. This poor agreement suggests that prescribed uncertainties might be too large (at least for WET) and supports the idea that the system yields an important adjustment of the WET signature in order to keep fossil/microbial flux partitioning unchanged. Notably, the isotopic signature in Canada is shifted upward from −70.0 ‰ (prior) to −59.4 ‰ (posterior), whereas that of Russia is shifted downward from −68.7 ‰ to −73.8 ‰. Although this demonstrates that the system is capable of applying offsets with different signs across different regions, it also appears to be unphysical, as the processes driving the source signatures in these high-latitude regions are similar (Ganesan et al., 2018; Oh et al., 2022). Nevertheless, the number of observations of WET source signatures is small and local uncertainties are considerable, especially in the tropics, where emissions from WET are the largest and where our system applies the most impactful adjustment. It is therefore difficult to invalidate the posterior adjustment. Further investigation, including a better assessment and prescription of the random and systematic uncertainties, is needed.

Figures B4 and B5 in Appendix B show the full temporal variations for the prior and posterior source signatures. For FF, high variations indicate a change in activities associated with fossil fuel extraction, e.g., switching from one location with a specific signature to another, transitioning from one fuel type (oil, gas, coal) to another, or a combination of both. For example, the substantial shift that occurred around 2009 in the United States was caused by a large increase in emissions from the extraction of natural gas. As we chose not to prescribe temporal error correlations between different years, the system is free to optimize each year independently to better fit δ(¹³C,CH₄) observations. For certain continental regions, such as Africa, temperate Asia, or South Asia, the interannual variability of the source signature adjustments is large and rather unrealistic, especially when compared to the emission adjustments in the same regions (see Fig. 7). It is unlikely that these changes occurred without detectable changes in emissions in the same areas, especially for temperate Asia, which exhibits larger emissions from FF than in the United States. These results suggest the need to prescribe yearly temporal error correlations to dampen this artificial interannual variability. However, the example from the United States also indicates that large changes can occur, and it is reasonable to assume, considering the lack of isotopic data, that the prior data might not contain any information about potential substantial changes. Therefore, while implementing stronger temporal correlations could be a way to mitigate unrealistic interannual variability for this category, it diminishes the likelihood of detecting such changes that remain undetected by the prior data. Nevertheless, it might be sufficient to reduce the prescribed uncertainties in the source signatures in order to balance out the pressure applied by the system on the emissions and the source signatures. Overall, the same reasoning applies to AGW, although there is no evidence from the prior data that AGW source signatures can change as rapidly as those of FF. Due to the scarcity of existing data on the temporal variability of source signatures, designing a data-driven methodology to estimate potential temporal correlations, especially at the regional scale, remains highly challenging. Investigating the correlations that the system creates between the uncertainties associated with source signatures and fluxes could offer a promising avenue for extending the analysis. Due to the high computational cost of an inversion performed with our inversion system, it is impossible to derive robust posterior uncertainties. This impossibility is a major drawback, and additional studies with this system cannot be performed in the future without tackling this issue.

Figure 7 Map of changes in posterior emissions from INV_REF between 2002–2007 and 2007–2014 (center panel) and time series of posterior emission estimates from INV_REF for multiple regions and all categories (panels around the map). For each panel, the time series show anomalies around the 2002–2014 mean value. For each category in each panel, the associated mean value is displayed in the same color as the solid line. The units of the variations and means are Tg  a⁻¹. The design of the figure is inspired by Fig. 1 in Chandra et al. (2021).

3.5 Attribution of the post-2007 CH₄ increase

We now present trend results comparing the time periods before (2002–2007) and after (2007–2014) the renewed increase. Posterior global emissions show a net increase of 24.0 Tg  a⁻¹ (see Fig. 7). This occurred in most of the regions aside from Europe (−1.8 Tg  a⁻¹), Canada (−0.7 Tg  a⁻¹) and Oceania (−0.2 Tg  a⁻¹). China, South Asia (mainly India), temperate Asia, Southeast Asia and Africa accounted for 40 %, 18 %, 18 %, 10 % and 9 % of the associated positive increase (+26.7 Tg  a⁻¹), respectively. These results are consistent with prior information that estimated a rise of 27.5 Tg  a⁻¹. The large contribution from China to the global increase since 2002 agrees well with the regional estimate (40 %) from Thompson et al. (2015). We also estimate that Central and South America did not contribute to the renewed growth, in contrast with Chandra et al. (2021), who suggest a large contribution from Brazil (11.5 %) to the increase in global emissions. In this region, we find that small increases in AGW emissions (+1.4 Tg  a⁻¹) and FFG emissions (+0.7 Tg  a⁻¹) are offset by decreases in WET emissions (−0.8 Tg  a⁻¹) and BB emissions (−0.8 Tg  a⁻¹).

Global AGW emissions increased by 14.2 Tg  a⁻¹ between 2002–2007 and 2007–2014. Europe is the only region where these emissions substantially decreased (−1.7 Tg  a⁻¹). 80 % of the net AGW increase occurred in Asia, and the rest occurred in Africa and South America. FFG emissions increased by 14.9 Tg  a⁻¹, with 50 % of this increase occurring in China. These emissions notably decreased in Africa (−0.7 Tg  a⁻¹), Europe (−0.2 Tg  a⁻¹) and Canada (−0.1 Tg  a⁻¹). By contrast, WET emissions decreased by 2.4 Tg  a⁻¹, with 71 % of the net decrease located in Central and South America (33 %), Canada (25 %), and Africa (13 %). BB emissions also decreased by 2.7 Tg  a⁻¹, mainly due to decreases in Southeast Asia (55 % of the net decrease) and South America (24 %). Note that the analysis period does not include the 2015 El Niño event.

Our results therefore suggest that the renewed CH₄ growth post-2007 (until 2014) was equally and mainly driven by increases in global AGW (49 %) and FFG (51 %) emissions. The decreases in global WET and BB emissions as well as the increase in global OH concentrations partially balanced this renewed growth. These findings are in partial agreement with recent studies (Chandra et al., 2021; Jackson et al., 2020; Thompson et al., 2018; Saunois et al., 2017), although only Jackson et al. (2020) explained the renewed growth with equal contributions from the AGW and FFG categories. The small decrease in BB emissions is consistent with other estimates (Thompson et al., 2018; Worden et al., 2017), but the decrease in WET emissions does not agree with recent findings that suggest either a constant trend or a positive trend (Chandra et al., 2021; Zhang et al., 2018; McNorton et al., 2018; Poulter et al., 2017; Bader et al., 2017).

However, posterior global WET emissions show negative anomalies between 1998 and 1999 and positive anomalies between 1999 and 2004. These anomalies are mainly located in South America, where about 30 % of WET emissions originate. Zhang et al. (2018) suggested that the 1998–2000 ENSO (El Niño–Southern Oscillation) caused negative anomalies in WET emissions between 1998 and 2000 because of El Niño and subsequent positive anomalies between 2000 and 2002 because of La Niña. The fact that positive anomalies persist until 2004 rather than 2002 in our posterior emissions cannot be easily explained. Also, the positive anomalies last 4–5 years in total, which is not consistent with the 2–3 years inferred by Zhang et al. (2018). As AGW emissions are also large in South America, the inversion system might be wrongly attributing large but decreasing emissions between 2002 and 2004 to WET emissions rather than AGW emissions. If the period 2002–2004, which exhibits large positive anomalies, is discarded, we find a small increase of 0.3 Tg  a⁻¹ in global WET emissions between 2004–2007 and 2007–2014, which is more consistent with the studies mentioned before.

3.6 Attribution of the post-2007 downward shift in δ(¹³C,CH₄)

The posterior global flux-weighted δ_source(¹³C,CH₄) decreased from −52.1 ‰ in 2002 to −53.1 ‰ in 2010 and then experienced an upturn to −52.5 ‰ in 2012–2014. Notably, between 2007 and 2010, the decline in the source signature was rapid (−0.2 ‰  a⁻¹), propagating into the atmosphere and contributing to the similar trend that appears in the observed globally averaged δ(¹³C,CH₄). The subsequent increase after 2012 led to a stabilization of the associated atmospheric signal.

Between 2002–2007 and 2007–2014, all regional isotopic signatures were shifted downward (about −0.5 ‰ over the globe), except in China (+0.6 ‰). This is mainly explained by an increase in emissions from ¹³C-depleted sources in most of the regions but also by a decrease in emissions from ¹³C-enriched sources and a decrease in AGW and FFG source signatures.

Additionally, we use a simple mathematical framework to attribute the shift in the global flux-weighted signature to the different emission categories. Here, $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ denotes the global flux-weighted δ_source(¹³C,CH₄). A first-order estimate of $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ is given by

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{\mathit{\delta }}=\sum _{i=\mathrm{1}}^N\left({\mathit{\delta }}_i\times {\displaystyle \frac{f_i}{F}}\right),\end{array}$$

where f_i denotes the global CH₄ emissions from a specific category, N is the number of emission categories (five here), $F={\sum }_{i=\mathrm{1}}^N{f}_i$ is the total CH₄ emissions, and $r_i=\frac{f_i}{F}$ is the contribution from each category to the total emissions. A small variation in the global flux-weighted source signature, $\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\delta }}$, can therefore be calculated using the derivatives dδ_i and df_i:

$$\begin{array}{}\text{(6)}& \mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\delta }}=\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\delta }}}{\partial {\mathit{\delta }}_i}}\times \mathrm{d}{\mathit{\delta }}_i\right)+\sum _{i=\mathrm{1}}^N\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\delta }}}{\partial f_i}}\times \mathrm{d}{f}_i\right),\end{array}$$

with

$$\begin{array}{}\text{(7)}& \left\{\begin{array}{ll}\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\delta }}}{\partial {\mathit{\delta }}_i}& =\frac{f_i}{F}\\ \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\delta }}}{\partial f_i}& ={\mathit{\delta }}_i\times \frac{\left(\mathrm{1}-r_i\right)}{F}-\sum _{\begin{array}{c}j=\mathrm{1}\\ j\ne i\end{array}}^N\left({\mathit{\delta }}_j\times \frac{f_j}{F^{\mathrm{2}}}\right).\end{array}\right.\end{array}$$

Using this simplified linear relationship, we find that the 0.34 ‰ decrease in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ (see Fig. 8) between 2002–2007 and 2007–2014 was due to

1. a decrease in the global AGW source signature (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of −0.22 ‰),

2. a small decrease in BB emissions (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of −0.15 ‰),

3. a large increase in AGW emissions (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of −0.14 ‰),

4. a decrease in the FFG isotopic signature (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of −0.09 ‰).

This decrease is partially offset by

1. a large increase in FFG emissions (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of +0.24 ‰),

2. a small decrease in wetland emissions (resulting in a shift in $\stackrel{\mathrm{‾}}{\mathit{\delta }}$ of +0.02 ‰).

Figure 8 Contributions of the changes in CH₄ emissions and δ_source(¹³C,CH₄) source isotopic signatures to the global source signature shift between 2002–2007 (left side) and 2007–2014 (right side) in INV_REF. The upper part of the figure shows the contributions from the individual emission categories to the total emissions. Associated bars are non-transparent and non-hashed. Percentages and emissions are displayed on top of the bars. The lower part of the figure shows the isotopic signatures of each category. Associated bars are slightly transparent and hashed. The lower center with a white background shows the contributions from changes in emissions (non-transparent and non-hashed) and source isotopic signatures (slightly transparent and hashed) to the total source signature shift (−0.34 ‰) between the two periods. This part is magnified (×100) for clarity. Results from the other sensitivity tests are given in Table 5.

Table 5 Upper part of the table: changes in CH₄ emissions (emi.) and isotopic signatures (sign.) between 2002–2007 and 2007–2014 for each emission category and each sensitivity test. Lower part of the table: contributions from changes in emissions and isotopic signatures to the global flux-weighted source signature ($\stackrel{\mathrm{‾}}{\mathit{\delta }}$) shift between 2002–2007 and 2007–2014 for each emission category and each sensitivity test.

n/a: not applicable.

3.7 Sensitivity of the results to isotopic constraints

The reference inversion assimilates both CH₄ and δ(¹³C,CH₄). To quantify the impact of assimilating δ(¹³C,CH₄) data, INV_CH4 assimilates CH₄ observations only and does not simulate the isotopic composition. Differences between INV_REF and INV_CH4 therefore provide insight into the influence of the isotopic constraint. Notably, these two inversions show similar results for CH₄ emissions. As both inversions are constrained by the same global sink, global emissions estimated by INV_CH4 are only 0.3 Tg  a⁻¹ lower on average over the 2002–2014 period. Tropical emissions are increased compared to INV_REF, and the contribution from tropical emissions to the total emissions is shifted from 59.8 % to 60.3 % (+2.5 Tg  a⁻¹), mainly due to an increase in WET tropical emissions (+1.6 Tg  a⁻¹), which is offset by decreases in the northern mid-latitudes (−1.2 Tg  a⁻¹) and high latitudes (−0.2 Tg  a⁻¹). WET emissions increase by 1.1 Tg  a⁻¹ in INV_CH4 between 2002–2007 and 2007–2014 (Table 5) instead of decreasing in INV_REF. The increases in AGW, FFG and WET emissions contribute 50 %, 47 % and 3 % of the post-2007 renewed growth, respectively, and are therefore slightly different from the results of INV_REF. To summarize, adding the isotopic constraint (INV_REF as compared to INV_CH4) slightly decreases the contribution from tropical emissions to global emissions, removes a very small contribution from WET emissions to the post-2007 renewed growth, and slightly changes the contributions from emission increases to the post-2007 renewed growth.

Differences between INV_REF and INV_CH4 are small, presumably as a result of the large prior uncertainties in source signatures, which allows the inverse system to adjust the atmospheric isotopic compositions at a low cost by changing signatures rather than emissions. To test this hypothesis, we run INV_LOCKED assuming a perfect knowledge of isotopic signatures (no uncertainties in the prior source isotopic signatures). Although the global total emissions obtained with INV_REF and INV_LOCKED are very similar, the individual contributions from each emission category are modified (Fig. 9c). On average over the 2002–2014 period, FFG emissions are increased by 24 % compared to INV_REF, mainly due to large relative increases in China and the Middle East (+30 % to 50 %). Global FFG emissions amount to 153 Tg  a⁻¹, making them more consistent with the large revisions (150–200 Tg  a⁻¹) derived by Schwietzke et al. (2016) with recent isotopic data. WET emissions located in boreal regions and in South America are decreased by around 30 %, whereas WET emissions from Central Africa are slightly increased, leading to a decrease in global WET emissions of 14 %. Finally, BB emissions are increased by 41 % and AGW emissions are slightly decreased by 7 %, with globally uniform changes (Fig. 9a and b). Furthermore, INV_LOCKED explains the renewed growth through contributions from enhanced FFG (46 %), AGW (37 %) and WET (17 %) emissions between 2002–2007 and 2007–2014 (Table 5). WET emissions therefore actively participate in the post-2007 CH₄ growth in INV_LOCKED, as opposed to INV_REF. The FFG, AGW and WET emission increases are, however, offset by a large decrease in BB emissions, nearly three times larger than in the other inversions. Also, emission IAVs are increased for all categories. In particular, the BB emission peaks in 2006 and 2009 are much higher in INV_LOCKED than in INV_REF relative to the mean over the period. Such variations are probably too large to be realistic when compared to prior data and other inversion studies (e.g., Chandra et al., 2021; Basu et al., 2022). However, they provide an upper bound for emission trends as constrained by isotopic values. In many inversion studies (e.g., Rice et al., 2016; Schwietzke et al., 2016; Turner et al., 2017; Rigby et al., 2017; Thompson et al., 2018; McNorton et al., 2018), isotopic signatures are fixed, and our results suggest that this may lead to significant errors in CH₄ emission trends. This stresses the importance of finding the right balance between overconstrained signatures, as in INV_LOCKED, and likely underconstrained signatures, as in INV_REF. At present, either isotopic constraints are too loose to yield critical information about sectorial and regional CH₄ emissions or our estimates of the associated uncertainties are overestimated in our methodology, which is also a possibility that we will address in future studies.

Finally, assimilating δ(D,CH₄) observations and optimizing δ_source(D,CH₄) source signatures in INV_DD has a very small influence on our posterior emission estimates, as indicated in Table 5. The most significant difference observed is a small positive shift of +0.5 ‰ in the BB posterior source signature compared to INV_REF. Consequently, with our setups, assimilating δ(D,CH₄) does not appear to provide any substantial additional constraint on the CH₄ budget estimate. Several factors may contribute to this result: (1) the existing network provides comparatively few δ(D,CH₄) observations in comparison to δ(¹³C,CH₄) observations; (2) δ(D,CH₄) observations only span the period from 2005–2010 and so the full run cannot be fully constrained by this data; and (3) the constraints may be too weak due to an overestimation of the prescribed uncertainties in δ_source(D,CH₄) source signatures. As including δ(D,CH₄) in the inversion doubles the computational cost compared to a setup like INV_REF, we recommend that δ(D,CH₄) should not be assimilated in our system until either the computational cost can be reduced, more observations become available or lower uncertainties are established. However, a hypothetical network of δ(D,CH₄) measurements obtained at a reasonable frequency and spanning a longer period of time could efficiently complement δ(¹³C,CH₄) observations and provide a wealth of information (Rigby et al., 2012). More specifically, reactions with OH, O¹D and Cl have fractionation coefficients that depend on the isotope. Therefore, incorporating δ(D,CH₄) constraints might help to disentangle the effects of the associated sinks and provide additional insights into the global sink and its mixture. However, optimizing the sinks introduces additional degrees of freedom and complexifies the inverse problem. With the current system and at such a high resolution for the optimized variables, we recommend against the simultaneous optimization of both the source signatures and the sinks. However, a coarser resolution for the optimized variables, or at least for the sink, might be able to accommodate a simultaneous optimization.

Figure 9 Comparison between INV_REF and INV_LOCKED results. Upper panels show the time series of emissions estimated by INV_REF (a) and INV_LOCKED (b). For these panels, the time series show the anomalies around the 2002–2014 mean value. For each category in these panels, the associated mean value is displayed in the same color as the solid line. Lower panels show the contributions from each emission category to the total emissions for 2002–2007 and 2007–2014 (c) and a map of the posterior total emission differences between INV_LOCKED and INV_REF averaged over the 2002–2014 period (d).

3.8 Sensitivity of the results to OH IAV

Last but not least, we have tested the impacts of OH trends on our results. INV_FLATOH and INV_REF show very similar results (Table 5). The main difference is that INV_FLATOH infers a smaller increase in total emissions (+18 Tg  a⁻¹) between 2002–2007 and 2007–2014 than INV_REF does (+24 Tg  a⁻¹). As a smaller sink is prescribed in INV_FLATOH compared to INV_REF, the increase in total emissions required to fit the observations of CH₄ amount fractions is also smaller. The contributions from the emission categories to the total emissions are barely affected (± 0.2 %). The contributions from AGW and FFG emissions to the increase in the total emissions between the two periods are exactly the same as in INV_REF. In addition, as the inter-hemispheric OH ratio is not modified, the contributions from the tropics, mid-latitudes and high latitudes are identical for each emission category. Overall, the differences between INV_REF and INV_FLAT are negligible and do not affect the conclusions deduced from the INV_REF results.

On the contrary, INV_TURNER infers a decline of 1.6 % in global emissions between 2002–2007 and 2007–2014 (Table 5), mainly driven by a large decrease in WET emissions (−15 Tg  a⁻¹) and a slightly larger decrease in BB emissions (−4 Tg  a⁻¹) than in INV_REF. With this prescribed OH IAV, the renewed CH₄ growth post-2007 is therefore entirely caused by a large decline in the global OH sink between the two periods. Changes in AGW and FFG emissions in INV_TURNER are still positive but are 2–3 times smaller than in INV_REF. Using all the information provided by the sensitivity tests inferring a net increase in emissions (i.e., without INV_TURNER), this increase is principally attributed to fossil sources (50 ± 3 %) and agriculture and waste sources (47 ± 5 %). Nevertheless, there is substantial variation in the results between configurations that optimize source signatures and those that do not.

The decline in the OH sink between the two periods affects the δ(¹³C,CH₄) atmospheric signal in two opposite ways:

1. If the OH sink is the only sink, a decline in OH concentrations has no effect on δ(¹³C,CH₄) in the long term (several decades) because the mean fractionation is not affected. However, in the short term (a decade), as OH concentrations decrease, ¹²CH₄ and ¹³CH₄ atmospheric lifetimes increase. Due to the fractionation effect, there is a time lag between increases in ¹²CH₄ and ¹³CH₄ amount fractions. ¹²CH₄ accumulates faster than ¹³CH₄, leading to a decrease in δ(¹³C,CH₄).

2. The total fractionation effect in the atmosphere is the result of averaging all of the fractionation effects associated with the different sinks (OH, O¹D, Cl, soils) weighted by their contributions to the total sink. Therefore, if the OH sink is reduced, the contributions from the other sinks (with larger fractionation effects) increase. Consequently, the total fractionation effect is also increased and δ(¹³C,CH₄) values are shifted upward.

As INV_TURNER infers a more depleted global source signature (−0.50 ‰) than INV_REF (−0.34 ‰), we can conclude that the downward shift in δ(¹³C,CH₄) induced by the first mechanism is smaller than the upward shift induced by the second mechanism, resulting in a net upward shift. The enhanced depletion of the global source signature counterbalances this net upward shift. In the inversion, such a depletion is mainly obtained by lowering the source signatures of AGW, FFG and WET sources between 2002–2007 and 2007–2014. Compared to INV_REF, the shifts in source signatures are almost identical for AGW but much larger for FFG and WET. The negative OH trend was obtained by Rigby et al. (2017) and Turner et al. (2017) with box modeling and methyl-chloroform constraints (Patra et al., 2021). However, Naus et al. (2019) suggested that inter-hemispheric transport, stratospheric loss and source and sink spatial distributions are not properly represented using box modeling, resulting in significant errors. They found a positive OH trend over the 1994–2015 period with a 3-D model, which is more consistent with the IAV of our OH-INCA field. Other studies agree with Naus et al. (2019) in finding a small or positive IAV for recent years (Montzka et al., 2011; Lelieveld et al., 2016; Nicely et al., 2018). Therefore, the results from the INV_TURNER inversion seem to be rather unlikely.

3.9 Comparison with Basu et al. (2022)

Basu et al. (2022) – referred to as “BA22” in this subsection – quantify the global CH₄ budget and investigate the post-2007 renewed growth, using the TM5-4DVAR inversion framework to assimilate both CH₄ and δ(¹³C,CH₄) measurements. In our opinion, their work is strongly relevant and tackles this complex topic with an appropriate and robust methodology. As our goals are similar, we compare our systems and methodologies here.

First, it is worth mentioning that our system is capable of assimilating δ(D,CH₄)) observational data and optimizing the associated source signatures. Although this feature has a small influence on our results in this work, we believe that its relevance will grow as more δ(D,CH₄)) data become available and the associated uncertainties decrease.

As already stated in their paper, the main difference between our systems is the optimization of source signatures. BA22 prefer to investigate the influence of the source signature uncertainties with different sensitivity tests that adopt various source signature maps. This choice relies on the fact that they can run a large number of inversions at low cost using a parallel configuration. It is a good strategy to assess the influence of systematic errors in source signatures. In our work, we did not investigate this influence, and we decided to optimize source signatures in order to consistently account for random errors in source signatures and emissions at the same time. As both interact and impact the atmospheric composition in very complex ways, it seemed important to us to perform at least one inversion combining all the uncertainties.

We agree with BA22 that the second major difference between the two studies lies in the construction of the prior CH₄ fluxes. However, they suggest that we constructed a prior that approximately matches the atmospheric CH₄ growth rate in Thanwerdas et al. (2022a). In Thanwerdas et al. (2022a) and in the present study, we derived our prior fluxes and source signatures solely on the basis of bottom-up estimates and literature data. For CH₄ emissions, the fact that prior simulations match the atmospheric CH₄ growth rate shows that bottom-up estimates are roughly consistent with atmospheric data, even before the inversion process.

We prefer not to adjust prior fluxes to match observational δ(¹³C,CH₄) data because, in this case, we assume that bottom-up estimates suffer from strong systematic uncertainties, which is difficult to demonstrate. When adjusting prior CH₄ fluxes, it is also assumed that CH₄ emissions derived by bottom-up estimates are more likely to be wrong than source signatures estimates. We believe that the opposite is more plausible because observational data on source signatures are very scarce at present. Therefore, we prefer to start from robust and validated data and let the inversion system combine them with the assimilated atmospheric observations and the random uncertainties. BA22 start with a flat prior, and the posterior results deviate significantly from the prior. As the prior data do not seem to have a strong influence on the posterior results, adjusting prior fluxes prior to the inversion should have no effect on the results. However, this is yet to be confirmed in our case, i.e., with a nonlinear observation operator.

BA22 are able to calculate posterior uncertainties using a large ensemble of inversions. It is a precious feature that we do not possess at present. This is made possible by the relatively low computational cost of their configuration (an adjusted prior and a linear formulation) but also by the fact that they divide the full assimilation window into shorter sub-windows (5 years) that are run in parallel. A 1-year overlap with previous and next sub-windows is applied. It would be interesting to compare the posterior results obtained with this parallelized configuration to an inversion with a complete assimilation window. As the relaxation time for isotopic composition in the atmosphere in response to a perturbation is much larger (decades; Tans, 1997) than that for CH₄ itself, we are concerned that using such short time periods might affect posterior results, especially if the observation operator is nonlinear. Modifying the prior data to fit the observed isotopic composition, as in BA22, might be a prerequisite for the success of this method.

It is clear that these setup differences propagate to posterior results. Using the additional δ(¹³C,CH₄) data, BA22 find that fossil CH₄ emissions and microbial emissions contributed about 15 % and 85 %, respectively, of the post-2007 CH₄ growth. As presented in the previous sections, our results are completely different. Most notably, they find a contribution of 30 % from fossil emissions to the total emissions on average over 1999–2016. While our reference inversion finds a much smaller number (21 %), our inversion with fixed source signatures (INV_LOCKED) gives a closer value (26 %). The small source partitioning discrepancy between INV_LOCKED and BA22's inversion results might also be caused by a difference in prescribed isotopic fractionation, as suggested by the sensitivity analysis of Lan et al. (2021).

We cannot fully explain why BA22's conclusions about the causes of the renewed CH₄ growth post-2007 differ so substantially from our own. However, it appears that BA22 also use a robust methodology to study the global CH₄ budget and the renewed growth. Despite significant differences, we find good complementarity between our approaches and hope to learn from each other in order to improve our systems and reconcile our results.

4 Conclusion and discussion

We used variational inversion modeling with the 3-D CTM LMDz-SACS to investigate the drivers of the post-2007 renewed growth of atmospheric CH₄. We assimilated CH₄, δ(¹³C,CH₄) and δ(D,CH₄) atmospheric observations and optimized both the fluxes and the source isotopic signatures of five independent emission categories for the period 1998–2017. Implementing multiple setups allowed us to investigate the influences of isotopic constraints and OH IAV on our results.

Most of our inversions find that the post-2007 renewed growth was caused by large increases in fossil fuel and geological emissions (FFG) as well as in agricultural and waste (AGW) emissions between 2002–2007 and 2007–2014. The contributions from these two categories were almost equal (51 % for FFG and 49 % for AGW). These were partially balanced by small decreases in wetland (WET) emissions and biofuel- and biomass-burning (BB) emissions and a small OH increase during this period.

Isotopic constraints, i.e., those achieved by assimilating δ(¹³C,CH₄) and δ(D,CH₄) observations, have little influence on the posterior emission estimates. Compared to a CH₄-only inversion, an inversion assimilating δ(¹³C,CH₄) observations and optimizing source signatures only slightly reduces tropical emissions (−2.5 Tg  a⁻¹), mainly those from wetlands. Notably, the global flux-weighted WET source signature is shifted upward (less negative) due to a shift in the tropics (+3.2 ‰) and in the northern mid-latitudes (+7.1 ‰). To fit the δ(¹³C,CH₄) observations, the system prefers to adjust the source signatures rather than the CH₄ emissions. Undoubtedly, the large uncertainties associated with source signatures make them less costly to modify. Our findings also reveal that the global downward shift in δ(¹³C,CH₄) between 2002–2004 and 2007–2014 was caused by an increase in ¹³C-depleted AGW emissions and a decrease in ¹³C-enriched BB emissions but also by decreases in the AGW and FFG source signatures. For example, a small change compared to the uncertainties in the AGW source signatures (−0.6 ‰) between the two periods results in a −0.24 ‰ (downward) shift of the global source signature in the reference inversion. These results might be very dependent on the prescribed Cl concentrations, especially in the troposphere, and we decided to use the most recent and consistent Cl concentration estimates to minimize the associated error.

If the δ_source(¹³C,CH₄) source signatures are considered to be perfectly known, i.e., with no uncertainty, the relative contributions of the different emission categories are significantly changed by the inversion. Contributions from FFG and BB emissions are increased and those from AGW and WET emissions are decreased. In addition, WET emissions are found to contribute (13 %) to the post-2007 renewed growth, along with AGW (37 %) and FFG (46 %) emissions. Such a partition between fossil and microbial sources is more consistent with recent inversion estimates based on isotopic data. However, none of these recent results account for random uncertainties in source signatures. This shows that reducing the prescribed uncertainties in source signatures is a necessary condition for providing more accurate emission estimates when assimilating isotopic data.

OH IAV also has an influence on the results when a negative trend consistent with the IAV inferred by Turner et al. (2017) is applied. In this case, the post-2007 renewed growth is entirely caused by the decline in OH concentrations, and AGW and FFG emissions only slightly increase over the 2002–2014 period. As recent findings suggest that such a decrease in OH concentrations is unlikely, the results from the other sensitivity inversions should be preferentially considered. Overall, using all the information provided by the sensitivity tests presented in this work, the net increase in global emissions is principally attributed to fossil sources (50 ± 3 %) and agriculture and waste sources (47 ± 5 %).

As this new inversion setup (with isotopic constraints) was used over a long time period for the first time, methods were deliberately simplified in order to provide a background for future inversions and improvements. For instance, our methods to prescribe error statistics in the matrix B obviously have room for improvement, even with the limited amount of data available at the present time. The uncertainties we prescribed for source signatures in the reference inversion might be slightly overestimated. A more robust estimate of current regional random uncertainties in source signatures is necessary before running other inversions with isotopic data.

Also, the main limitation of our inversion system is the associated computational cost and the absence of posterior uncertainties. Formally, posterior uncertainties are given by the Hessian of the cost function (Meirink et al., 2008). This matrix is difficult to compute at an achievable cost considering the size of the inverse problem. Other means must be implemented to obtain the posterior uncertainty, such as estimating a lower-rank approximation of the Hessian using Monte Carlo ensembles of the variational inversion to represent the prior uncertainties (Chevallier, 2007). However, the amount of time required to run a single inversion is too large at present, preventing the derivation of robust posterior statistics as well as attempts to account for systematic uncertainties. Recent developments in the CIF (Chevallier et al., 2023; Chevallier, 2013) may help us to significantly reduce our computational costs and run Monte Carlo ensembles. While these new features have not been tested with realistic configurations yet, preliminary results are promising.

The inversion system proposed in this work benefits from the advantages of both 3-D modeling and variational inversion methods, and it also includes the optimization of the source isotopic signatures. Additionally, it accounts for the observation operator nonlinearity, which is an important component of isotopic data assimilation, particularly when source signatures are also optimized. To our knowledge, such a system is unique and allows us to reconcile emissions and source signatures with the limitation of still-large random uncertainties in the isotopic signatures. More developments are necessary to improve the robustness of the estimates and the relevance of such a system, but we believe that this study represents a significant step towards better quantification of the CH₄ sectorial and regional emissions and the global CH₄ budget.

Appendix A: Surface in situ CH₄, δ(¹³C, CH₄) and δ(D, CH₄) observation sites

Table A1 List of surface in situ CH₄ observation sites that provided measurements assimilated in the inversions between 1998 and 2018. AOC, PAO, POC and WPC are mobile stations. Their characteristics are compiled into a single line that provides latitude and longitude ranges of the measurements. Stations that retrieved samples consisting mainly of well-mixed marine boundary layer (MBL) air are indicated in bold.

Table A2 List of surface in situ δ(¹³C,CH₄) observation sites that provided measurements assimilated in the inversions between 1998 and 2018. WPC is a mobile station. Its characteristics are compiled into a single line that provides latitude and longitude ranges of the measurements. Stations that retrieved samples consisting mainly of well-mixed MBL air are indicated in bold.

n/a: not applicable.

Table A3 List of surface in situ δ(D,CH₄) observation sites that provided measurements assimilated in the inversion INV_DD between 2005 and 2010. Stations that retrieved samples consisting mainly of well-mixed MBL air are indicated in bold.

Appendix B: Additional results

Figure B1 Panel (a) shows a comparison between δ(D,CH₄) observations and prior and posterior simulations. Panel (b) shows linear regressions applied on the monthly and globally averaged δ(D,CH₄) observations. We performed three linear regressions: (1) one over the full data period from May 2005 to July 2009 (brown line), (2) one over the period from May 2005 to January 2007 (violet line), and one over the period from January 2007 to July 2009 (blue line). For each linear regression, the coefficient (a), its standard error and the Pearson's correlation coefficient (r) are displayed in the legend. Note that the x axis stops before 2010, as we have only selected months with sufficient data for the average to be representative of the whole globe.

Figure B2 Same as Fig. 4 but with prior data. Note that the scale for panel (c) has been modified and is not centered on zero anymore because the prior agreement with δ(¹³C,CH₄) data is too low. The large red-shaded area in panel (d) is caused by a change in the OH sink (INV_TURNER).

Figure B3 Comparison between prior–posterior and observed wetland δ_source(¹³C,CH₄). Observations are taken from Supplementary Data 1 provided by Oh et al. (2022). For each observation, prior and posterior values are sampled using the grid cell corresponding to the latitude and longitude provided in the dataset. Error bars for each observation point represent the observation uncertainty. Each panel shows a comparison with observations located in a selected region: the tropics (<30^∘ N/S), temperate (30–50^∘ N/S), western boreal (50–90^∘ N and <15^∘ W) and eastern boreal (50–90^∘ N and >15^∘ W). For each panel, the identity line and two linear fitting lines (prior in green and posterior in red) are displayed. The parameters of the fitting lines, the Pearson's correlation coefficients (r) and the root-mean-square error (RMSE) are given in the legend box.

Figure B4 Same as Fig.7 but for prior source signatures δ_source(¹³C,CH₄). For each panel, the time series show the anomalies around the 2002–2014 mean value. The units of variations and means are ‰. Note that x axis ranges from 1998–2018 to illustrate the effects of the spin-up and spin-down mentioned in Sect. 2.8. Also note that the regions used here are slightly different from the regions selected for the optimization.

Figure B5 Same as Fig. B4 but for posterior source signatures δ_source(¹³C,CH₄) from INV_REF. Note that the green (WET) and blue (NAT) lines are flat because (1) the prior signatures are constant over time, (2) these categories do not result from the aggregation of multiple subcategories and (3) we optimize only one scaling factor per region for the entire period. Therefore, these values do not vary with time. Also, note that BB source signatures vary only because the regions used here are slightly different from the regions selected for the optimization. Therefore, the flux-weighted average produces some temporal variability.

Table B1 Posterior CH₄ emissions for the globe and three different latitudinal bands averaged over 2002–2014 for all inversions and all categories. The units of the values shown are Tg  a⁻¹.