Estimating emissions of methane consistent with atmospheric measurements of methane and δ¹³C of methane

We have constructed an atmospheric inversion framework based on TM5-4DVAR to jointly assimilate measurements of methane and δ¹³C of methane in order to estimate source-specific methane emissions. Here we present global emission estimates from this framework for the period 1999–2016. We assimilate a newly constructed, multi-agency database of CH₄ and δ¹³C measurements. We find that traditional CH₄-only atmospheric inversions are unlikely to estimate emissions consistent with atmospheric δ¹³C data, and assimilating δ¹³C data is necessary to derive emissions consistent with both measurements. Our framework attributes ca. 85 % of the post-2007 growth in atmospheric methane to microbial sources, with about half of that coming from the tropics between 23.5^∘ N and 23.5^∘ S. This contradicts the attribution of the recent growth in the methane budget of the Global Carbon Project (GCP). We find that the GCP attribution is only consistent with our top-down estimate in the absence of δ¹³C data. We find that at global and continental scales, δ¹³C data can separate microbial from fossil methane emissions much better than CH₄ data alone, and at smaller scales this ability is limited by the current δ¹³C measurement coverage. Finally, we find that the largest uncertainty in using δ¹³C data to separate different methane source types comes from our knowledge of atmospheric chemistry, specifically the distribution of tropospheric chlorine and the isotopic discrimination of the methane sink.

1 Introduction

Current atmospheric levels of methane (CH₄) are about 2.5 times pre-industrial levels, primarily due to anthropogenic emissions (Dlugokencky et al., 2011). The main sources of CH₄ to the atmosphere today are known, which are periodically summarized by the Global Carbon Project (GCP, Saunois et al., 2020). In brief, they include anthropogenic sources from agriculture (ruminants, manure, and rice), waste management (landfills and waste treatment), fossil fuel production and use (coal, oil, and natural gas), and biomass burning (including biofuels). The remainder is from natural processes, predominantly tropical and high-northern-latitude wetlands, with smaller contributions from termites, wild animals, and geologic seeps. In the latest GCP report, however, a large disparity of ∼160 Tg yr⁻¹ remains between the bottom-up budget constructed from inventories and the top-down budget derived from atmospheric CH₄ measurements (Saunois et al., 2020), signifying gaps in our understanding of the CH₄ budget.

As shown in Fig. 1, CH₄ levels have been rising rapidly since 2007 after a period of relatively slow growth in 1999–2006 (Dlugokencky et al., 2011; Saunois et al., 2020). The mechanisms behind the relative stability of 1999–2006 and growth thereafter, however, are not yet fully understood. Possible mechanisms suggested in the literature include an approach to a steady state in the early 2000s (Dlugokencky et al., 2003), followed by an increase in either agricultural (Schaefer et al., 2016) or fossil (Worden et al., 2017) emissions or both (Saunois et al., 2020; Jackson et al., 2020), an increase in global (Schwietzke et al., 2016) or tropical (Nisbet et al., 2016, 2019) microbial emissions, a decrease in methane uptake by upland soils (Ni and Groffman, 2018), and decadal changes in the atmospheric sinks of methane (Rigby et al., 2017; Turner et al., 2017). It is difficult to choose between these competing explanations based on atmospheric CH₄ measurements alone. However, measurements of the ¹³C:¹²C ratio of CH₄, denoted δ¹³C−CH₄ or δ¹³C for short, provide some additional information to distinguish between these hypotheses (Lan et al., 2021).

Different CH₄ sources have distinct δ¹³C signatures over large spatial scales, and different sinks consume ¹²CH₄ and ¹³CH₄ at slightly different rates, imposing different signals on atmospheric δ¹³C (Miller, 2004). Therefore, atmospheric δ¹³C measurements can help constrain and refine the CH₄ budget. In an earlier publication, we described the simulation of atmospheric CH₄ and δ¹³C using the model TM5 (Krol et al., 2005) and its use for evaluating competing hypotheses about renewed CH₄ growth since 2007 (Lan et al., 2021). In this work, we construct and apply a variational inversion framework based on TM5 to assimilate CH₄ and δ¹³C measurements and estimate space- and time-varying emissions of CH₄ disaggregated by source type. With this framework, we perform atmospheric inversions from 1997 to 2016 to infer large-scale methane emissions from different sources, assess the added value of δ¹³C measurements compared to traditional CH₄-only inversions, and investigate the possible factors behind the post-2007 growth in atmospheric CH₄.

Several previous studies have used the information provided by δ¹³C measurements to infer mechanisms behind the behavior of atmospheric methane over the past 2 decades. However, many of these studies approximated the global atmosphere as a small number of connected boxes, with homogeneous emissions and chemistry in each box (e.g., Schwietzke et al., 2016; Schaefer et al., 2016; Nisbet et al., 2016, 2019; Worden et al., 2017; Zhang et al., 2021b). They were therefore susceptible to biases inherent in box models (Naus et al., 2019) and were unable to use the information contained in spatial gradients of atmospheric CH₄ and δ¹³C. Moreover, by construction, box models have to simplify the complexity of δ¹³C source signatures, transport variability, and loss processes and cannot extract information from spatial gradients in atmospheric measurements. In Sect. 4.3, we explore some of the limitations of box models and their differences from a 3D model by comparing our conclusions to those from the recently published work of Zhang et al. (2021b), which used a box model in combination with atmospheric CH₄ and δ¹³C data to infer causes behind the post-2007 growth in atmospheric CH₄. Some studies have used 3D atmospheric circulation models to estimate CH₄ emissions consistent with observed δ¹³C (e.g., Bousquet et al., 2006; Rice et al., 2016). However, they have generally used globally uniform δ¹³C source signatures, when in reality signatures of some of the most important sources such as wetlands and fossil fuels have strong latitudinal gradients and spatial variations. In this study, we confront our best estimate of spatiotemporally varying methane emissions and source signatures with a newly constructed multi-laboratory dataset of atmospheric CH₄ and δ¹³C measurements in the TM5-4DVAR framework. Our technique is analogous to a recently submitted paper by Thanwerdas et al. (2022), and in Sect. 4.4 we discuss the similarities and differences between our methods.

Figure 1 Global average CH₄ (blue circles, left axis) and δ¹³C (red diamonds, right axis) from NOAA marine boundary layer (MBL) and other background sampling sites. The gray box denotes the period from 1999 to 2006 when the atmospheric CH₄ burden was relatively stable, in contrast to the periods of growth before and after. Regular δ¹³C measurements started at NOAA background sites in 1998, which is the first year with an estimate of the global δ¹³C. The selection of marine boundary layer sites and the construction of global averages are described in detail by Masarie and Tans (1995).

2 Method

2.1 Formulation of the mass balance equations

The atmospheric mass balance of ¹²CH₄ can be written as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{}^{\mathrm{12}}\mathrm{C}& =\sum _s{}^{\mathrm{12}}{F}_{\mathrm{s}}-\sum _l{}^{\mathrm{12}}{k}_l\times {}^{\mathrm{12}}\mathrm{C}\times \left[l\right]\\ & =\sum _s{}^{\mathrm{12}}{F}_{\mathrm{s}}-{}^{\mathrm{12}}\mathrm{C}\sum _l{L}_l,\end{array}\end{array}$$

where F_s denotes the surface flux from source category s, and l denotes species contributing to chemical loss, namely Cl, OH and O(¹D). We classify all sources into three isotopically distinct source categories, namely pyrogenic, fossil, and microbial CH₄. While the upland soil sink of methane is included in the sources in conventional methane modeling, for reasons described in Sect. 2.4 we have modeled it as a fourth loss mechanism. The combination ¹²k_l×[l] can be denoted as a loss rate or inverse lifetime L_l due to species l. For ¹³CH₄, we can write a corresponding equation,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{}^{\mathrm{13}}\mathrm{C}=\sum _s{}^{\mathrm{13}}{F}_{\mathrm{s}}-{}^{\mathrm{13}}\mathrm{C}\sum _l{\mathit{\alpha }}_l{L}_l,\end{array}$$

where ${\mathit{\alpha }}_l={}^{\mathrm{13}}{k}_l{/}^{\mathrm{12}}{k}_l$. Using the definition of δ and denoting the atmospheric isotope ratio (signature of source s) as δ_a (δ_s), respectively, we can substitute ${}^{\mathrm{13}}\mathrm{C}={}^{\mathrm{12}}\mathrm{C}{r}_{\text{std}}({\mathit{\delta }}_{\mathrm{a}}+\mathrm{1})$ and ${}^{\mathrm{13}}{F}_{\mathrm{s}}={}^{\mathrm{12}}{F}_{\mathrm{s}}{r}_{\text{std}}({\mathit{\delta }}_s+\mathrm{1})$ and get

$$\begin{array}{}\text{(3)}& {}^{\mathrm{12}}\mathrm{C}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathit{\delta }}_{\mathrm{a}}=\sum _s\left({\mathit{\delta }}_s-{\mathit{\delta }}_{\mathrm{a}}\right){}^{\mathrm{12}}{F}_{\mathrm{s}}-{}^{\mathrm{12}}\mathrm{C}({\mathit{\delta }}_{\mathrm{a}}+\mathrm{1})\sum _l{\mathit{ϵ}}_l{L}_l,\end{array}$$

where ${\mathit{ϵ}}_l={\mathit{\alpha }}_l-\mathrm{1}$ and r_std=0.0112372 is a predefined standard ratio¹ (Craig, 1957). While Eqs. (1) and (3) are mathematically complete descriptions of the ¹²CH₄ and ¹³CH₄ budgets, they are not the most convenient form for constructing a dual-tracer CH₄ and δ¹³C inversion. This is because it is total CH₄ that is measured and not the two isotopologues separately. We therefore construct an alternate formulation in terms of ${\mathit{\delta }}^{\prime }=({}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CH}}_{\mathrm{4}})/r_{\text{std}}-\mathrm{1}$, which can be related to the more traditional $\mathit{\delta }=({}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}/{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}})/r_{\text{std}}-\mathrm{1}$ by

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\delta }}^{\prime }& {\displaystyle }={\displaystyle \frac{\mathrm{1}+\mathit{\delta }}{\mathrm{1}+r_{\text{std}}(\mathrm{1}+\mathit{\delta })}}-\mathrm{1},\text{(5)}& {\displaystyle }\mathit{\delta }& {\displaystyle }={\displaystyle \frac{\mathrm{1}+{\mathit{\delta }}^{\prime }}{\mathrm{1}-r_{\text{std}}(\mathrm{1}+{\mathit{\delta }}^{\prime })}}-\mathrm{1}.\end{array}$$

