Variational inverse modeling within the Community Inversion Framework v1.1 to assimilate δ¹³C(CH₄) and CH₄: a case study with model LMDz-SACS

Atmospheric CH₄ mole fractions resumed their increase in 2007 after a plateau during the 1999–2006 period, indicating relative changes in the sources and sinks. Estimating sources by exploiting observations within an inverse modeling framework (top-down approaches) is a powerful approach. It is, nevertheless, challenging to efficiently differentiate co-located emission categories and sinks by using CH₄ observations alone. As a result, top-down approaches are limited when it comes to fully understanding CH₄ burden changes and attributing these changes to specific source variations. δ¹³C(CH₄)_source isotopic signatures of CH₄ sources differ between emission categories (biogenic, thermogenic, and pyrogenic) and can therefore be used to address this limitation. Here, a new 3-D variational inverse modeling framework designed to assimilate δ¹³C(CH₄) observations together with CH₄ observations is presented. This system is capable of optimizing both the emissions and the associated source signatures of multiple emission categories at the pixel scale. To our knowledge, this represents the first attempt to carry out variational inversion assimilating δ¹³C(CH₄) with a 3-D chemistry transport model (CTM) and to independently optimize isotopic source signatures of multiple emission categories. We present the technical implementation of joint CH₄ and δ¹³C(CH₄) constraints in a variational system and analyze how sensitive the system is to the setup controlling the optimization using the LMDz-SACS 3-D CTM. We find that assimilating δ¹³C(CH₄) observations and allowing the system to adjust isotopic source signatures provide relatively large differences in global flux estimates for wetlands (−5.7 Tg CH₄ yr⁻¹), agriculture and waste (−6.4 Tg CH₄ yr⁻¹), fossil fuels (+8.6 Tg CH₄ yr⁻¹) and biofuels–biomass burning (+3.2 Tg CH₄ yr⁻¹) categories compared to the results inferred without assimilating δ¹³C(CH₄) observations. More importantly, when assimilating both CH₄ and δ¹³C(CH₄) observations, but assuming that the source signatures are perfectly known, these differences increase by a factor of 3–4, strengthening the importance of having as accurate signature estimates as possible. Initial conditions, uncertainties in δ¹³C(CH₄) observations, or the number of optimized categories have a much smaller impact (less than 2 Tg CH₄ yr⁻¹).

1 Introduction

Methane (CH₄) is a powerful greenhouse gas and is responsible for 23 % (Etminan et al., 2016) of the radiative forcing induced by the well-mixed greenhouse gases (CO₂, CH₄, N₂O). Atmospheric CH₄ mole fractions have increased quasi-continuously since the preindustrial era and by about 9 ppb yr⁻¹ (parts per billion per year) from 1984 to 1998 (Dlugokencky, 2021). After a plateau between 1999 and 2006 that still generates attention and controversy (e.g., Fujita et al., 2020; Thompson et al., 2018; McNorton et al., 2018; Turner et al., 2017; Schaefer et al., 2016; Schwietzke et al., 2016; Rice et al., 2016), the mole fractions resumed their increase at a large rate, exceeding 10 ppb yr⁻¹ in 2014 and 2015. Trends in atmospheric CH₄ are caused by a small imbalance between large sources and sinks. Assessing their spatiotemporal characteristics is particularly challenging, considering the variety of CH₄ emissions. Yet, identifying and quantifying the processes contributing to these changes is mandatory to formulate relevant CH₄ mitigation policies that would contribute to meeting the target of the 2015 UN Paris Agreement on Climate Change and to limit climate warming to 2 ^∘C.

Thanks to continuous efforts of surface monitoring networks, the spatial coverage and the accuracy of the atmospheric CH₄ measurements provided to the scientific community have increased over the last decades. Consequently, top-down estimates using inversion methods emerged and became relevant, along with bottom-up estimates, for explaining and quantifying the recent sources and sinks variations. The first inverse modeling techniques were designed in the late 1980s and early 1990s for inferring greenhouse gas sources and sinks from atmospheric CO₂ measurements (Enting and Newsam, 1990; Newsam and Enting, 1988). Without a regularization of the problem, e.g., providing prior information, the inverse problem is ill-conditioned (or ill-posed). It means that there is no unique solution to the problem but also that a small error in the assimilated data (here atmospheric observations) can result in large errors in the derived solution. Several inversion methods have been designed over the years, among which there are analytical (e.g., Bousquet et al., 2006; Gurney et al., 2002), ensemble (e.g., Zupanski et al., 2007; Peters et al., 2005), and variational methods (e.g., Chevallier et al., 2005). The variational formulation uses the adjoint equations of a specific model to compute the gradient of a cost function and then minimize it, for example, using a gradient descent method. Computational times and memory costs do not scale with the number of measurements and the number of variables to control, contrary to the analytical and ensemble methods, which can hardly accommodate very large observational datasets and control vectors at the same time. Thus, the variational formulation is preferred to the others when optimizing emissions and sinks at the pixel scale using large volumes of observational data, although its main limitation is the numerical cost for accessing posterior uncertainties when there is nonlinearity in the inversion problem (Berchet et al., 2021).

Inversion systems generally assimilate measurements from ground-based stations and/or satellites to constrain the global sources and sinks of CH₄, starting from a prior knowledge of these. These systems are very effective for providing total emission estimates (e.g., Saunois et al., 2020; Bergamaschi et al., 2018, 2013; Saunois et al., 2017; Houweling et al., 2017, and references therein). However, differentiating the contributions of multiple co-located CH₄ source categories is challenging as it only relies on different seasonality cycles and on applied spatial distributions and error correlations (e.g., Bergamaschi et al., 2013, 2010). The atmospheric isotopic signal contains additional information on CH₄ emissions that can help to separate emission categories based on their source origin. The atmospheric isotopic signal δ¹³C(CH₄) is defined as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left({\mathrm{CH}}_{\mathrm{4}}\right)={\displaystyle \frac{R}{R_{\text{std}}}}-\mathrm{1},\end{array}$$

where R and R_std denote the sample and standard ¹³CH₄:¹²CH₄ ratios. We use the Vienna Pee Dee Belemnite (V-PDB) scale with R_std = 0.00112372 (Craig, 1957) throughout this paper. The isotopic source signatures of CH₄, here denoted by δ¹³C(CH₄)_source, notably differ between emission categories ranging from ¹³C-depleted biogenic sources (−61.7 ± 6.2 ‰; 1 standard deviation) and thermogenic sources (−44.8 ± 10.7 ‰) to ¹³C-enriched thermogenic sources (−26.2 ± 4.8 ‰) (Sherwood et al., 2017; Schwietzke et al., 2016), although the distributions are very large and overlaps exist between the extreme values. Consequently, δ¹³C(CH₄) depends on both CH₄ emissions and their isotopic signatures. Saunois et al. (2017) pointed out that many emission scenarios inferred from atmospheric inversions are not consistent with δ¹³C(CH₄) observations and that this constraint must be integrated into the inversion systems to avoid such inconsistencies. In addition, they highlighted the sensitivity of the atmospheric isotopic signal to the source partitioning and prescribed isotopic ratios. Since the 1990s, δ¹³C(CH₄) has been monitored at multiple sites, providing opportunities to use this constraint within an inversion framework. In addition, these values have been shifting towards more negative values since 2006 (Nisbet et al., 2019) when CH₄ trends resumed their increase, suggesting that this isotopic data can help to understand the processes that contributed to the renewed growth. However, implementing the assimilation of such measurements into an inversion system is not straightforward and introduces additional complexity.

Hereinafter, the assimilation of δ¹³C(CH₄) observations to constrain the estimates of an inversion is referred to as the isotopic constraint. The implementation of such a constraint in an inversion system has already been attempted in previous studies focusing on CH₄ (e.g., Thompson et al., 2018; McNorton et al., 2018; Rigby et al., 2017; Rice et al., 2016; Schaefer et al., 2016; Schwietzke et al., 2016; Rigby et al., 2012; Neef et al., 2010; Bousquet et al., 2006; Fletcher et al., 2004) but, to our knowledge, never in a variational system associated to a 3-D chemistry transport model (CTM). Adding this isotopic constraint to a variational inversion system is challenging as, in contrast to an analytic inversion in which the response functions of the model are precomputed, the isotopic constraints have to be considered both in the forward (simulated isotopic values) and the adjoint (sensitivity of isotopic observations to optimized variables) versions of the model.

This new system was implemented in the Community Inversion Framework (CIF), supported by the European Union H2020 project VERIFY (http://community-inversion.eu, last access: 3 June 2022) and required to implement new forward, tangent linear and adjoint operations. The forward operations were previously used to estimate the impact of the Cl sink on the modeling of CH₄ and δ¹³C(CH₄) in LMDz-SACS (Thanwerdas et al., 2019). The purpose of this study is to present the technical implementation of the isotopic constraint in a variational inversion system and to investigate the sensitivity of this new configuration to different parameters. Our aim is not to estimate trends in sectoral emissions over the last 2 decades: future studies will address the estimation of CH₄ emissions over longer periods of time using this new system. The technical implementation and the various tested configurations are presented in Sect. 2. We analyze the results in Sect. 3. Section 4 presents our conclusions and recommendations on using such a multi-constraint variational system.

2 Methods

2.1 Theory of variational inversion

The notations introduced here follow the convention defined by Ide et al. (1997) and Rayner et al. (2019). The observation vector is called y^o. It includes all available observations here, namely CH₄ and δ¹³C(CH₄) measurements retrieved by surface stations, over the full simulation time window (see Sect. 2.4.2). The associated errors are 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. x is the control vector and includes all the variables (here CH₄ surface fluxes, initial CH₄ mole fractions, source signatures δ¹³C(CH₄)_source, and initial δ¹³C(CH₄) values) optimized by the inversion system. Hereinafter, these variables will be referred to as the control variables. Prior information about the control variables are provided by the vector x^b. Its associated errors are also assumed to be unbiased and Gaussian and are described within the error covariance matrix B. ℋ is the observation operator that projects the control vector x into the observation space. This operator mainly consists of the 3-D CTM (here LMDz-SACS, as introduced in Sect. 2.2). Nevertheless, the CTM is followed by spatial and time operators, which interpolate the simulated fields to produce simulated equivalents of the assimilated observations at specific locations and times, making the simulations and observations comparable. An additional transformation operator, implemented in the new system, enables comparison between distinct simulated tracers, e.g., ¹²CH₄ and ¹³CH₄, and observations, e.g., δ¹³C(CH₄) (see Sect. 2.3).

In a variational formulation of the inference problem that allows for ℋ nonlinearity, the cost function J is defined as follows:

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

The cost function is therefore a sum of the following two parts:

1. The first part is induced by the differences between the posterior and prior variables (J_b).

2. The second is induced by the differences between simulations and observations (J_o).

The minimum of J can be reached iteratively with a descent algorithm that requires several computations of the gradient of J with respect to the control vector x, as follows:

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

ℋ^* denotes the adjoint operator of ℋ. As in analytical and ensemble methods, the variational formulation necessitates the inversion of both error matrices R and B. In most applications, R is considered diagonal as point observations are distant in time and space (i.e., uncorrelated observation errors), allowing for the inverse to be calculated easily, although that assumption should be revised with the increasing availability of satellite sources (Liu et al., 2020). B is rarely diagonal due to spatial and temporal correlations of errors in the fluxes. However, B is often decomposed as combinations of smaller matrices, e.g., using Kronecker products of sub-correlation matrices, which allows us to compute its inverse by blocks.

2.2 The chemistry transport model

The LMDz general circulation model (GCM) is the atmospheric component of the Institut Pierre-Simon Laplace Coupled Model (IPSL-CM) developed at the Laboratoire de Météorologie Dynamique (LMD) (Hourdin et al., 2006). The version of LMDz we use is an offline version dedicated to the inversion framework created by Chevallier et al. (2005). The precomputed air mass fluxes provided by the online version of LMDz are given as inputs to the transport model, significantly reducing the computational time. The model is set up 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 up to about 75 km. About 20 levels are dedicated to the stratosphere and the mesosphere. The model time step is 30 min, and the output mole fractions are 3 h snapshots. The horizontal winds are nudged towards ECMWF meteorological analyses (ERA-Interim) in the online version of the model and then fed to the offline version. Vertical diffusion is parameterized by a local approach from Louis (1979), and deep convection processes are parameterized by the Tiedtke (1989) scheme.

The offline model LMDz is coupled with the Simplified Atmospheric Chemistry System (SACS; Pison et al., 2009). This chemistry system was previously used to simulate the oxidation chain of hydrocarbons, including CH₄, formaldehyde (CH₂O), carbon monoxide (CO), and molecular hydrogen (H₂). For the purpose of this study, this system was converted into a chemistry parsing system. It follows the same principle as the one used by the regional model CHIMERE (Menut et al., 2013) and therefore allows for user-specific chemistry reactions. As a result, it generalizes the previous SACS module to any possible set of reactions. The adjoint code has also been implemented to allow variational inverse modeling. The different species are either prescribed (here OH, O(¹D), and Cl) or simulated (here ¹²CH₄ and ¹³CH₄). The prescribed species are neither transported in LMDz, nor are their mole fractions updated through chemical production or destruction. Such species are only used to calculate reaction rates to update the simulated species at each model time step. In this study, the isotopologues ¹²CH₄ and ¹³CH₄ are simulated as separate tracers, and CH₄ is defined as the sum of both isotopologues. Cl + CH₄ oxidation was implemented to complete the chemical removal of CH₄, which previously only accounted for OH + CH₄ and O(¹D) + CH₄ in the SACS scheme.

In the atmosphere, radicals (OH, O(¹D), or Cl) react faster with ¹²CH₄ than with ¹³CH₄. This effect is called the kinetic isotope effect (KIE) or the fractionation effect. Fractionation values are prescribed to the different sinks in SACS. Here, this value is defined by $\text{KIE}=k_{\mathrm{12}}/k_{\mathrm{13}}$, where k₁₂ is the rate constant of a reaction between a radical and ¹²CH₄. k₁₃ is the rate constant of the reaction between the same radical and ¹³CH₄. Additional information and prescribed KIE values are provided in Sect. S2 in the Supplement.

The chemistry transport LMDz-SACS is used to test the new variational inverse modeling system that is described in the next section.

2.3 Technical implementation of the isotopic constraint

The isotopic multi-constraint system was implemented in the CIF. The CIF has been designed to allow comparison of different approaches, models, and inversion systems used in the inversion community (Berchet et al., 2021). Different atmospheric transport models, regional and global and Eulerian and Lagrangian, are implemented within the CIF. The system presented in this paper was originally designed to run and be tested with LMDz-SACS but can theoretically be coupled with all models implemented in the CIF framework. The system is able to do the following:

• assimilate δ¹³C(CH₄) and CH₄ observations together,

• independently optimize fluxes and isotopic signatures for multiple emission categories, and

• optimize δ¹³C(CH₄) and CH₄ initial conditions.

Figure 1 The minimization iteration process in the newly designed system. The black circles with a gold border indicate the reading direction to follow. Step 1 (blue rectangle) refers to a forward run. Step 2 (orange rectangle) refers to the forward and adjoint operations required to compare observations and simulated values. Step 3 (green rectangle) refers to an adjoint run (this step must be read from the right to the left). Step 4 (red ellipse) refers to the minimization of the cost function operated by the dedicated minimization algorithm. Note that results of Step 2 are used in both the minimization process (red ellipse) and as inputs for Step 3. The minimization iteration process followed by the previous system is also illustrated in Fig. S1 in the Supplement.

Figure 1 shows the different steps of a minimization iteration of the cost function. Each iteration performed with the descent algorithm can be decomposed into the four main steps presented below. For clarity, here we only present the optimization of CH₄ fluxes and associated source signatures, but CH₄ and δ¹³C(CH₄) initial conditions can also be optimized by the system following the same process.

1. The process starts with a forward run. The different flux variables are extracted and converted into ¹²CH₄ and ¹³CH₄ mass fluxes for each category following Eqs. (5) and (6) below.

   $$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{12}}^i={\displaystyle \frac{M_{\mathrm{12}}}{M_{\mathrm{T}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+A^i}}\cdot F_{\mathrm{T}}^i\text{(6)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{13}}^i={\displaystyle \frac{M_{\mathrm{13}}}{M_{\mathrm{T}}}}\cdot {\displaystyle \frac{A^i}{\mathrm{1}+A^i}}\cdot F_{\mathrm{T}}^i,\end{array}$$

   with

   $$\begin{array}{}\text{(7)}& A^i=(\mathrm{1}+{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}({\mathrm{CH}}_{\mathrm{4}}{)}_{\mathrm{source}}^i)\cdot R_{\text{std}}.\end{array}$$

   $F_{\mathrm{T}}^i$, $F_{\mathrm{12}}^i$ and $F_{\mathrm{13}}^i$ are the CH₄, ¹²CH₄, and ¹³CH₄ mass fluxes of a specific category i, respectively. M_T, M₁₂, and M₁₃ are the CH₄, ¹²CH₄, and ¹³CH₄ molar masses, respectively. ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}({\mathrm{CH}}_{\mathrm{4}}{)}_{\mathrm{source}}^i$ is the isotopic signature of the category i. M_T should preferably depend on M₁₂ and M₁₃ when converting the mass fluxes, as follows:

   $$\begin{array}{}\text{(8)}& M_{\mathrm{T}}={\displaystyle \frac{M_{\mathrm{12}}+A^i\cdot M_{\mathrm{13}}}{\mathrm{1}+A^i}}.\end{array}$$

   However, the complexity of the forward, tangent linear, and adjoint codes would be largely enhanced by such a relationship. The code structure would also be less generic, i.e., it could not be used for a joint assimilation of multiple isotopologues of CH₄, such as both δ¹³C(CH₄) and δD(CH₄). We choose to implement M_T as a constant that can be prescribed freely by the user, therefore without considering any influence of the M₁₂ and M₁₃ values, which are also prescribed by the user. As the observed isotopic source signatures roughly vary between −70 ‰ and −10 ‰, a maximum variation of 0.004 % in M_T could be expected. It will very likely not affect the results of our study or that of any other inversion performed with our system.

   The ¹²CH₄ and ¹³CH₄ total fluxes are then calculated by summing all categories and are used by the model LMDz-SACS to simulate the ¹²CH₄ and ¹³CH₄ atmospheric mole fractions over the time window considered. After the simulation, the simulated values are converted to the CH₄- and δ¹³C(CH₄)-simulated equivalent of the assimilated observations using Eqs. (9) and (10) below, as follows:

   $$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\left[{\mathrm{CH}}_{\mathrm{4}}\right]=\left[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}\right]+\left[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\right]\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left({\mathrm{CH}}_{\mathrm{4}}\right)={\displaystyle \frac{\left[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\right]}{\left[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}\right]}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\text{std}}}}-\mathrm{1}.\end{array}$$

   [CH₄], [¹²CH₄], and [¹³CH₄] are CH₄, ¹²CH₄, and ¹³CH₄ atmospheric mole fractions simulated by the model in moles per mole (hereafter mol mol⁻¹), respectively.