In terms of this δ^′, the mass balance equations become

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}C=\sum _s{F}_{\mathrm{s}}-C\sum _l{L}_l+r_{\text{std}}C({\mathit{\delta }}_{\mathrm{a}}^{\prime }+\mathrm{1})\sum _l{L}_l(\mathrm{1}-{\mathit{\alpha }}_l),\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}C{\mathit{\delta }}_{\mathrm{a}}^{\prime }& =\sum _s{\mathit{\delta }}_s^{\prime }{F}_{\mathrm{s}}-C{\mathit{\delta }}_{\mathrm{a}}^{\prime }\sum _l{\mathit{\alpha }}_l{L}_l+C\sum _l{L}_l(\mathrm{1}-{\mathit{\alpha }}_l)\\ & -r_{\text{std}}C({\mathit{\delta }}_{\mathrm{a}}^{\prime }+\mathrm{1})\sum _l{L}_l(\mathrm{1}-{\mathit{\alpha }}_l),\end{array}\end{array}$$

where $C={}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}+{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}$ and $F_{\mathrm{s}}={}^{\mathrm{12}}{F}_{\mathrm{s}}+{}^{\mathrm{13}}{F}_{\mathrm{s}}$ are total methane moles and fluxes, respectively. This reformulation of ¹³CH₄ abundance in terms of total carbon is similar to that by Tans et al. (1993). In Eq. (6), if we consider the coefficients of any L_l, then the second term supplies C∼1800 ppb, while the third term supplies $r_{\text{std}}C({\mathit{\delta }}_{\mathrm{a}}^{\prime }+\mathrm{1})(\mathrm{1}-{\mathit{\alpha }}_l)\sim \mathrm{0.086}$ ppb, approximating r_std=0.01, ${\mathit{\delta }}_{\mathrm{a}}^{\prime }=-\mathrm{0.05}$ and α_l=0.995. In Eq. (7), with the same approximations, the coefficients of L_l in the last three terms are respectively 89.5, 9, and 0.086 ppb. So in both equations, we ignore $r_{\text{std}}C({\mathit{\delta }}_{\mathrm{a}}^{\prime }+\mathrm{1}){\sum }_l{L}_l(\mathrm{1}-{\mathit{\alpha }}_l)$, leading to

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}C& {\displaystyle }\simeq \sum _s{F}_{\mathrm{s}}-C\sum _l{L}_l,\text{(9)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}C{\mathit{\delta }}_{\mathrm{a}}^{\prime }& {\displaystyle }\simeq \sum _s{\mathit{\delta }}_s^{\prime }{F}_{\mathrm{s}}-C{\mathit{\delta }}_{\mathrm{a}}^{\prime }\sum _l{\mathit{\alpha }}_l{L}_l+C\sum _l{L}_l(\mathrm{1}-{\mathit{\alpha }}_l).\end{array}$$

In this formulation, the two tracers to be simulated are total CH₄ (which is measured) and an artificial tracer $C{\mathit{\delta }}_{\mathrm{a}}^{\prime }$. All measurements of δ¹³C are converted to ${\mathit{\delta }}_{\mathrm{a}}^{\prime }$ before assimilation. Note that the tracer $C{\mathit{\delta }}_{\mathrm{a}}^{\prime }$ does not have any surface flux of its own. There is “production” at the surface proportional to the CH₄ surface flux and loss in the atmosphere. The loss reactions of Eqs. (8) and (9) are coupled, and the loss of the tracers from time t to t+δt is calculated by solving the differential equation to give

$$\begin{array}{}\text{(10)}& {\displaystyle }C(t+\mathit{\delta }t)=C\left(t\right)e^{-\mathit{\delta }t\sum _l{L}_l},\text{(11)}& {\displaystyle }\begin{array}{rl}C{\mathit{\delta }}^{\prime }(t+\mathit{\delta }t)& =\left[C\left(t\right)+C{\mathit{\delta }}^{\prime }\left(t\right)\right]e^{-\mathit{\delta }t\sum _l{\mathit{\alpha }}_l{L}_l}\\ & -C\left(t\right)e^{-\mathit{\delta }t\sum _l{L}_l}.\end{array}\end{array}$$

2.2 Inversion framework

We use the TM5-4DVAR inversion framework (Meirink et al., 2008), which has been used to estimate surface fluxes of CO, CO₂, and CH₄ (Hooghiemstra et al., 2011; Bergamaschi et al., 2013; Krol et al., 2013; Basu et al., 2013, 2014) in single-tracer inversions, as well as source-specific CO₂ fluxes in multi-tracer inversions (Basu et al., 2020, 2016; Ma et al., 2021). At the heart of the framework is the TM5 offline tracer transport model (Krol et al., 2005) and its adjoint, driven by ECMWF ERA-Interim reanalysis winds and run globally at 3^∘ × 2^∘ with 25 vertical layers defined by sigma–pressure hybrid coordinates. Two tracers are simulated in TM5: total methane or C of Eq. (8) and the artificial tracer $C{\mathit{\delta }}_{\mathrm{a}}^{\prime }$ of Eq. (9). Measurements of CH₄ are directly compared to modeled values of C, while measurements of δ¹³C are first converted to ${\mathit{\delta }}_{\mathrm{a}}^{\prime }$ and then to $C{\mathit{\delta }}_{\mathrm{a}}^{\prime }$ by multiplying with values of CH₄ mole fractions measured in the same air samples.

TM5-4DVAR minimizes a cost function J as a function of surface fluxes x (the set of all F_s of Sect. 2.1) by balancing fits to atmospheric observations y with deviations from the prior fluxes x₀,

$$\begin{array}{}\text{(12)}& J\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(Hx-y{)}^T{R}^{-\mathrm{1}}(Hx-y)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(x-x_{\mathrm{0}}\left)B^{-\mathrm{1}}\right(x-x_{\mathrm{0}}),\end{array}$$

where H is the transport, chemistry, and observation operator connecting surface fluxes with atmospheric measurements, and R and B are the error covariances of Hx−y and prior fluxes, respectively. Our formulation of R contains both the analytical measurement uncertainty and a model representativeness error proportional to local tracer gradients (Meirink et al., 2008). For each source type (pyrogenic, fossil, and microbial), the diagonal elements of B per time step and lateral grid cell are proportional to the prior flux, or f×x₀. Off-diagonal elements of B are constructed assuming an exponential decay of the prior error correlation in space and time with source-specific scales λ and τ, respectively. The values of f, λ, and τ for the different source types are given in Table 1. While there is no unique way of specifying these parameters, our choices yield reasonable 1σ prior uncertainties in global total microbial, fossil, and pyrogenic emissions of ∼25 Tg yr⁻¹ (∼7 %), ∼30 Tg yr⁻¹ (∼17 %), and ∼2 Tg yr⁻¹ (∼6 %), respectively. The 1σ uncertainty in the annual global total CH₄ emission comes to ∼40 Tg yr⁻¹ (∼7 %) with these choices. The cost function J of Eq. (12) is minimized over 50 iterations by a conjugate gradient minimizer utilizing the Lanczos algorithm (Lanczos, 1950; Courtier et al., 1998).

Table 1 Parameters for constructing the prior flux error covariance.

In TM5-4DVAR, we calculate J(x) of Eq. (12) with TM5 and $\partial J/\partial x$ with its adjoint. A traditional variational estimation would require us to run the forward and adjoint models multiple times over the entire period over which we want to estimate fluxes. However, these model runs require a significant amount of time, and iterations must be performed in succession. For example, at our 3^∘ × 2^∘ resolution, TM5 simulates a decade in 8 h. So to perform an inversion over 2 decades with 50 iterations (one iteration is one forward and one adjoint model run), it would take $\mathrm{8}\times \mathrm{2}\times \mathrm{2}\times \mathrm{50}=\mathrm{1600}$ h or 67 d just for the model runs, not counting time spent in the computing queue. This is impractical given the need to do tests required of any new inversion system. Therefore, we split up our target period into several inversions that were run in parallel as shown in Fig. 2a. A single forward run from 1984 to 2017 produced initial C and $C{\mathit{\delta }}_{\mathrm{a}}^{\prime }$ fields for all inversions. This forward run was identical to scenario C_WL+ of Lan et al. (2021) and matched the long-term atmospheric CH₄ and δ¹³C trends over that period. Six 5-year inversions were run simultaneously with 2 years of overlap (red bars) between inversions, starting in 1997, 2000, 2003, 2006, 2009, and 2012. After all six inversions finished, the fluxes from the middle 3-year period of each inversion (blue bars) were considered for analysis. For simulating prior and posterior mole fractions, fluxes from the non-overlapping periods (1997–2001, 2001–2004, 2004–2007, … 2013–2017) were stitched together and a single forward run was done with those fluxes.

Figure 2 A schematic of the time splitting of our inversions. Red bars denote spin-up and spin-down periods, and blue bars denote periods from which fluxes were considered in our analysis. Panel (a) denotes the time splitting used in most of our inversions as described in Sect. 2.2, whereas panel (b) denotes the time splitting used specifically with climatological priors as described in Sect. 3.4. In the splitting scheme in (a), each inversion spans 5 years, and the entire time span is covered with six inversions running simultaneously, starting from initial fields produced by a 1984–2017 model run with prior fluxes. In the splitting scheme in (b), each inversion spans 6 years, and the entire span is covered with four inversions. However, except for the 1997–2003 inversion, all other inversions are started from the optimized mole fraction fields at the end of year 5 of the previous inversion, and therefore the inversions cannot be run in parallel.

2.3 Prior fluxes and δ¹³C source signatures

The prior fluxes and their δ¹³C source signatures for the different categories of methane emissions are described in detail as scenario C_WL+ in Lan et al. (2021). Briefly, the prior fluxes are based on bottom-up emission estimates with adjustments to match global atmospheric CH₄ increases and to satisfy the global mass balance of δ¹³C over 1984–2017. For biomass burning or pyrogenic emissions, we use the Global Fire Emission Database (GFED) 4.1s for 1997–2016 (van der Werf et al., 2017) and estimates from the Reanalysis of Tropospheric chemical composition (RETRO) project before 1997 (Schultz et al., 2008). Other anthropogenic emissions are based on the EDGAR 4.3.2 inventory (Janssens-Maenhout et al., 2019). We use natural fossil emissions reported by Etiope et al. (2019). Emission estimates from wild animals and termites are adopted from Houweling et al. (1999) and Sanderson (1996), respectively. Wetland emissions and upland soil consumption of methane are estimated by a process-based model (Zhuang et al., 2004; Liu et al., 2020), after which the soil sink is modeled as a first-order loss process as explained in Sect. 2.4.

The δ¹³C source signatures used in our study are mainly spatially resolved maps based on the Global δ¹³C Source Signature Inventory 2020 for coal, oil and gas (ONG), biomass and biofuel burning, ruminant and wild animal sources (Sherwood et al., 2021; Lan et al., 2021), spatial maps for geological seeps (Etiope et al., 2019), and wetland sources (Ganesan et al., 2018). Globally averaged values are used for waste, landfills, termites, rice, and other energy and industry, given insufficient data to develop spatial distributions for their δ¹³C signatures (Lan et al., 2021).

The sum of the bottom-up methane emission estimates described above is not consistent with top-down estimates of global total emissions based on observed atmospheric CH₄ growth and estimated loss, which requires a 46 Tg yr⁻¹ increase in the annual global emission in 2016 compared to the 1999–2006 quasi-stable period. In addition, the δ¹³C mass balance requires 167 Tg yr⁻¹ emissions from fossil sources (including natural geological seeps) to be consistent with modeled sinks and the δ¹³C source signatures described above. Therefore, we (i) scale the ONG emissions from EDGAR 4.3.2 uniformly using annual scaling factors to reach a total of 167 Tg yr⁻¹ from all fossil sources, (ii) impose a linear trend on wetland emissions to achieve an increase of 46 Tg yr⁻¹ in total 2016 emissions compared to 1999–2006, and (iii) adjust emissions from agricultural and waste sectors to match the year-on-year global CH₄ growth rate derived from marine boundary layer observations (Dlugokencky et al., 2011). This ensures that our global CH₄ and ¹³CH₄ budgets approximate the long-term trends in atmospheric CH₄ and δ¹³C over 1984–2017.

2.4 Methane loss mechanisms and fractionation

Atmospheric methane has four loss mechanisms: atmospheric oxidation by OH and Cl throughout the atmosphere, destruction by O(¹D) in the stratosphere, and surface uptake by upland soils. In all our inversions, these sinks are prescribed and not optimized. Monthly climatological CH₄ loss rates in the stratosphere due to OH, Cl, and O(¹D) were constructed from a run of the ECHAM5/MESSy1 chemistry transport model (Steil et al., 2003; Jöckel et al., 2006). Loss due to tropospheric Cl is simulated using a recent model-derived estimate of tropospheric Cl (Hossaini et al., 2016). For tropospheric OH, we use the monthly OH climatology of Spivakovsky et al. (2000) after scaling by 0.9 to match the declining atmospheric abundance of methyl chloroform in the early 2000s (Montzka et al., 2011).

In most CH₄ inversions, upland soil sinks are folded into the net wetland flux. However, the soil sink fractionates strongly between ¹²CH₄ and ¹³CH₄ (King et al., 1989), and therefore we keep it separate from wetland fluxes. We model the soil sink as a first-order reaction at the surface, in which the loss rates of ¹²CH₄ and ¹³CH₄ in the lowest model layer are $\mathrm{d}{}^{\mathrm{12}}\mathrm{C}/\mathrm{d}t=-k_{\text{ss}}{F}_{\text{ss}}{}^{\mathrm{12}}\mathrm{C}$ and $\mathrm{d}{}^{\mathrm{13}}\mathrm{C}/\mathrm{d}t=-{\mathit{\alpha }}_{\text{ss}}{k}_{\text{ss}}{F}_{\text{ss}}{}^{\mathrm{13}}\mathrm{C}$, respectively. F_ss is the prior soil sink map from the Terrestrial Ecosystem Model (TEM), and k_ss is an arbitrary constant tuned to a value such that in a forward run with prior fluxes, the global total soil sink matches the prior total.

The fractionation between ¹²C and ¹³C for each of the loss reactions is modeled as $k_{\mathrm{12}}/k_{\mathrm{13}}=\mathrm{1}/\mathit{\alpha }=C{e}^{D/T}$ (Saueressig et al., 2001), where T is the air temperature in Kelvin. The soil sink fractionation is cast in a similar form for convenience. Coefficients C and D that we used are tabulated in Table 2.

Saueressig et al. (2001)Saueressig et al. (2001)Saueressig et al. (2001)King et al. (1989)

Table 2 Fractionation parameters for CH₄ loss.

2.5 CH₄ and δ¹³C measurements

To maximize the spatiotemporal coverage of in situ CH₄ and δ¹³C data, we have developed a new database by harmonizing measurements from NOAA/INSTAAR with those from 30 other laboratories around the world (Lan et al., 2021). All CH₄ data have been quality-checked and converted to a common CH₄ scale, namely the World Meteorological Organization (WMO) X2004A scale maintained at NOAA's Global Monitoring Laboratory (Dlugokencky et al., 2005). For data not on the WMO X2004A scale, we applied lab-specific scale multipliers estimated based on (i) comparisons of measurements of common air samples during the WMO/IAEA Round Robin Comparison Experiment (Crotwell et al., 2020) and (ii) comparisons of co-located atmospheric measurements made by NOAA and other laboratories. We constructed the uncertainty in the assimilated CH₄ measurements from a combination of (i) measurement repeatability of a single sample (hereafter called the single measurement precision), (ii) lab-specific long-term reproducibility based on analyzer type and sampling frequency reported in the literature, and (iii) each lab's realization of the calibration scale. If a scale conversion was needed to bring measurements onto the WMO X2004A scale, the mole fraction uncertainty due to the scale multiplier uncertainty was added in quadrature. The final uncertainties are typically less than 9 ppb for all CH₄ measurements.

We used δ¹³C data from the Institute for Arctic and Alpine Research (INSTAAR) as well as other isotope laboratories making precise measurements of atmospheric methane with isotope ratio mass spectrometers. The INSTAAR δ¹³C data were measured in a subset of air samples collected from NOAA's Global Greenhouse Gas Reference Network (GGGRN). Because different labs have independent ties to primary reference materials which do not agree, we calculated offsets to bring the δ¹³C data onto the INSTAAR realization of the Vienna Pee Dee Belemnite (VPDB) scale (Miller et al., 2002). These offsets were based on measurements of cylinders, flasks filled from cylinders, or co-located sample data, and are all described in Umezawa et al. (2018). When there was not a direct comparison, e.g., between INSTAAR and TU or INSTAAR and NIPR, we used comparisons between each of these labs and the Institute for Marine and Atmospheric research Utrecht (IMAU). Each comparison had an uncertainty associated with it, and these were combined in quadrature to account for uncertainty in the offset correction. The total uncertainty in assimilated δ¹³C measurements was typically less than 0.15 ‰. The final database of assimilated CH₄ and δ¹³C measurements is available at https://doi.org/10.15138/64w0-0g71 (Lan et al., 2022).

With the following exceptions, we assimilate all the observations from this database including marine boundary layer sites, surface and tower sites over continents (Andrews et al., 2014), and vertical profiles from routine aircraft measurements (Sweeney et al., 2015). Intermittent aircraft profiles such as from the HIPPO (Wofsy, 2011) and ATom (Thompson et al., 2022) campaigns are not assimilated. CH₄ data from flasks taken aboard routine flights between Japan and Australia as part of the CONTRAIL program have been assimilated (Machida et al., 2008; Matsueda et al., 2015). A subset of the CONTRAIL flasks was also analyzed for δ¹³C (Umezawa et al., 2012), which were not assimilated. For continental tower sites with multiple intake heights, only data from the highest intake are considered in inversions to minimize local influence. For sites with continuous CH₄ analyzers, the CH₄ data are averaged hourly, and only hourly averages between 11:00 and 16:00 local solar time are assimilated; these are the times when planetary boundary layer heights are likely to be best represented by transport models. For continuous CH₄ analyzers on mountain tops, we only assimilate hourly averages between 00:00 and 05:00 local solar time to avoid possible up-slope contamination. Site-specific statistical filtering based on a nonparametric curve fitting routine (Thoning et al., 1989) is further applied, with the exception of vertical profiles, to remove large outliers with potential local or other contamination. The number of CH₄ and δ¹³C measurements assimilated each year is summarized in Table 3, and their locations are plotted in Fig. 3.

Table 3 The number of CH₄ and δ¹³C observations assimilated in our inversions broken down by year.

Figure 3 Locations of assimilated CH₄ and δ¹³C measurements. The symbol sizes represent the number of measurements between 1997 and 2017 assimilated from each location. Overlapping symbols over some of the locations are due to multiple agencies measuring at those locations.

2.6 Uncertainty estimation and sensitivity tests

The uncertainty of surface emission estimates is a combination of random and systematic uncertainties. Random uncertainties are associated with components of the inversion system whose errors are assumed to be zero on average. In the formulation of the cost function (12), the prior flux x₀ is assumed to have a probability density function (PDF) centered on the true flux with variance around the truth given by the prior covariance matrix B. Similarly, the model–observation mismatch Hx–y is assumed to have a PDF centered around the mismatch between the true atmospheric mole fraction and true fluxes propagated through an unbiased transport model, with variance around this mean given by R. The random uncertainty in the optimal estimate is given by $\widehat{B}$, where