2. These simulated values are then compared to the available observations in order to compute ℋ(x)−y^o, which is further used to infer the cost function and generate CH₄ and δ¹³C(CH₄) adjoint forcings (indicated by the superscripted asterisk) that compose the vector δy^* as follows:

   $$\begin{array}{}\text{(11)}& \mathit{\delta }{\mathit{y}}^{*}={\mathbf{R}}^{-\mathrm{1}}\left(\mathcal{H}\right(\mathit{x})-{\mathit{y}}^{\mathrm{o}}).\end{array}$$

   Although this vector is normally used directly as input to the adjoint model (see Eq. 4), the CH₄ and δ¹³C(CH₄) adjoint forcings must first be converted into the ¹²CH₄ and ¹³CH₄ adjoint forcings in the new system.

3. The newly designed adjoint code that converts CH₄ and δ¹³C(CH₄) adjoint forcings into ¹²CH₄ and ¹³CH₄ adjoint forcings is based on Eqs. (12)–(14), depending on the type of the initial observation.

   $$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathrm{CH}}_{\mathrm{4}}}^{*}=[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathrm{CH}}_{\mathrm{4}}}^{*}=[{\mathrm{CH}}_{\mathrm{4}}{]}^{*}\text{(13)}& {\displaystyle }& {\displaystyle }[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{*}=-{\displaystyle \frac{\left[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}\right]}{[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}{]}^{\mathrm{2}}}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\text{std}}}}\cdot {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left({\mathrm{CH}}_{\mathrm{4}}\right)}^{*}\text{(14)}& {\displaystyle }& {\displaystyle }[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{*}={\displaystyle \frac{\mathrm{1}}{\left[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}\right]}}\cdot {\displaystyle \frac{\mathrm{1}}{R_{\text{std}}}}\cdot {{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\left({\mathrm{CH}}_{\mathrm{4}}\right)}^{*}.\end{array}$$

   $[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathrm{CH}}_{\mathrm{4}}}^{*}$ and $[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathrm{CH}}_{\mathrm{4}}}^{*}$ are adjoint forcings associated with CH₄ observations. $[{}^{\mathrm{12}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{*}$ and $[{}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}{]}_{{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}}^{*}$ are adjoint forcings associated with δ¹³C(CH₄) observations. The adjoint code of the CTM is then run with these adjoint forcings as inputs.

   Outputs of the adjoint run provide the sensitivities of the adjoint forcings to the ¹²CH₄ and ¹³CH₄ mass fluxes of a specific category i, denoted by $F_{\mathrm{12}}^{*,i}$ and $F_{\mathrm{13}}^{*,i}$. Equations (15) and (16) convert them back to sensitivities to the initial control variables, denoted by $F_{\mathrm{T}}^{*,i}$ and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}({\mathrm{CH}}_{\mathrm{4}}{)}_{\mathrm{source}}^{*,\mathrm{i}}$.

   $$\begin{array}{}\text{(15)}& {\displaystyle }{F}_{\mathrm{T}}^{*,i}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+A}}\cdot [{\displaystyle \frac{M_{\mathrm{12}}}{M_{\mathrm{T}}}}\cdot F_{\mathrm{12}}^{*,i}+{\displaystyle \frac{M_{\mathrm{13}}}{M_{\mathrm{T}}}}\cdot A\cdot F_{\mathrm{13}}^{*,i}]\text{(16)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}({\mathrm{CH}}_{\mathrm{4}}{)}_{\mathrm{source}}^{*,i}=& R_{\text{std}}\cdot {\displaystyle \frac{F_{\mathrm{T}}}{(\mathrm{1}+A{)}^{\mathrm{2}}}}\cdot [{\displaystyle \frac{M_{\mathrm{13}}}{M_{\mathrm{T}}}}\cdot F_{\mathrm{13}}^{*,i}\\ & -{\displaystyle \frac{M_{\mathrm{12}}}{M_{\mathrm{T}}}}\cdot F_{\mathrm{12}}^{*,i}].\end{array}\end{array}$$