$$\begin{array}{}\text{(13)}& {\widehat{B}}^{-\mathrm{1}}={\displaystyle \frac{{\partial }^{\mathrm{2}}J}{\partial x^{\mathrm{2}}}}=H^T{R}^{-\mathrm{1}}H+B^{-\mathrm{1}}.\end{array}$$

Variational inversion systems such as TM5-4DVAR can construct a low-rank approximation of $\widehat{B}$ during the optimization. However, for large state vectors the $\widehat{B}$ thus constructed is an overestimation of the true posterior uncertainty (Meirink et al., 2008; Bousserez et al., 2015). We therefore construct an estimate of $\widehat{B}$ by performing an ensemble of 100 independent inversions for each of the 5-year inversions of Fig. 2a, with prior fluxes and observations perturbed according to the covariances specified by B and R, respectively. With 100 ensemble members, our estimate of $\widehat{B}$ is expected to be within 10 % of the exact analytical solution for $\widehat{B}$ (Bousserez et al., 2015). Furthermore, our ensemble of inversions allows us to compute any posterior covariance and correlation between estimated fluxes, such as between large regions or different CH₄ source types.

Systematic uncertainties are associated with aspects of the inversion system that are assumed to be fixed and perfectly known in principle but might in fact be biased in practice. In our inversion system, such aspects include, but are not limited to, atmospheric transport and chemistry, isotope source signatures, and the wetland inundation maps used to construct the prior wetland emissions. Because the posterior covariance estimate does not include systematic errors, we explore the impact of such errors by performing inversions with different realizations of potentially biased inputs in the following sensitivity tests.

2.6.1 Tropospheric chlorine

The magnitude and distribution of the Cl sink in the troposphere are uncertain, with estimates as high as 13–37 Tg yr⁻¹ based primarily on Southern Hemisphere background observations (Allan et al., 2007). However, more recent studies have found a more limited role of tropospheric Cl as a methane oxidant (Gromov et al., 2018). Consequently, most CH₄ inverse models neglect tropospheric Cl as a methane oxidant. However, due to the strong isotopic fractionation in the CH₄ + Cl reaction, Cl plays an important role in determining atmospheric δ¹³C (Strode et al., 2020; Lan et al., 2021). It is therefore important to test the sensitivity of our conclusions to the imposed tropospheric Cl sink within the range of realism. The Cl estimate by Hossaini et al. (2016) we use in this study is on the higher side of the range posited by Gromov et al. (2018). We perform an inversion with the tropospheric Cl field reported by Wang et al. (2021) as an alternative lower specification. In order to keep the global CH₄ lifetime unchanged between the two scenarios of tropospheric Cl, we scale the tropospheric OH field by 0.9 and 0.92, respectively, when we use the Cl fields of Hossaini et al. (2016) and Wang et al. (2021). Since the two scenarios lead to slightly different sink fractionation in the atmosphere, prior ONG and ruminant fluxes are adjusted to match the long-term atmospheric δ¹³C trend for both cases.

2.6.2 OH fractionation

We use chemical fractionation factors reported by Saueressig et al. (2001) since they provide factors for all atmospheric sink processes from a consistent set of laboratory measurements. While these are the most recent and generally accepted, for CH₄ oxidation by OH another set of coefficients (C=1.0054, D=0) have previously been reported by Cantrell et al. (1990). To the best of our knowledge, this earlier result has not been refuted in the literature, nor is there any independent evidence supporting one set of coefficients over another. Instead, the most recent evaluation of atmospheric reaction rates (Burkholder et al., 2019) recommends using the Saueressig et al. (2001) rates with increased uncertainty in the OH fractionation to include Cantrell et al. (1990) as a possibility. Since the sink fractionation plays a significant role in determining atmospheric δ¹³C, we perform an additional inversion with the OH fractionation of Cantrell et al. (1990) to gauge its impact. Since the two OH fractionation factors lead to different sink fractionation in the atmosphere, prior ONG and ruminant fluxes are adjusted to match the long-term atmospheric δ¹³C trend for both cases.

2.6.3 δ¹³C source signatures

In principle, it is possible to estimate both CH₄ fluxes and δ¹³C source signatures in a dual-tracer inversion (Thanwerdas et al., 2022). However, this makes the problem nonlinear and the inversion convergence slow. It is also difficult to construct a prior covariance for δ¹³C source signatures since much of the uncertainty stems from extrapolating a limited number of δ¹³C signature measurements to the entire domain of CH₄ sources, resulting in errors that are systematic and non-Gaussian. Therefore, we explore the impact of δ¹³C signature uncertainty on our results by running inversions with alternate specifications of δ¹³C signature maps as follows.

Source signature maps for biomass burning were calculated by multiplying C₃ and C₄ signatures of −26.7 ‰ and −12.5 ‰, respectively (Cerling et al., 1998), with the C₃ $/$ C₄ fraction for each 1^∘ × 1^∘ latitude–longitude grid cell (Lan et al., 2021). For ruminants and wild animals, C₃ and C₄ signatures were taken to be −54.5 ‰ and −67.8 ‰, respectively, from the Global δ¹³C Source Signature Inventory 2020 (Sherwood et al., 2021; Lan et al., 2021). In this way, the C₃ $/$ C₄ vegetation distribution determines the source signatures of both biomass burning and ruminant emissions. Our default inversion averages the C₃ $/$ C₄ distributions of Still et al. (2003) and its modified version as used by Randerson et al. (2012). To explore the uncertainty from the assumed C₃ $/$ C₄ map, we perform two additional inversions with δ¹³C source signature maps derived separately from the two individual C₃ $/$ C₄ distributions. In addition, country-level ruminant emission signatures were compiled by Chang et al. (2019), including their temporal changes due to shifting ruminant diet and due to the downward trend in atmospheric δ¹³CO₂ that is photosynthesized by the vegetation. We use the ruminant CH₄ source signatures of Chang et al. (2019) in a third inversion. The three instances of source signatures related to the C₃ $/$ C₄ distribution described here were significantly different, requiring us to adjust the prior flux apportionment to meet our goal of matching long-term CH₄ and δ¹³C trends. Specifically, we changed the prior fossil CH₄ emissions from the default of 167 to 158 Tg yr⁻¹ for the inversions using δ¹³C signatures derived from Still et al. (2003) and Chang et al. (2019). For the inversion using δ¹³C signatures derived from Randerson et al. (2012), we adjusted the prior fossil emission to 175 Tg yr⁻¹. In all cases, this was achieved by globally scaling the ONG and ruminant emissions to achieve long-term CH₄ and ¹³CH₄ mass balance.

For the global maps of ONG and coal emission signatures, our default inversion assumes time-invariant maps over the study period. However, considering the rapid development of US shale gas production and a shift in production from conventional to shale gas in the past decades, we estimate that the mean US ONG signature (production-weighted mean of shale and conventional gas) increased by 2.7 ‰ from 2006 to 2016 (Lan et al., 2021). We incorporate this in an alternate specification of fossil CH₄ source signatures and perform an inversion with this new map.

Finally, our default inversion setup uses the latitude-based source signature specification of Ganesan et al. (2018) for wetland emissions. Over the past several years we have implemented carbon isotopes in the TEM land surface model (Zhuang et al., 2004), making it possible to derive process-based δ¹³C wetland source signatures consistent with wetland emissions (IsoTEM, Oh et al., 2022). We perform an inversion with wetland source signatures from the IsoTEM model as an alternative to our default wetland source signatures.

2.6.4 Wetland inundation extent

Wetland inundation extent is a leading driver of uncertainty in bottom-up estimates of wetland CH₄ emissions and therefore in the global CH₄ budget. We explore this uncertainty by performing inversions with prior wetland CH₄ fluxes derived from the TEM model (Zhuang et al., 2004) driven by two different inundation maps. Our default setup uses a time-varying or dynamic inundation map based on the satellite-based Surface WAter Microwave Product Series (SWAMPS, Schroeder et al., 2015) combined with the Global Lakes and Wetlands Dataset (GLWD, Lehner and Döll, 2004; Poulter et al., 2017). In addition, we also drive the TEM model with the static inundation map of Matthews and Fung (1987), in which case meteorology is the only source of seasonal and interannual variation of prior wetland emissions. These two inundation maps produce significantly different atmospheric CH₄ and δ¹³C gradients (scenarios C_WL+ and Q_static_WL of Lan et al., 2021) in a forward run and therefore serve as a robust test of our inversion results with different inundation extents.

2.6.5 Initial δ¹³C gradients

Large-scale gradients of atmospheric δ¹³C take significantly longer to respond to changes in emissions compared to gradients of CH₄ (Tans, 1997), requiring multi-decade spin-ups for models trying to simulate atmospheric δ¹³C (Lan et al., 2021). Inverse models, on the other hand, take significantly less time to be spun up since fluxes during the spin-up period are modified to fit observed atmospheric δ¹³C. The exact spin-up duration required depends on the accuracy of the initial modeled δ¹³C gradients and the inversion setup. To test if a 1-year spin-up for our inversions as depicted in Fig. 2a is sufficient, we perform two additional inversions with different starting δ¹³C large-scale gradients. Specifically, of the flux scenarios simulated by Lan et al. (2021), we choose scenarios H_mean_sig and Q_static_WL, which produced the flattest and steepest north–south gradients in δ¹³C, respectively (see Lan et al., 2021, Fig. 5). We perform inversions starting from CH₄ and δ¹³C fields provided by forward simulations of those scenarios at each of the starting points in Fig. 2a. The resultant spread in fluxes provides an estimate of the sensitivity of our setup to erroneous initial δ¹³C gradients.

Figure 4 Observed (gray circles) and posterior modeled (colored symbols) CH₄ time series at three NOAA baseline observatories. Inversions both with and without δ¹³C data fit the CH₄ data equally well.

Figure 5 Observed (gray circles) and posterior modeled (colored symbols) δ¹³C time series at three NOAA baseline observatories. The inversion with δ¹³C data fits the observations throughout the inversion period, but the inversion without δ¹³C data – a traditional CH₄ inversion – drifts away from the observations with time. Note that both inversions were started with the same CH₄ and δ¹³C fields in 1997, but by the time δ¹³C data were available in mid-1998 they had already drifted apart, leading to the apparent initial offset in the plots above.

Figure 6 Modeled minus observed δ¹³C as a function of latitude and altitude from the ATom (top), HIPPO (middle), and CONTRAIL (bottom) aircraft campaigns. Mismatches are shown for the prior flux, the CH₄-only inversion, and the dual-tracer CH₄ + δ¹³C inversion. Altitudes have been binned in 1 km bins, while latitudes have been binned either in 10^∘ (ATom, HIPPO) or 5^∘ (CONTRAIL) bins. The rightmost panels show the number of samples averaged per bin.

3 Results

3.1 Fit to atmospheric CH₄ and δ¹³C data

Both the CH₄-only and the CH₄ + δ¹³C inversions fit the atmospheric CH₄ data, while only the latter is consistent with atmospheric δ¹³C data. This is demonstrated both at surface sites from which data were assimilated and data from aircraft campaigns that were withheld for validation. Figure 4 shows that both inversions fit the observed CH₄ time series at three NOAA baseline observatories. However, despite starting from realistic atmospheric CH₄ and δ¹³C fields, the CH₄-only inversion moves progressively farther from observed δ¹³C with time at those same locations, as demonstrated in Fig. 5. Only the CH₄ + δ¹³C inversion fits both atmospheric CH₄ and δ¹³C data. This is also demonstrated in Fig. 6, which compares modeled δ¹³C to δ¹³C measured by the HIPPO and ATom aircraft campaigns, and from regular flights between Japan and Oceania as part of the CONTRAIL program. ATom and HIPPO campaigns sampled primarily background air over the oceans at multiple latitudes and altitudes, and neither CH₄ nor δ¹³C data from those campaigns were assimilated. CONTRAIL primarily sampled the marine background at multiple altitudes as well, except for a small number of samples taken during takeoff and touchdown in Japan. CH₄ flask samples from CONTRAIL were assimilated in both inversions, but their δ¹³C measurements were not assimilated. The CH₄-only inversion compares far less favorably to the δ¹³C measurements than the joint inversion. Therefore, it is reasonable to conclude that our CH₄-only inversion, and very likely most traditional CH₄-only inversions, does not yield a CH₄ emission distribution consistent with atmospheric δ¹³C observations. We therefore expect our CH₄ + δ¹³C inversion to provide more accurate emission estimates and source partitioning than our CH₄-only inversion.

3.2 Large-scale fluxes from CH₄ and CH₄ + δ¹³C inversions

The top row of Fig. 7 shows the global total annual emissions from two inversions, a CH₄-only or “traditional” methane inversion without δ¹³C data, and a joint CH₄ + δ¹³C inversion developed in this work. The shaded regions in Fig. 7 denote 2σ random errors derived from 100-member Monte Carlo ensembles of inversions described in Sect. 2.6. Annual averages of the emissions and random errors are summarized in Table 4. Since many methane studies calculate emissions between and outside the ±30^∘ latitude band, in Table 4 we report those numbers as well for ease of comparison. The global total emission from all categories is unaffected by the addition of δ¹³C data, since δ¹³C does not place any additional constraint on the total CH₄ emission. However, the partitioning between microbial and fossil sources is changed significantly with the addition of δ¹³C data. Based on comparison to atmospheric data as noted in Sect. 3.1, we expect the source partitioning from our CH₄ + δ¹³C inversion to be more accurate compared to our traditional CH₄-only inversion.

Figure 7 also shows the total and source-disaggregated CH₄ emissions from our CH₄ and CH₄ + δ¹³C inversions over three latitude bands; the tropics are bounded between 23.5^∘ S and 23.5^∘ N. Relative to the prior, tropical (extratropical) total emissions are adjusted upward (downward) by both inversions, and there is little sensitivity of the tropical versus extratropical partitioning to the assimilation of δ¹³C data. In the northern extratropics, the partitioning of CH₄ emissions between the different source types does not change significantly with the addition of δ¹³C data. However, in the tropics the inversion with δ¹³C data shows significantly higher fossil (and lower microbial) emissions than the inversion without δ¹³C data. Fossil CH₄ emissions in the southern extratropics are significantly different for most years in the presence of δ¹³C data, but similarly significant differences do not exist for the other source types. Finally, our estimate of pyrogenic emissions does not change significantly in Fig. 7 in the presence and absence of δ¹³C.

Figure 7 Total and source-specific annual emissions of CH₄ globally and from three latitudinal bands. “Tropics” in this context refers to the region between 23.5^∘ N and 23.5^∘ S, while the northern and southern extratropics are to the north and south, respectively. The shaded regions denote 2σ prior and posterior error bars.

Table 4 Annual averages of CH₄ emissions between 1999 and 2016 as well as their 2σ uncertainties shown in Fig. 7 and three additional zonal bands (in  Tg yr⁻¹). The CH₄-only and CH₄+δ¹³C inversions of Fig. 7 have been abbreviated to CH₄ and “joint”, respectively. Note that (a) the 2σ uncertainties tabulated are the averages of the uncertainties across 18 years, not the uncertainties in the 18-year average emissions, and (b) zeroes below are due to rounding to the nearest integers, not because the concerned quantities are exactly zero.

3.3 Systematic errors in emission estimates

As explained in Sect. 2.6, we estimate possible biases in our flux estimates by running the inversion with different choices of non-optimized input. The spread in annual emissions due to alternate specifications of atmospheric chemistry (tropospheric chlorine of Sect. 2.6.1 and OH fractionation of Sect. 2.6.2) is shown in Fig. 8. Analogous spreads due to different specifications of δ¹³C source signatures (Sect. 2.6.3), wetland inundation maps (Sect. 2.6.4), and initial atmospheric δ¹³C fields (Sect. 2.6.5) are shown in Fig. 9. Note that the y-axis ranges in Figs. 8 and 9 are different. The average spread in annual emissions from different latitude bands and source types is summarized in Table 5 for each sensitivity test. The average of the annual posterior uncertainties as depicted in Fig. 7 is also provided in Table 5 as “MC-derived (2σ)” for reference, with the caveat that 2σ uncertainties are not directly comparable to the range across a few inversions.

Figure 8 Total and source-specific annual emissions of CH₄ globally and from three latitudinal bands as in Fig. 7. The shaded regions denote the spread (max to min) of annual emissions from sensitivity tests described in Sect. 2.6.1 and 2.6.2. A stronger OH fractionation makes the atmosphere heavier, requiring a larger (smaller) fraction of microbial (fossil) emissions to match the same observations. A smaller chlorine sink makes the atmosphere lighter, requiring a smaller (larger) fraction of microbial (fossil) emissions to match the same observations. That is, a stronger OH fractionation and a smaller chlorine sink affect the fossil–microbial partitioning in opposite ways, resulting in the structure seen here.

Figure 9 Total and source-specific annual emissions of CH₄ globally and from three latitudinal bands as in Fig. 7. The shaded regions denote the spread (max to min) of annual emissions from sensitivity tests described in Sect. 2.6.3, 2.6.4, and 2.6.5.

Most of the sensitivity tests have little impact on the global total CH₄ emission, and the spread in the total CH₄ emission from different latitude bands is generally smaller than the posterior uncertainty of our base inversion. However, by far the largest source of error in partitioning the total emission into fossil and microbial sources comes from our representation of atmospheric chemistry, namely the distribution of tropospheric chlorine and the kinetic isotope effect of CH₄ destruction by OH. Unless the uncertainty in these two factors can be reduced, our ability to use δ¹³C measurements to partition different source types will be seriously hampered. The uncertainty arising from our limited knowledge of δ¹³C source signatures, to the extent represented by the different signature maps used, is lower than the uncertainty due to atmospheric chemistry. Lastly, the uncertainty due to an incorrect specification of the initial atmospheric δ¹³C field is minimal, in line with our expectation that an inversion will rapidly correct for it by adjusting emissions during its spin-up period. However, we note here that the “incorrect” initial fields we constructed for the last test still satisfied the global δ¹³C mass balance by construction. The sensitivity to an incorrect initial condition will likely be higher if the initial field does not satisfy global δ¹³C mass balance.

Table 5 Average annual Monte Carlo-derived random uncertainty and possible bias in CH₄ emissions, separated by source type (Tot: total, Mic: microbial, Fos: fossil, Pyr: pyrogenic) and latitude bands as in Fig. 7. For each source type and region, the mechanism behind the largest possible bias has been demarcated by highlighting the bias in bold. The MC-derived numbers are 2σ posterior uncertainties, and all other numbers represent the range between maximum and minimum estimates. All numbers are in teragrams (Tg) of CH₄ per year.

3.4 Attribution of the post-2007 methane growth

As discussed earlier (Fig. 1 and discussion in Sect. 1), the atmospheric methane burden has been steadily growing since 2007 after a period of quasi-stability during 1999–2006. We use our CH₄ + δ¹³C inversion to ask whether the addition of δ¹³C data can provide information on the sources of the additional methane. Figure 7 suggests that the trend in CH₄ emissions comes largely from microbial emissions in a CH₄ + δ¹³C inversion. However, it is possible that this attribution to microbial emissions comes from our prior – which had a trend in the microbial emissions and a temporally flat fossil contribution – instead of the atmospheric data. To assess the robustness of our inferred microbial and fossil emission trends, we perform a second set of inversions with the following modifications.

1. We construct climatological prior fluxes and source signatures by averaging our prior emissions and signatures from 2000 to 2006. Neither the resulting priors nor the source signatures have any time trend.

2. Since the methane budget from climatological priors is no longer in balance with the atmospheric growth, we cannot use the overlapping inversions of Fig. 2a to run multiple periods in parallel. Instead, we run four 6-year inversions in sequence, spanning 1997–2003, 2002–2008, 2007–2013, and 2012–2018, following the scheme shown in Fig. 2b. The first inversion used the same initial field in 1997 as our default inversion. Every successive inversion used the previous inversion's fifth year mole fraction field as an initial condition. The last year of each inversion is discarded in the end, and the first 5 years' fluxes are stitched together and analyzed.