4. The minimization algorithm utilizes these sensitivities to compute the gradient of the cost function. It then finds an optimized control vector which reduces the cost function and is used for the next iteration.

In order to confirm that the several adjoint operations have been correctly implemented, we also provide the results of multiple adjoint tests in Sect. S4 in the Supplement.

Table 1 Emissions and flux-weighted isotopic signatures of the CH₄ sources averaged over 2012–2017 for different categories and their subcategories. Prior uncertainties in fluxes are set to 100 % for all categories and subcategories. Note: Unc. is the prior uncertainty in the isotopic signature prescribed to the category or the subcategory as a percentage of the signature. CH19 is Chang et al. (2019), GA18 is Ganesan et al. (2018), SH17 is Sherwood et al. (2017), WA16 is Warwick et al. (2016), ZA16 is Zazzeri et al. (2016), TO12 is Townsend-Small et al. (2012), KL10 is Klevenhusen et al. (2010), BO06 is Bousquet et al. (2006), BR01 is Bréas et al. (2001), SA01 is Sansone et al. (2001), CH00 is Chanton et al. (2000), HO00 is Holmes et al. (2000), CH99 is Chanton et al. (1999), BE98 is Bergamaschi et al. (1998), and LE93 is Levin et al. (1993).

2.4 Setup of the reference simulation

The reference configuration (REF) is a variational inversion that optimizes the CH₄ emission fluxes and δ¹³C(CH₄) isotopic source signatures of the following five different categories: biofuels–biomass burning (BB), agriculture and waste (AGW), fossil fuels (FFs), wetlands (WET), and other natural sources (NAT). CH₄ and δ¹³C(CH₄) initial conditions are also optimized. The assimilation time-window is the 2012–2017 period. The five categories originate from an aggregation of 10 subcategories (see Table 1) and are chosen to be as isotopically consistent as possible. Sinks are not optimized here.

2.4.1 Control vector x and B matrix