The posterior uncertainties of the emissions derived from this modified setup are calculated by performing a Monte Carlo suite of 100 inversions as described in Sect. 2.6. The Monte Carlo runs follow the geometry of Fig. 2b as well, with the ith inversion (i=1 to 100) of each period initialized from the fifth year mole fraction field of the ith inversion of the previous period. This allows us to calculate not only annual uncertainties but also uncertainties in long-term averages.

To study the transition around 2007, we considered two periods: 2000–2006 and 2008–2014. Average total and source-specific emissions over the two periods are shown in Fig. 10, as is the change in the average emissions between the two periods. The prior fluxes do not change between the two periods, and therefore the estimated change must be driven by the atmospheric observations. Both the CH₄-only and the CH₄ + δ¹³C inversions estimate a change in the total emission of 27.1±0.6 Tg yr⁻¹ to match the increase in the atmospheric burden. However, while the CH₄-only inversion attributes ∼70 % of that to fossil CH₄ emissions and only ∼29 % to microbial emissions, the addition of δ¹³C data switches the balance to ∼15 % fossil and ∼85 % microbial emissions. This change in the allocation of the methane emission in the presence of δ¹³C data is significant compared to the uncertainties in the changes as depicted in Fig. 10. The contribution of pyrogenic emissions to the change is small in both inversions, and its change between the two inversions is not significant compared to its uncertainty. This is consistent with the downward trend in the global average δ¹³C in Fig. 1, since microbial sources are the lightest of the three source types.

Figure 10 Change in global CH₄ emissions between the periods 2000–2006 and 2008–2014, total (top left) and disaggregated by source type. The gray bars denote prior emissions, and the colored bars denote two inversions, one with and the other without assimilated δ¹³C data. For each source type, the first two columns show the average emission over the two periods in question, and the third column shows the change between the two periods. The 1σ error bars are derived from a 100-member Monte Carlo ensemble of inversions following the configuration of Fig. 2b.

Geographically, the change between the two periods is driven almost equally by the tropics and the northern extratropics (Fig. 11). In the tropics, the addition of δ¹³C data results in higher microbial emissions in both periods. The change between the two periods is also attributed to microbial emissions, unlike a CH₄-only inversion, which attributes the change primarily to fossil methane. In the northern extratropics, although the presence of δ¹³C data points to an increase in microbial emissions between the two periods, the relative apportionment of the increase between fossil and microbial emissions does not differ significantly from the CH₄-only emission if we consider the respective uncertainty estimates. Therefore, the largest contribution to the global increase in microbial emissions between the two periods (Fig. 10) comes from the tropics.

Figure 11 Similar to Fig. 10, but disaggregated into different latitude bands. As in Fig. 7, “tropics” refers to the region between 23.5^∘ S and 23.5^∘ N. Pyrogenic emissions have not been plotted because of their small contribution in all latitude bands.

It is worth noting here that a change in emission strengths is not the only possible mechanism for an increase in atmospheric CH₄; a reduction in the sink strength could also induce a positive trend in atmospheric CH₄ post-2007. However, Lan et al. (2021) have shown that the changes in sinks proposed so far in the literature to explain the post-2007 CH₄ growth are not consistent with the observed δ¹³C trend post-2007. We therefore do not consider those alternatives here.

3.5 Separating microbial and fossil emissions

The CH₄ observations assimilated in a CH₄-only inversion constrain the total CH₄ emission, and any source disaggregation relies on spatiotemporal separation of emissions as encoded in the prior emissions and their uncertainties. Since the two largest CH₄ source types (microbial and fossil) have different δ¹³C source signatures, assimilating δ¹³C observations should provide additional information to separate the two sources compared to a CH₄-only inversion. We can evaluate this additional information by looking at the posterior correlation between microbial and fossil emissions both globally and regionally. Posterior correlations between global annual microbial and fossil CH₄ emissions, calculated from our 100-member ensemble of independent inversions as described in Sect. 2.6, are shown in Fig. 12. Error bars on the correlations shown in Fig. 12 represent the 95th percentile range of 20 000 evaluations of the correlation by randomly sampling the 100-member inversion ensemble with replacement (Efron and Tibshirani, 1994).

For all the years shown in Fig. 12, a CH₄-only inversion results in a strong negative correlation between global microbial and fossil emissions, consistent with the idea that atmospheric CH₄ measurements constrain the total CH₄ budget much better than individual source types. The addition of δ¹³C data reduces this negative correlation, implying that δ¹³C provides additional information to disentangle different CH₄ source types. The degree of disentanglement, represented by the reduction in the negative correlation, is determined by the δ¹³C measurement coverage in a particular year and atmospheric transport connecting the emissions to those measurements. The correlation reduction is limited in our inversions by the relative sparsity of δ¹³C measurements; even in the most recent 2012–2017 inversion period, only 2.8 % of CH₄ measurements have corresponding δ¹³C measurements, overwhelmingly in locations far removed from significant CH₄ emissions (Fig. 3). Having more δ¹³C measurements in general, and specifically closer to emissive regions, should allow further disentangling of the different CH₄ source types.

Figure 12 Posterior correlation between global annual microbial and fossil CH₄ emissions (a) as well as over the northern extratropics (b) and Asia (c) for the two inversions of Fig. 7.

Over smaller regions, only the northern extratropics and Asia show significant decorrelation between annual fossil and microbial emissions (Fig. 12b and c) with the addition of δ¹³C data. While several other regions show similar reductions, the reductions are typically not significant compared to the 95th percentile error bars. The significant decorrelation seen for northern extratropical and Asian emissions may be because most δ¹³C measurements are in the northern extratropics and downwind of Asia in the Pacific. To see similar significant decorrelation over other regions we will likely need increased δ¹³C coverage closer to those regions. Although Fig. 3 shows some δ¹³C measurements over North America and Europe, the majority of those measurements are from the background air sampling sites Niwot Ridge and Jungfraujoch, respectively, and therefore do not contribute to significant decorrelation of fossil and microbial emissions from those continents.

4 Conclusions and discussion

We have constructed a variational atmospheric inversion system capable of assimilating CH₄ and δ¹³C measurements to estimate source-specific methane emissions within the TM5-4DVAR framework. We have assimilated CH₄ and δ¹³C measurements from a multi-agency air sampling network in this framework to estimate fossil, microbial, and pyrogenic emissions of atmospheric CH₄ globally. We have derived Bayesian uncertainty estimates for our emissions (random error) and investigated the impact of biases from non-optimized aspects of our inversion (systematic error). In Sect. 4.1 we summarize the main lessons learned from our work. In Sect. 4.2 to 4.5 we compare our work and results with several other estimates in the literature. In Sect. 4.6, we outline several planned additions and improvements to our framework in the future. Finally, in Sect. 4.7 we present areas of progress that we think are required in order to better use atmospheric δ¹³C data to disentangle the methane budget.

4.1 Enumerated conclusions

First, Figs. 5 and 6 show that our inversion assimilating only CH₄ does not yield a CH₄ emission distribution consistent with atmospheric δ¹³C. This is very likely true of CH₄ inversions in general, since they have no constraints forcing them to match atmospheric δ¹³C gradients and trends. Starting from a prior emission distribution consistent with atmospheric δ¹³C trends does not ensure that the posterior emission estimates will remain consistent. Our CH₄-only inversion started from a prior that reproduced the global mean atmospheric δ¹³C trend (scenario C_WL+ of Lan et al., 2021), yet the posterior deviated from it as shown in Figs. 5 and 6. We conclude that the only way to guarantee a posterior emission distribution consistent with both atmospheric CH₄ and δ¹³C data is to assimilate them simultaneously.

Second, given an atmospheric sink scenario, our current observational coverage allows us to estimate the global total CH₄ emission with a 2σ random uncertainty of ∼3.8 Tg yr⁻¹, which is less than 1 % of the total emission. Microbial, fossil, and pyrogenic emission uncertainties are around 3.5 %, 8 %, and 10 %, respectively, at the global scale. Given these posterior uncertainties, there are significant differences between inversions with and without δ¹³C data in the apportionment of the total CH₄ emission between microbial and fossil sources both globally and in the tropics (Fig. 7 and Table 4). In both regions, the inclusion of δ¹³C data in an inversion results in a significantly higher proportion of fossil emissions compared to microbial emissions, which we consider realistic since it matches both atmospheric CH₄ and δ¹³C data (Figs. 4 and 5). Pyrogenic emissions are relatively insensitive to the inclusion of δ¹³C data.

Third, we tested the sensitivity of our results to several factors that can lead to biases or systematic errors, as detailed in Sect. 2.6. This included different maps of the δ¹³C isotopic source signatures, static and dynamic maps of the wetland inundation extent, different initial δ¹³C fields, different fractionation factors for the CH₄+OH oxidation mechanism, and different fields of tropospheric Cl. The last two factors had by far the largest impacts on the large-scale apportionment between microbial and fossil emissions, even though their impact on the total CH₄ budget was nil or negligible. With the OH fractionation of Cantrell et al. (1990), the global microbial emission increases to 414 Tg yr⁻¹ and the fossil emission drops to 131 Tg yr⁻¹. With the lower estimate of tropospheric Cl from Wang et al. (2021), the global microbial emission decreases to 345 Tg yr⁻¹, while the fossil emission increases to 199 Tg yr⁻¹. Since some CH₄ inversions in the literature do not simulate a tropospheric Cl sink of CH₄, we tested the impact of this limiting case as well. In the absence of a tropospheric Cl sink, the global microbial emission drops further to 331 Tg yr⁻¹ and the fossil emission increases to 213 Tg yr⁻¹. Most of these shifts in the global partitioning are accompanied by shifts in the latitudinal partitioning. All of these are significant revisions to the partitioning of Table 4, suggesting that the ability of atmospheric δ¹³C measurements to partition the total CH₄ emission into different source types, at least over large regions, is limited by our knowledge of these two critical chemical processes. The uncertainty in our knowledge of δ¹³C source signatures, long considered a limitation on the use of δ¹³C data, is almost never a leading driver of uncertainty in Table 5, although it is usually more significant than either inundation extent or the initial δ¹³C field. Finally, our tests suggest that the impact of an incorrect initial δ¹³C field can be ameliorated by a relatively short spin-up of 1 year in an inversion, in contrast to a multi-decadal spin-up necessary for a forward model run.

Fourth, atmospheric δ¹³C data strongly suggest that the rise in microbial emissions is the primary driver of the post-2007 growth in atmospheric CH₄. While a CH₄-only inversion starting from priors without a time trend attributes ∼70 % of the growth to fossil emissions, the addition of δ¹³C data shifts that to microbial emissions being responsible for ∼85 % of the growth. Since the latter inversion is consistent with atmospheric δ¹³C data while the former is not (Fig. 5), we consider a majority microbial contribution to the post-2007 growth to be more realistic. A disaggregation of the growth by latitude bands suggests that a significant majority of the increase in tropical methane emissions is due to microbial and not fossil emissions. Moreover, although some of the sensitivity tests of Sect. 2.6 lead to different partitioning between microbial and fossil emissions, they all suggest a steeper trend in microbial compared to fossil emissions in Figs. 8 and 9.

Fifth, the ability of δ¹³C data to disentangle different CH₄ source types can be quantified by the reduction in the posterior correlation between emissions from those sources owing to the addition of δ¹³C data compared to a CH₄-only inversion. Considering the two largest source types of methane (microbial and fossil), we see significant reductions in their posterior correlation over the globe as well as the northern extratropics and Asia. The degree of decorrelation, however, is limited, and we do not see significant decorrelation over other regions. We hypothesize that this is not a limitation of our understanding of δ¹³C but rather of its limited observational coverage. Even in the most recent years less than 3 % of assimilated CH₄ measurements were accompanied by δ¹³C measurements, almost exclusively from background sites. It is very likely that an increase in the observational coverage of δ¹³C, preferably close to source regions, will improve the capability of δ¹³C measurements to distinguish between different CH₄ source types.

Sixth, while it is difficult to compare our emission budget directly with GCP due to different partitioning schemes, we note that our fossil fuel emissions for both the 2000–2009 and 2008–2016 periods are higher than the GCP top-down and bottom-up emissions. However, our estimate of the change in fossil fuel emissions between the two periods is significantly lower than the GCP estimates. Concurrently, our estimate of the change in microbial emissions over the same time is significantly higher than the GCP top-down estimate. Both of these discrepancies are driven by atmospheric δ¹³C data, since our CH₄-only inversion provides changes that are consistent with GCP estimates. We therefore conclude that the microbial and fossil emission change estimates in the GCP budget are consistent with atmospheric CH₄ data but not with δ¹³C data. Finally, our pyrogenic emission estimates are consistent with or close to the GCP estimates for both periods.

4.2 Comparison to the GCP methane budget

The Global Carbon Project (GCP) periodically publishes top-down and bottom-up budgets of methane emissions from a suite of models. However, a meaningful comparison between our emissions and the 2020 GCP budget (Saunois et al., 2020) is not straightforward. The GCP bottom-up (BU) budget for 2008–2017, with 737 Tg yr⁻¹ emissions and 625 Tg yr⁻¹ sinks, significantly overestimates the atmospheric growth rate. This is primarily due to an overestimate of both microbial (159 Tg yr⁻¹ freshwater sources) and fossil (45 Tg yr⁻¹ geologic sources) methane in the GCP budget, making a direct comparison with our microbial and fossil estimates meaningless. The GCP top-down (TD) estimates do not provide a fossil–microbial split of “other natural” emissions, also making a direct comparison with our estimates difficult. However, it is possible to calculate emissions for certain GCP categories from our inversions for some limited comparisons. For the period 2000–2009, we compare directly by computing the 2000–2009 average from our inversions shown in Fig. 7. However, since our inversions stop in 2016 and Table 3 of Saunois et al. (2020) reports the 2008–2017 average, we first calculate the 2008–2016 mean from the GCP budget using their reported 2008–2017 mean and their 2017 emissions. For computing the GCP means quoted here, we use the spreadsheet at https://doi.org/10.18160/GCP-CH4-2019 in Saunois et al. (2020), which may result in small differences from the numbers quoted in Table 3 of Saunois et al. (2020) due to rounding.

Our estimate of fossil fuel CH₄ emissions is significantly higher than GCP estimates, consistent with earlier work by Schwietzke et al. (2016) on δ¹³C-based source apportionment. Subtracting 35 Tg yr⁻¹ of geologic emissions from our fossil emission estimates, we arrive at 137±2 Tg yr⁻¹ of fossil fuel emissions during both the 2000–2009 and 2008–2016 periods, which is significantly higher than the GCP BU (TD) estimates of 111 and 127 Tg yr⁻¹ (99 and 109 Tg yr⁻¹), respectively. Recent measurements of ¹⁴CH₄ in ice cores suggest that geologic methane emissions are much lower than we have assumed, with an upper 95th percentile confidence limit of 5.4 Tg yr⁻¹ and mean of 1.6 Tg yr⁻¹ (Hmiel et al., 2020). If true, this would push our fossil fuel emission estimates even higher to 170±2 Tg yr⁻¹, indicating a significant underestimate in the GCP budget.

Assuming that methane emissions from geological seeps do not change significantly over decadal timescales, we estimate a change of 0±2 Tg yr⁻¹ in fossil fuel CH₄ emissions from 2000–2009 to 2008–2016. This is markedly different from the BU (TD) GCP budget, which estimates an increase of 16 Tg yr⁻¹ (10 Tg yr⁻¹) between the two periods. While our baseline estimates for the two periods may be influenced by systematic biases (Sect. 2.6), the change between the two periods is relatively robust. With the alternate specification of tropospheric chlorine (Wang et al., 2021) and alternate fractionation due to the OH oxidation (Cantrell et al., 1990), which are the two biggest sources of bias in source apportionment by δ¹³C, the change in our fossil fuel emission estimate between the two periods is 1 and −1.7 Tg yr⁻¹, respectively, which is well within our uncertainty estimate of 2 Tg yr⁻¹ and significantly lower than both the GCP BU and GCP TD estimates. The GCP BU (TD) estimate of an increase of 16 Tg yr⁻¹ (10 Tg yr⁻¹) between the two periods is closer to our estimate of 8.4±5.6 Tg yr⁻¹ from a CH₄-only inversion, which is not consistent with δ¹³C data.

For reasons mentioned above, we cannot directly compare our microbial emission estimates to GCP emission estimates. However, if we assume that methane from termites, wild animals, and oceans does not change over decadal timescales, we can compare the change in the GCP TD estimate of wetlands, agriculture, and waste from 2000–2009 to 2008–2016 with the change in microbial emissions in our inversion estimates. The GCP TD budget estimates a change of 12.6 Tg yr⁻¹ between those two periods compared to our estimate of 26±2 Tg yr⁻¹ from a joint δ¹³C and CH₄ inversion and of 18±6 Tg yr⁻¹ from a CH₄-only inversion. Thus, the change in microbial emissions in the GCP TD budget is at the lower end of but consistent with our estimate from a CH₄-only emission, while it is not consistent with our budget after incorporating δ¹³C data. We cannot perform a similar analysis with the GCP BU budget because freshwater emissions cannot be assumed to be static over decadal timescales.

Finally, our pyrogenic emission estimates for both the 2000–2009 and 2008–2016 periods are 30±0.6 Tg yr⁻¹, with a change of 0.3±0.5 Tg yr⁻¹. These are close to the GCP BU (TD) estimates of 31 and 30 Tg yr⁻¹ (29 and 31 Tg yr⁻¹), respectively. Neither the GCP budgets nor our inversion show significant changes in pyrogenic methane emissions between the two periods.

4.3 Comparison with Zhang et al. (2021b)

Zhang et al. (2021b) derive sector-specific changes in methane emissions from a two-box model by first creating an ensemble of emission scenarios consistent with known emissions of methane and their uncertainties, then comparing the resultant hemispheric mean CH₄ and δ¹³C in the atmosphere with measurements. Based on the ensemble members that agree with hemispheric mean CH₄ and δ¹³C time series, they come up with sector-specific emission changes. We note here that according to Figs. 1b and 2a of Zhang et al. (2021b), none of their choices fit the atmospheric δ¹³C data very well. We hypothesize that this is due to both the biases in their bottom-up emission estimates, which could not be sufficiently corrected due to tight assumed uncertainties in them, and the inability of a box model to interpret spatial gradients of CH₄ and δ¹³C. Specifically, Fig. 1b of Zhang et al. (2021b) shows that none of the ensemble members can reproduce the recent downward trend in δ¹³C. Therefore, their estimated fractional contribution of microbial (fossil) emissions to the recent growth in CH₄ is likely to be too low (high).

Zhang et al. (2021b)Zhang et al. (2021b)

Table 6 Comparison of the attribution of the recent growth of atmospheric CH₄ between Zhang et al. (2021b) and this work. The numbers show the change in emissions from a baseline period of 2000–2006 to two different periods, 2007–2012 and 2013–2017. All numbers are in teragrams (Tg) of CH₄ per year. Uncertainties in the Zhang et al. (2021b) estimates were obtained by personal communication with Zhen Zhang.

Table 6 summarizes the changes in methane emissions from the “quasi-stable” period of 2000–2006 to two later periods during the renewed growth, 2007–2012 and 2013–2017. For both the periods, we estimate significantly smaller contributions from fossil emissions compared to Zhang et al. (2021b). While the microbial contributions look similar in magnitude, they constitute a smaller fraction of the total change for Zhang et al. (2021b). Lastly, while Zhang et al. (2021b) estimate a downward trend in pyrogenic emissions, we do not estimate a trend significantly different from zero. To summarize, we estimate a larger (smaller) influence of increasing microbial (fossil) emissions behind the recent CH₄ growth, which is likely to be more accurate because it better reproduces the observed δ¹³C trend in the atmosphere (compare Fig. 5 here with Fig. 1b of Zhang et al., 2021b). We cannot separate the anthropogenic from natural microbial contributions in our framework, and we note that such a separation by Zhang et al. (2021b) relies on their modeled wetland fluxes and not on atmospheric CH₄ or δ¹³C data.

4.4 Comparison with Thanwerdas et al. (2022)