We adopt the CH₄ emissions compiled for inversions performed as part of the global methane budget (Saunois et al., 2020). Anthropogenic (including biofuels) and biomass burning emissions are based on the EDGARv4.3.2 database (Janssens-Maenhout et al., 2017) and the GFED4s databases (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 the EDGARv4.3.2 database, which ended 2012, until 2017. The natural source emissions are based on averaged literature values; see Poulter et al. (2017) for wetlands, Kirschke et al. (2013) for termites, Lambert and Schmidt (1993) for oceanic sources, and Etiope (2015) for geological (onshore) sources. Emissions from geological sources have been scaled down to 15 Tg CH₄ yr⁻¹ in the prior emissions adopted in Saunois et al. (2020). All prior fluxes are prescribed at a monthly resolution and at the spatial resolution of LMDz. Globally averaged emissions over the 2012–2017 period are listed in Table 1.

Prior estimates of isotopic source signatures are provided either at the pixel scale (for wetlands), at the regional scale based on Atmospheric Tracer Transport Model Intercomparison Project (TransCom) regions (Patra et al., 2011) or at the global scale. The wetlands signature map is taken from Ganesan et al. (2018). Livestock isotopic source signatures are taken from Chang et al. (2019) and aggregated into the 11-region map by selecting region-specific values. Livestock source signatures have been likely decreasing over time since the 1990s due to changes in the C₃/C₄ diet within the major livestock-producing countries, and therefore, annual values are prescribed. However, these estimates end in 2013, and we set the years 2014 to 2017 equal to the year 2013. Consequently, only the year 2012 has a different prescribed value from the other years. Coal and oil, gas, and industry (OGI) isotopic signature values are inferred from Sherwood et al. (2017) and Zazzeri et al. (2016) and aggregated into the same 11-region map. The EDGARv4.3.2 categories PRO_OIL and PRO_GAS (fugitive emissions during oil and gas exploitation) largely contribute (∼ 90 %) to the total of the oil, gas, and industry subcategory. Therefore, we chose to neglect the influence of other sub-subcategories (such as industry) on the isotopic signature of the category. As for the biofuels–biomass burning category, we use region-specific signatures over the 11 regions. A global signature value is prescribed for each of the other categories. Except for the livestock category, all prior signatures are set constant over time. To infer the δ¹³C(CH₄)_source map of a category based on the subcategories, the ¹²CH₄ and ¹³CH₄ fluxes for each emission subcategory within a category are derived based on Eqs. (5) and (6) and added up. The resulting fluxes are then converted back to a δ¹³C(CH₄)_source map representing the aggregated isotopic signature of the category. Additional information regarding the chosen isotopic signatures and their references is provided in Sect. S1 in the Supplement.

In total, three values per month (10 d, 10 d, and the rest) for fluxes and their associated isotopic signatures are included in the control variables. Although the time variations in isotopic signatures are poorly constrained in the literature, we choose to include the same number of variables for fluxes and isotopic signatures in order to illustrate the full capabilities of the system and have it ready when more isotopic constraints will appear.

The portion of the diagonal of B associated to prior CH₄ emission fluxes is filled in with the variances set to 100 % of the square of the maximum of emissions over the cell and its eight neighbors during each month. Off-diagonal terms of B (covariances) are based on correlation e-folding lengths (500 km over land and 1000 km over sea). The same method is applied for isotopic source signatures, although a specific percentage of uncertainties deduced from the global values of Sherwood et al. (2017) is used to infer each category's diagonal term (see Table 1). No temporal correlations are considered here. Finally, prior uncertainties in initial conditions are set to 10 % for CH₄ (∼ 180 ppb) and 3 % for δ¹³C(CH₄) (∼ 1.4 ‰).

Figure 2 Locations of CH₄ and δ¹³C(CH₄) surface stations. Affiliated networks are not displayed. More information can be found in Tables S3 and S4 in the Supplement.

2.4.2 Observation vector y and R matrix

CH₄ observations are taken from the data archived at the World Data Centre for Greenhouse Gases (WDCGG) of the World Meteorological Organization Global Atmospheric Watch (WMO-GAW) program. We selected 66 stations from 13 surface monitoring networks providing in situ measurements of CH₄ mole fractions. The stations are displayed in Fig. 2.

δ¹³C(CH₄) observations are taken from 18 surface stations from the Global Greenhouse Gas Reference Network (GGGRN), part of the NOAA Earth System Research Laboratory's (ESRL) Global Monitoring Laboratory (NOAA GML). Air samples were collected on an approximately weekly basis during the 2012–2017 period and analyzed by the Institute of Arctic and Alpine Research (INSTAAR) to provide δ¹³C(CH₄) isotope ratio measurements. The analytical uncertainty of the isotopic measurements, based on a surveillance cylinder, is 0.06 ‰. In this study, we focused on estimating monthly and annual flux variations rather than investigating daily or weekly variations. Prescribing error correlations in the R matrix (introduced in Sect. 2.1) can be used to ensure that the inversion preferentially constrains the components we are interested in (i.e., long-term trend and seasonal cycle). In order to keep the R matrix diagonal and to focus on monthly and annual variations of the signal, we chose to use δ¹³C(CH₄) observational data based on a curve fitting the original δ¹³C(CH₄) observations. The fitting curve is a function including three polynomial parameters (quadratic) and eight harmonic parameters, as in Masarie and Tans (1995). After the fitting, the pseudo-observations were sampled at the same time as the original observations. We also hypothesized that the convergence would be slightly faster if a smooth curve fitting the real observations was used instead of the real observations, which appeared to be false (see Sect. 3.1). One sensitivity inversion aims at estimating the error introduced by this simplification (simulation S2 in Table 2).

The R matrix for both CH₄ and δ¹³C(CH₄) is defined as diagonal, assuming that observation errors are not correlated, neither in space nor in time. This diagonal matrix can be decomposed into two parts, i.e., measurement and model error variances. Measurement errors account for instrumental errors, whereas model errors encompass transport and representativity errors induced by the model as follows:

$$\begin{array}{}\text{(17)}& \mathbf{R}={\mathbf{R}}_{\text{measurement}}+{\mathbf{R}}_{\text{model}}.\end{array}$$

Here, we use the provided observation errors to fill the R_measurement diagonal matrix. GlobalView Methane (CH₄; GLOBALVIEW-CH4, 2009) values are used to represent model errors and prescribe variances at each station for CH₄ mixing ratio measurements in order to fill the R_model diagonal matrix. This simple approach has been used previously in atmospheric inversions (Locatelli et al., 2015, 2013; Yver et al., 2011; Bousquet et al., 2006; Rodenbeck et al., 2003). Errors in GlobalView CH₄ are computed at each site as the root mean square error (RMSE) of the measurements on a smooth curve fitting them. As GlobalView CH₄ does not provide errors for GlobalView CH₄ measurements, the same method has been applied here. The RMSE of the measurements on a smooth curve, fitting them over the 2012–2017 period, is prescribed as the standard deviation for each site providing δ¹³C(CH₄) measurements. These errors range between 3–19 ppb for CH₄ observations and 0.11–0.20 ‰ for δ¹³C(CH₄) observations. Mean prescribed errors for each station are provided in Tables S3 and S4 in the Supplement.

2.4.3 Spin-up

Before starting the inversion, the model has been spun-up over 30 years, using constant emissions and recycling meteorology from the year 2012 in order to consider the long timescales for isotopic changes (Tans, 1997). At the end of the spin-up, δ¹³C(CH₄) values have been offset (+1.4 ‰) to fit the global mean δ¹³C(CH₄) in January 2012, and CH₄ mole fractions have been scaled to fit the global mean CH₄ mole fraction in January 2012. Due to the nonlinearity of transport and mixing, offsetting δ¹³C(CH₄) initial values in a forward run can generate errors. This impact is discussed later, using a configuration where δ¹³C(CH₄) initial conditions have not been offset (S1).

2.4.4 Sensitivity tests

A set of nine different configurations, including REF, has been designed to assess the impact of assimilating δ¹³C(CH₄) observations in addition to CH₄ observations and also to evaluate the sensitivity of the inversion results to the system's setup.

Multiple parameters have been tested throughout the following configurations:

1. NOISO has no isotopic constraint. Therefore, this configuration only simulates CH₄ and assimilates CH₄ observations.

2. S1 uses δ¹³C(CH₄) initial conditions that are not offset and are therefore directly taken from the spin-up.

3. S2 assimilates the real δ¹³C(CH₄) observations instead of the fitting curve data.

4. In S3, the δ¹³C(CH₄) model uncertainties are divided by a factor 2.

5. T1 uses 10 subcategories instead of 5 aggregated categories, increasing the degrees of freedom.

6. In theory, the system is capable of optimally adjusting two source signatures if the assimilated information is sufficient. For instance, the system can choose to shift one signature downward and another upward in a given pixel in order to improve the fitting in this specific pixel. The configuration T2 has been specifically designed to investigate whether the system would be able to retrieve a realistic distribution (similar to REF) starting from globally averaged signatures for each category.

7. In T3, the δ¹³C(CH₄) source signature uncertainties are set to a very low value (1 %) in order to prevent the system from optimizing them. In other words, all changes are put on CH₄ emissions.

8. Finally, T4 applies both changes from T2 and T3.

Table 2 Nomenclature and characteristics of the configurations. Details are provided in Sect. 2.4.4.

^* Prior uncertainties in initial δ¹³C(CH₄) conditions have been set to 10 %.

Table 2 summarizes the different configurations and the associated changes. The configurations have been grouped into two sets to facilitate the analysis of results. On the one hand, S-group configurations (REF + S1–S4) have setup variations that are not expected to largely influence the results compared to REF. On the other hand, T-group configurations (T1–T4) alter parameters that are very likely to impact the results.

Figure 3 Minimization of the cost function for all configurations. (a) Value of the cost function with respect to the number of iterations. (b) CH₄ contribution to J_o. (c) δ¹³C(CH₄) contribution to J_o. (d) RMSE associated to observed–simulated CH₄. (e) RMSE associated to observed–simulated δ¹³C(CH₄). For clarity reasons, S1 and S3 initial values are not displayed because they are much larger than those of REF.

3 Results

3.1 Minimization of the cost function

The minimization process is performed using the M1QN3 algorithm (Gilbert and Lemaréchal, 1989). One full simulation (forward + adjoint) with the isotopic constraint necessitates about 170 central processing unit (CPU) hours to run for 6 years, i.e., 2.4 CPU hours per month simulated. The computational burden is increased by a factor of 2 in comparison to an inversion without the isotopic constraint due to the doubling of simulated tracers (¹²CH₄ and ¹³CH₄). One full simulation is generally enough to complete one iteration of the minimization process, but two or three simulations are sometimes required by M1QN3. Therefore, the number of simulations is slightly larger than the number of iterations. Figure 3 displays the minimization process of the cost function for all configurations.

Except for S1 and T1, the inversions were stopped when the gradient norm reduction exceeded 96 % for the third consecutive iteration. The number of iterations are compared to investigate the sensitivity of the computational cost to the setup. In total, 32 iterations (37 simulations) for NOISO, 43 iterations (47 simulations) for REF, and about 50 iterations for the others were necessary. Consequently, although assimilating δ¹³C(CH₄) observations requires at least 11 additional iterations, the setup has little influence on the number of iterations if the same convergence criteria is used. Also, using curve-fitted data instead of real observations do not reduce the computational burden as we first speculated.

S1 and T1 inversions were extended until their cost function reached the same reduction as REF in order to estimate the additional computational burden required to reach similar results when initial conditions are not offset (S1) and the number of categories is increased (T1). There were 10 and 21 additional iterations necessary for T1 and S1, respectively. For T1, it shows that increasing the degrees of freedom also increases the computational burden. For S1, it highlights the benefits of offsetting δ¹³C(CH₄) initial conditions.

As we assume no correlation of errors in R, J_o (see Eq. 3) can be divided into CH₄ and δ¹³C(CH₄) contributions. Figure 3 shows that all configurations lead to a fast reduction of the δ¹³C(CH₄) contribution. During the first 10 iterations, it decreased from 50 %–90 % (depending on the configuration) to 10 %–20 %. Conversely, the CH₄ contribution increased from 10 %–50 % to 80 %–90 %. By adjusting the isotopic source signatures (all configurations besides T3–T4), the system was able to efficiently and rapidly reduce the discrepancies between simulated and observed δ¹³C(CH₄). As a result, the δ¹³C(CH₄) RMSE decreased very rapidly during the first 10 iterations, while the CH₄ RMSE decreased at a roughly constant rate. Consequently, the system is preferentially adjusting δ¹³C(CH₄) over CH₄ values to reduce the cost function, presumably because the ratio of RMSE to the prescribed observational error for δ¹³C(CH₄) is, on average, about twice as large as for CH₄. In other terms, it is simpler for the system to adjust δ¹³C(CH₄) before attempting to modify CH₄. The ratio of the number of δ¹³C(CH₄) observations to the number of CH₄ observations is not expected to play a significant role in the convergence process, although we did not study this influence rigorously. This ratio is only expected to affect the contribution of a component (δ¹³C(CH₄) or CH₄) to the total cost function.

The decrease rate associated with δ¹³C(CH₄) RMSE can be increased by reducing the model uncertainties prescribed to the δ¹³C(CH₄) observations. S3 is an example of such an adjustment, as the model uncertainties have been divided by two. With this configuration, the system requires five fewer iterations than REF to reach a similar δ¹³C(CH₄) RMSE reduction but seven additional iterations to reach a similar CH₄ RMSE reduction. T3 and T4 configurations constrain the isotopic signatures; thus, the reduction in the δ¹³C(CH₄) contribution necessitates 25 more iterations than REF to reach similar RMSE reduction. To summarize, the decrease rate associated with δ¹³C(CH₄) RMSE is highly dependent on the prescribed uncertainties in δ¹³C(CH₄) observations and the ability of the system to adjust source signatures.

Figure 4 Global monthly δ¹³C(CH₄) and CH₄ means between 2012 and 2017. The dashed black and solid blue lines in each panel denote the observed and prior simulated values (REF), respectively. The red and green ranges show the maximum and minimum values of the T group and S group, respectively. The thick and dashed purple line denotes the posterior REF values. Globally averaged values are computed using a method similar to Masarie and Tans (1995), where a function including three polynomial parameters (quadratic) and eight harmonic parameters is fitted to each time series at available sites. The final value is obtained by performing a latitude-band weighted average over the marine boundary layer (MBL) sites. The latitude band width was set at 30^∘. The posterior NOISO lines were not included because (1) the posterior NOISO global source signature is −54.1 ‰, and the line would therefore reach lower values than the REF PRIOR, affecting the visual clarity of the upper plot. (2) The posterior NOISO CH₄ values are extremely close to the REF values, and including it would also affect the clarity of the lower plot.

3.2 CH₄ and δ¹³C(CH₄) fitting

As expected, the assimilation process greatly improves the agreement between simulated and observed values for both CH₄ and δ¹³C(CH₄). Figure 4 shows the globally averaged time series of CH₄ and δ¹³C(CH₄).

CH₄ RMSE using prior estimates is 19.4 ppb and drops to 14.3 ± 0.2 ppb (1σ) on average over all the configurations using posterior estimates. Prior estimates capture the observed CH₄ well, and the improvement is therefore relatively small. In addition, all configuration results regarding CH₄ are very similar. In particular, NOISO is not performing much differently than the other configurations, indicating that the additional isotopic constraint does not affect the fitting to CH₄ observations.

Prior δ¹³C(CH₄) prescribed in REF are continuously decreasing from −47.2 ‰ to −48.2 ‰ and thus agrees very poorly (RMSE is 0.47 ‰) with observed values. This can be due to an underestimation (values that are too negative) of some isotopic source signatures, an underestimation of the KIE values associated with the various sinks, an underestimation of the various sinks intensities (mostly Cl and OH), and/or a poor prior estimation of the source partitioning, i.e., an underestimation of ¹³C-enriched sources (FFs or BB) or an overestimation of ¹³C-depleted sources (WET or AGW). The data assimilation process reconciles simulated and observed δ¹³C(CH₄) (RMSE is 0.086 ± 0.008 ‰) for all configurations, albeit that small differences, depending on the setup, emerge.

The S group provides a better match to δ¹³C(CH₄) observations than the T group (0.081 ± 0.003 ‰ versus 0.091 ± 0.007 ‰). The fit is very similar within the S group. In contrast, the spread in the T group is larger, with δ¹³C(CH₄) RMSE being equal to 0.093 ‰, 0.091 ‰, and 0.099 ‰, respectively, for T2, T3, and T4. These results suggest that giving more freedom to the system to adjust the isotopic signatures and providing region-specific estimates of prior source signatures instead of global values may be key elements for reaching a better agreement. Best results (i.e., smallest RMSE) are obtained with T1 (0.079 ‰). However, this configuration necessitates 10 additional iterations to reach better results than REF. Without these additional iterations, REF would have the best results (0.081 ‰).

Figure 5 RMSE distribution over the surface stations. Panels (a, c) show the prior CH₄ and δ¹³C(CH₄) RMSE, and panels (b, d) show the posterior RMSE. Panels (a, b) show δ¹³C(CH₄) RMSE, and panels (c, d) show CH₄ RMSE. For clarity reasons, S1 prior is not shown for δ¹³C(CH₄) because the associated prior misfit is much larger than that of the other configurations. The box-and-whisker plots are covering the whole range of data. In panel (c), all station labels are identical; therefore, most of them are removed to improve the clarity.

Figure 5 shows the RMSE distribution at all measurement sites for each configuration. All sites exhibit a RMSE reduction (from prior to posterior) for both CH₄ and δ¹³C(CH₄), except for Bukit Kototabang, Indonesia (BKT), with T3 and T4 configurations. Furthermore, BKT, Moody, USA (WKT), Ulaan Uul, Mongolia (UUM), Anmyeon-do, Republic of Korea (AMY), and Pondicherry, India (PON), exhibit a posterior CH₄ RMSE above 25 ppb, showing that CH₄ measurements retrieved at these stations are not properly reproduced by the model, despite the optimization. Prescribed observation errors are likely not the main cause because mean values for these stations are large (10–15 ppb) but not the largest among all the assimilated stations. It can also be due to transport error or misrepresentation of sources close to the sites. Addressing this misfit is beyond the scope of this study, although the configuration influences the results. BKT and UUM fitting are notably deteriorated with T3 and T4 configurations. For example, BKT appears to be influenced by biomass burning sources in South East Asia, which are strongly dependent on the configuration (see Sect. 3.3). Moreover, T3 provides the poorest δ¹³C(CH₄) fitting at AMY (0.24 ‰). Therefore, using global values for source signatures and preventing the system from optimizing them lead to poorer fitting. On the contrary, T1 improves the results, indicating that additional degrees of freedom can help to reconcile simulations with observations, especially in South East Asia and East Asia where these stations are located.

3.3 Global and regional emission increments

We are primarily interested in the additional information provided by the assimilation of δ¹³C(CH₄) data. Rather than discussing the regional and global CH₄ emissions and comparing these results to previous estimates, we investigate the differences between emissions inferred from configurations with and without the additional isotopic constraint. Long-term inversions will be run in the future with this system to provide more robust estimates of CH₄ emissions and compare them to the existing literature.

The inversion time window is the 2012–2017 period. However, flux and source signature estimates in the 2012–2013 and 2016–2017 periods are not interpreted, as the system appears to require a 2-year spin-up (2012–2013) and a 2-year spin-down (2016–2017), over which the inversion problem is not sufficiently constrained, and isotopic signatures vary widely over time. Therefore, only the 2014–2015 estimates are analyzed in Sects. 3.3 and 3.4. Figure S2 in the Supplement shows the time series of isotopic signatures and illustrates this choice. These long effects are certainly caused by the relatively long relaxation timescales of isotopic ratios in the atmosphere (Tans, 1997) compared to that of total CH₄. Fully understanding this would require a lot of time and running multiple inversions (or possibly only tangent linear simulations), starting from different initial conditions spanning the prescribed uncertainty envelope, to infer until when the initial atmospheric isotopic ratios and/or isotopic source signatures can influence the time series of atmospheric isotopic ratios. This was too much work for this study but will certainly be addressed in future studies.

Figure 6 REF and NOISO emission increments for the 2014–2015 period. Prior estimates (PRIOR) are identical for both configurations. The color-filled bars show the differences between REF posterior and prior estimates (REF increment). The hatched bars show the differences between NOISO posterior and prior estimates (NOISO increment). The upper panel refers to the global emissions. The lower panels refer to multiple regions of the globe. The regions are shown on the lower right panel. Red and blue error bars represent the minimum and maximum of the T group and S group, respectively. Circles on the red error bar show the results from the T group.

Figure 6 shows global and regional increments from the NOISO and REF inversions relative to prior estimates. Hereinafter, these differences will be referred to as REF increment (REF – PRIOR) and NOISO increment (NOISO – PRIOR). The difference between both increments will be called an increment difference. Note that prior emissions are identical for all configurations. Posterior total emissions is 594.6 ± 1.2 Tg CH₄ yr⁻¹ over all configurations, indicating that the isotopic constraint and setup configurations do not significantly affect the posterior global emissions. A higher discrepancy between the budgets would have indicated a malfunction in the system, as the prescribed sinks are identical. The small associated standard deviation is likely caused by a slight difference in the fitting to the observations and/or by the spatial variability in the prescribed sink coupled with a small relocation of emissions, depending on the configuration. For instance, OH concentrations are larger in the tropics, and a relocation of emissions from the tropics to higher latitudes would be compensated for by larger global emissions. Between REF and NOISO, there is only a difference of 0.02 Tg CH₄ yr⁻¹. We can therefore conclude that the additional isotopic constraint either relocates the emissions or reallocates them between categories, as intended. All but one of the emission categories exhibit large changes between NOISO and REF, namely WET, FFs, AGW, and BB categories.

Overall, increments are large in regions with high emissions. Global increment differences (between REF and NOISO) in AGW (−6.4 Tg CH₄ yr⁻¹) and FF emissions (+8.6 Tg CH₄ yr⁻¹) are mainly due to regional increment differences in China and temperate Asia. AGW regional increment differences are equal to −2.1 Tg CH₄ yr⁻¹ in temperate Asia and in China. Similarly, FF regional increment differences are equal to +1.5 Tg CH₄ yr⁻¹ in temperate Asia and +5.0 Tg CH₄ yr⁻¹ in China. The WET global increment difference (−5.7 Tg CH₄ yr⁻¹) is mainly due to differences in Canada (−2.0 Tg CH₄ yr⁻¹) and South America (−2.3 Tg CH₄ yr⁻¹), but other regions such as Russia, temperate Asia, and South East Asia are involved. BB emissions are also modified when implementing the isotopic constraint. Their global increment difference is equal to +3.2 Tg CH₄ yr⁻¹, principally owing to regional increment differences in South East Asia (+1.7 Tg CH₄ yr⁻¹), Canada (+0.4 Tg CH₄ yr⁻¹), and Africa (+0.4 Tg CH₄ yr⁻¹). The NAT category exhibits very few changes (less than 1 Tg CH₄ yr⁻¹), even in relative values.

S group configurations infer posterior results that are consistent with REF, with only small variations, depending on the category and the region (see Table S5 in the Supplement). In particular, S1 provides roughly the same results as REF, but with more iterations, highlighting again that offsetting the initial conditions can help to reduce the computational burden without affecting the results. On the contrary, T-group configurations are affecting the increments, although T1 and T2 configurations are generally much closer to REF than T3 and T4. T1 (yellow dot; Fig. 6) and T2 (blue dot; Fig. 6) exhibit differences with the S group mainly in China, where WET and FF increments are modified (∼ −3 Tg CH₄ yr⁻¹). More importantly, almost freezing the isotopic signatures to their prior values (T3 and T4) results in increment differences 3 to 4 times larger than with REF, i.e., more than 10 Tg CH₄ yr⁻¹ at the global scale. It highlights the dependence of the inferred CH₄ emissions on the prior source signatures estimates. In other words, the quality of isotopic source signatures (values and uncertainties) appears to be critical for the robustness of emissions estimates.

Table 3 Global posterior CH₄ emissions by source category and region (Tg CH₄ yr⁻¹) for the REF configuration. Uncertainties are reported as the min–max range of all configurations in square brackets.

Figure 7 REF flux-weighted source signature increments for the 2014–2015 period. The color-filled bars show the differences between REF posterior and prior estimates (REF increment). The upper panel refers to the global emissions. The lower panels refer to multiple regions of the globe. The regions are shown on the lower right panel. Red and blue error bars represent the minimum and maximum of the T group and S group, respectively. Circles on error bars show the results from the T group.

3.4 Global and regional source signature increments

Isotopic source signatures are also optimized by the system. Figure 7 provides the differences of flux-weighted source signatures between REF posterior and prior estimates for different regions and each emission category.

With configurations that allow the source signatures to be optimized, all source signatures are shifted upwards by the inversions in order to correct the excessively strong negative trend in δ¹³C(CH₄). At the global scale, the flux-weighted source signatures of WET, FFs, AGW, BB, and NAT are increased by 1.7, 0.5, 0.9, 0.5, and 0.1 ‰, respectively. The global source signature is increased from −53.9 ‰ (prior) to −52.6 ± 0.2 ‰ (posterior, with a standard deviation over the configurations). More information is provided in Table S6 in the Supplement. The posterior global signature is strongly dependent on the KIE of atmospheric oxidation. This effect tends to deplete air in ¹³CH₄, shifting the δ¹³C(CH₄) to more positive values as the CH₄ molecules emitted by the sources are removed from the atmosphere. The mean KIE in our simulations depends on (1) the prescribed OH, O(¹D), and Cl concentrations and (2) the prescribed KIE values associated to the individual sinks. As the mean KIE is the same for all configurations, the posterior global source signatures are very close.

The WET global source signature, associated with REF posterior estimates, exhibits the larger upward shift compared to prior estimates, from a value of −60.8 ‰ to −59.1 ‰. Large upward WET source signature shifts are located in boreal regions (North America and Russia) but also in South America and temperate Asia. The AGW source signature is increased by 0.9 ‰, mainly due to changes in Asia. The FF source signature is increased by 0.5 ‰ due to a large increment in China (+1.2 ‰). Finally, the BB source signature is modified in South East Asia (+1.4 ‰) and Canada (+0.8 ‰).

These changes are consistent within the S group (see blue error bars in Fig. 7), although small variations are visible (e.g., ± 0.3 ‰ for WET in Canada). The source signature is therefore modified nearly to the same extent in all regions, no matter which configuration in the S group is analyzed. T1 (see the yellow dot in Fig. 7), with more optimized categories than the others, shows small differences at the global scale (less than 0.3 ‰ for all categories), although differences of more than 1 ‰ are found in China. Therefore, increasing the number of degrees of freedom leads to similar flux estimates but can affect the signatures at a regional scale. T2 estimates are shifted upward to reach a less negative global isotopic source signature without coming closer to the regional distribution of the S group. This is likely caused by the scarcity of δ¹³C(CH₄) stations, and correcting this behavior seems challenging without additional observations. The problem might be circumvented by using the regional scale rather than the pixel scale to optimize isotopic signature values. Future inversions will test this assumption.

These results must be interpreted with caution because the input data suffer from high uncertainties. The artificial increase in the source signatures by our system can hardly be related to the literature and former investigations. Consequently, it is challenging to conclude whether an increase in the source signatures would be more realistic (i.e., supported by observational data) than, for instance, only increasing the emissions of ¹³C-enriched sources such as BB. This system is only based on a mathematical and physical framework connecting the several groups of uncertainties (observational, prior fluxes, prior source signatures, and prior sinks) and finding the most likely solution. Better estimates of these uncertainties must be prescribed before obtaining robust results. In particular, the uncertainties on KIE values and sink intensities have not been tested here and could largely influence the results. Also, the uncertainties on source signatures are relatively smooth in REF compared to recent country-specific estimates (Sherwood et al., 2017). Assessing these uncertainties should be a key aspect for future studies using this new inversion system to quantify the global CH₄ budget.

3.5 Posterior uncertainties

Formally, posterior uncertainties are given by the Hessian of the cost function. This matrix can hardly be computed 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, or computing multiple configurations covering a given range of possibilities. Here, using multiple configurations provides insight into the posterior uncertainty associated with the posterior fluxes. We calculated the full uncertainty range using the minimum and maximum values among all the configurations, as in Saunois et al. (2020). WET, AGW, FF, and BB flux estimates (Table 3) exhibit an uncertainty of 10 %, 7 %, 19 %, and 38 %, respectively. BB is the most uncertain estimate relative to its intensity, although FF shows the largest absolute uncertainty (23 Tg CH₄ yr⁻¹). These uncertainties are unlikely to be affected by the assimilation of additional δ¹³C(CH₄) data because we expect the uncertainties on the isotopic source signatures to have a much larger influence. However, this remains to be tested in future work if posterior uncertainties can be calculated.

At present, M1QN3 is not the only optimization algorithm that can be utilized to perform variational inversions in the CIF. The CONGRAD algorithm (Fisher, 1998), which follows a conjugate gradient method combined with a Lanczos algorithm, is also implemented. In particular, it facilitates the computation of posterior uncertainties considerably. Any change in algorithm is very easy and accessible to any CTM embedded in the CIF. However, CONGRAD has not yet been tested with δ¹³C(CH₄) data. As CONGRAD is only designed for linear problems, using this algorithm could radically change the results of the inversions performed with the isotopic constraints, and future work will focus on using CONGRAD to perform the inversions with isotopic constraints.

4 Conclusions and perspectives

We present here a new variational inversion system designed to assimilate the observations of both a specific trace gas and its isotopic data. This system allows us to optimize both tracer emissions and associated isotopic signatures for multiple source categories. To test this system, we have assimilated CH₄ and δ¹³C(CH₄) data retrieved at different measurement sites over the globe.

Different configurations have been tested in order to assess the sensitivity of the system to the setup. We have shown that offsetting the δ¹³C(CH₄) initial conditions before the inversion (S1), using δ¹³C(CH₄) curve fitting data instead of the original observations (S2) and reducing the prescribed uncertainties in the δ¹³C(CH₄) observations (S3), has very little effect on the inferred fluxes (less than 2 Tg CH₄ yr⁻¹ for each category at the global scale). However, offsetting δ¹³C(CH₄) initial conditions before the inversion results in a reduced computational time (21 less iterations).

Other setup choices have more influence on the results. Increasing the number of source categories (T1) requires more computational time (10 more iterations) to reach a cost function (and RMSE) reduction similar to REF. Moreover, although the global posterior emissions with an increased number of categories are very close to those inferred with REF (less than 1 Tg CH₄ yr⁻¹), the posterior isotopic signatures can be modified in some regions (more than 1 ‰ in China). Also, starting from globally averaged values for the source signatures (T2) makes the system unable to retrieve the region-specific isotopic signatures from REF. Increasing the number of δ¹³C(CH₄) observations could help us to cope with this issue. Finally, configurations constraining the source signatures (T3–T4) show differences in global flux estimates of more than 10 Tg CH₄ yr⁻¹, compared to REF. This emphasizes the need for good prior source signature estimates.

The major drawback of this inversion system is undoubtedly the large computational burden of a full minimization process. At least 40 iterations appear to be necessary to reach a satisfying convergence state at the regional scale. For the LMDz-SACS model, a maximum of eight CPUs can be run in parallel, resulting in an elapsed time of 5–6 weeks to run one of the inversions of this study. A new generation of transport models, such as DYNAMICO (Dubos et al., 2015), could help to address this problem in the future by allowing more processors to run in parallel. Also, further developments will implement some parallelization methods to enable computational burden reduction (e.g., Chevallier, 2013). In addition, variational inversions, as implemented in the CIF, do not provide a quantification (even approximated) of the posterior uncertainties. Dedicated efforts need to be made to address this issue in the future and at an achievable numerical cost. In particular, using the CONGRAD algorithm instead of M1QN3 could be a solution, as both algorithms can be easily selected in the CIF. However, additional work is needed to ensure that switching the optimization algorithm does not affect the results inferred with our new system.

This system is implemented within the CIF framework and can therefore be used for inversions with the various CTMs embedded in the CIF, provided the adjoint codes of the models exist. As the operations developed for the purpose of this study are performed outside the model structure, forward, tangent linear, and adjoint codes from other CTMs do not require any modifications as long as the model is capable of simulating both ¹²CH₄ and ¹³CH₄ simultaneously. The prior input must be adapted to the new model (spatial and time resolution), but the format of the observational data and of the prescribed errors can be preserved. Also, due to the variational method benefits, the efforts dedicated to the preparation of inputs do not scale with either the size of the observational datasets or the length of the simulation time window. Therefore, this system is very powerful and is particularly relevant to study, in a consistent way, the influence of multiple physical parameters on atmospheric isotopic ratios such as the transport, the isotopic signatures, the emission scenarios, the KIE values, etc. We did not try to assess the sensitivity of the system to these parameters here as only the technical aspects of the system were tested. This will be part of future analyses.

As mentioned in the introduction, future work will address the estimation of CH₄ emissions over longer periods of time using this new system. For instance, the 2000–2006 CH₄ stabilization period and the subsequent renewed growth are particularly interesting to study using the isotopic constraint as global δ¹³C(CH₄) started to decrease after 2006. These periods of time have already attracted considerable critical attention from many inversion studies (with or without the isotopic constraint) and comparing the results derived from such a complete 3-D variational inversion system with other recent estimates should be highly relevant. The most important limitation of assimilating δ¹³C(CH₄) lies in the fact that very limited δ¹³C(CH₄) data are available, and therefore, evaluating the posterior simulated δ¹³C(CH₄) is often challenging, if not impossible. However, satellite and balloon/AirCore data can easily be used to evaluate the posterior simulated CH₄.

δ¹³C(CH₄) is not the only kind of isotopic data that can be assimilated in such a system. Many δD(CH₄) observations have also been retrieved during the 2004–2010 period at many different locations. These isotopic values can provide additional information that can further help to discriminate the co-emitted CH₄ fluxes (Rigby et al., 2012). Moreover, ethane (C₂H₆) is co-emitted with CH₄ by fossil fuel extraction and distribution (Kort et al., 2016; Smith et al., 2015), and observations have been available at a multitude of sites since the early 1980s. Therefore, assimilating these data can provide additional constraints. The system will therefore be improved in the future in order to assimilate δ¹³C(CH₄), δD(CH₄) and C₂H₆ observations together.