Thanwerdas et al. (2022) describe an alternative variational inversion framework using the LMDz-SACS model to assimilate CH₄ and δ¹³C measurements. We find it heartening that others have decided to tackle this complicated problem. Since they reserve decadal dual-tracer inversions for future work, we will compare their technique with ours to highlight the similarities and differences. The biggest difference lies in the decision of Thanwerdas et al. (2022) to optimize δ¹³C source signatures compared to our choice of keeping them fixed for a specific inversion. While δ¹³C source signatures are uncertain for many methane sources, we explain our reasons for not optimizing them in Sect. 2.6.3. Instead, we explore the impact of source signature uncertainty with different constructions of the source signature map as detailed in Sect. 2.6.3. In the end, at least for large geographical regions, the uncertainty from source signatures did not prove to be a leading uncertainty (Table 5). The second major difference between the two inversion frameworks lies in the construction of the prior CH₄ fluxes. While Thanwerdas et al. (2022) use a prior that approximately matches the atmospheric CH₄ growth rate, we construct our priors to match both the CH₄ growth rate and the δ¹³C trend over 2 decades. We suspect this and the linearity of our formulation due to not optimizing source signatures to be the reasons why our inversion required a shorter spin-up time compared to Thanwerdas et al. (2022).

There are also a few differences in implementation between the two frameworks. Most notably, Thanwerdas et al. (2022) estimate the posterior uncertainty as the spread between different inversion configurations, correctly stating that an evaluation of the posterior covariance matrix would require significantly more computing resources. We evaluate that posterior covariance matrix for both CH₄ + δ¹³C and CH₄-only inversions and present both types of uncertainty, namely the systematic uncertainty as the spread between multiple inversion configurations and the random (Bayesian) uncertainty as the spread of an ensemble of 100 independent inversions. The configurations we explore for the systematic uncertainty are also different from Thanwerdas et al. (2022) and include alternate specifications of the Cl oxidant and the isotopic discrimination of the CH₄ + OH reaction. We find the latter two to be the most significant drivers of uncertainty for partitioning CH₄ emissions using δ¹³C data.

While our implementation of the inversion is different from Thanwerdas et al. (2022), our goals are very similar. We look forward to long-term inversions of CH₄ and δ¹³C data using LMDz-SACS so that we may compare and contrast with our results presented here and figure out how best to use isotopic measurements to solve the atmospheric methane puzzle.

4.5 Comparison with other top-down studies using CH₄ and δ¹³C

In addition to the publications compared in detail above, some others have used combinations of CH₄ and δ¹³C observations to infer sector-specific sources of methane and the causes behind the recent increase, either in 2D (Thompson et al., 2018) or 3D (McNorton et al., 2018) atmospheric models. In terms of attributing the recent growth of CH₄, the most significant difference between our work and either of those studies is that our inferred change in biomass burning emissions is smaller. For example, our GFED 4.1s prior has a reduction of 1.35 Tg yr⁻¹ in pyrogenic emissions from 2003–2006 to 2007–2015 (van der Werf et al., 2017). This is changed to a reduction of 0.81 Tg yr⁻¹ by our CH₄ + δ¹³C inversion (Fig. 7) and to a reduction of 0.19±0.72 Tg yr⁻¹ when we start from a prior without trends (Fig. 10). This suggests that the estimate of reduced pyrogenic emissions is driven primarily by its presence in the prior and not by δ¹³C data, which are sparse near regions with large fire emissions. To contrast, McNorton et al. (2018) infer a reduction of 2.9 Tg yr⁻¹ in pyrogenic emissions between the same periods, allowing for a larger increase in fossil emissions to balance CH₄ and δ¹³C trends. Finally, we note that even though the INV_FIXED inversion of McNorton et al. (2018) with climatological priors allows for variations in OH and is therefore not directly comparable to our analogous inversion, their INV_FIXED estimates no reduction in pyrogenic emissions (Table 6 of McNorton et al., 2018). This supports our hypothesis that a reduction in pyrogenic emissions is driven primarily by the prior and not by atmospheric δ¹³C data.

In agreement with our estimates, Thompson et al. (2018) also attribute the majority of the growth between 2007 and 2014 to microbial and not fossil sources. They infer a 3±2 Tg yr⁻¹ reduction in pyrogenic emissions during that period compared to a reduction of 0.58 Tg yr⁻¹ inferred by our CH₄+δ¹³C inversion (0.50±1.66 Tg yr⁻¹ reduction if we start with a climatological prior). We note that our prior pyrogenic emission for 2014 is 3.45 Tg yr⁻¹ higher than 2007, and a CH₄-only inversion estimates an increase of 1.65 Tg yr⁻¹ between the two years (Fig. 7); therefore, even the small reduction we estimate must be driven by δ¹³C data. However, since this reduction is smaller than the 3±2 Tg yr⁻¹ estimated by Thompson et al. (2018), our fossil emission increase is also smaller to balance the δ¹³C trend. Finally, Thompson et al. (2018) do not estimate a significant role of the chemical sink behind the recent CH₄ growth, consistent with our previous work (Lan et al., 2021).

4.6 Future work

While we feel confident in the CH₄ emission estimates reported here, there are several areas which we plan to explore and improve in future work.

4.6.1 Alternate OH

The atmospheric CH₄ budget is determined by the balance between its sources and sinks, the latter primarily driven by the OH radical. While there have been some efforts to optimize atmospheric OH in concert with CH₄ emissions (e.g., Zhang et al., 2018, 2021a; Yin et al., 2021), we do not think in situ CH₄ samples provide sufficient information to constrain the sink independently. Moreover, estimates of OH abundance and variability over the past decades, either from CH₄ inversions (Yin et al., 2021) or otherwise (Bousquet et al., 2005; Montzka et al., 2011; Nicely et al., 2018), are consistent with a limited role of OH variability in recent trends in atmospheric CH₄. This is why, similar to the vast majority of CH₄ inversions, we have chosen to keep the OH sink fixed to a field consistent with observed trends and gradients of methyl chloroform (MCF, Spivakovsky et al., 2000; Patra et al., 2014, 2020). Nonetheless, we acknowledge that our knowledge of atmospheric OH is imperfect and uncertain, and in future work we plan to explore alternate specifications of OH that are consistent with our knowledge of atmospheric chemistry and MCF trends and gradients.

4.6.2 Alternate optimizer and source signature uncertainty

Errors in the specification of the δ¹³C source signatures can have a significant impact on the inferred methane emissions (Thanwerdas et al., 2022). While we have explored alternate specifications, it is possible that the true uncertainty in δ¹³C source signatures is larger than the range we have explored. Optimizing the δ¹³C source signatures with a realistic prior covariance structure may yield larger but more realistic error bounds on source-specific methane emissions. We plan to explore that option in the future, which will require an alternate to the conjugate gradient optimizer (Lanczos, 1950) we currently use. We have tested the M1QN3 optimizer used by Thanwerdas et al. (2022) and have found its convergence to be slow and inefficient for our system. Therefore, we plan to explore and implement alternate optimizers that can work efficiently on nonlinear problems in order to have the option of estimating δ¹³C source signatures. Concurrently, we will work on a more complete characterization of the δ¹³C source signature uncertainty, which will be required in order to derive a prior error covariance matrix for δ¹³C.

4.6.3 OSSEs

We have tested the ability of existing δ¹³C observations to infer mechanisms behind the recent CH₄ growth and separate different CH₄ source types, and we found that the ability to distinguish fossil from microbial emissions – as reflected by the posterior correlation between them – is limited at policy-relevant scales (Sect. 3.5). We strongly suspect that this is a limitation of the existing δ¹³C observational coverage and not of the inversion technique. If we consider expanding the δ¹³C measurement network to improve that ability in the future, we need to quantify the added value of different expansion strategies. We plan to do this with observation system simulation experiments (OSSEs) simulating different observational networks, as we have done for ¹⁴C of CO₂ in the past (Basu et al., 2016).

4.6.4 Satellite CH₄ retrievals

Several satellites have been launched by various space agencies in the past decades to estimate atmospheric CH₄ from space, and several more are slated to go up over the next decade. As the technique to use δ¹³C in CH₄ inversions matures, we hope to eventually add satellite CH₄ data to such inversions to provide stronger regional constraints.

4.7 Future needs

While our effort to quantify source-specific methane emissions from atmospheric δ¹³C measurements has yielded several significant results, it has also uncovered areas for further progress. First, Fig. 3 shows a lack of δ¹³C data over the tropics and close to source regions. This is due to the historical choice of preferentially measuring δ¹³C at background sites and the general lack of tropical samples. More δ¹³C measurements in historically under-sampled regions would let us better separate tropical emissions and attribute continental emissions to specific sources. Second, while the soil sink has a small impact on CH₄-only studies, its impact on δ¹³C studies is large owing to the strong fractionation. Recent studies have suggested that the soil sink may be underestimated in current biogeochemical models (Oh et al., 2020) and may have trends due to a changing climate (Ni and Groffman, 2018; Murguia-Flores et al., 2021), which are both aspects that need better characterization. Third, our knowledge of the OH fractionation of ¹³CH₄ and the distribution of tropospheric Cl is a source of significant uncertainty in partitioning methane sources using δ¹³C data (Sect. 3.3). There is a wide range of estimates for OH fractionation in the literature (Rust and Stevens, 1980; Cantrell et al., 1990; Saueressig et al., 2001; Whitehill et al., 2017) and an analogously wide range of distributions of tropospheric Cl (Allan et al., 2007; Hossaini et al., 2016; Gromov et al., 2018). The uncertainty in both of these needs to be reduced in order to use δ¹³C data more effectively. Finally, while we have ruled out a sink-only explanation for the recent growth of methane (Lan et al., 2021), a small contribution from sink variations is possible. Since direct measurements of OH are few and far between and atmospheric levels of MCF are too low for quantifying recent variations in OH, we need alternate ways of separating source and sink influences. One possibility is to use measurements of δD in methane to constrain the sink (Quay et al., 1999; Sowers, 2006; Whiticar and Schaefer, 2007). We plan to explore the potential of δD measurements in future work.