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:

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:

ℋ^* 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 1The 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.

Download

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 1Emissions 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).

Download Print Version | Download XLSX

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 2Locations 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.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 2Nomenclature 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 %.

Download Print Version | Download XLSX

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 3Minimization 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.

Download

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 4Global 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.

Download

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 5RMSE 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.

Download

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 6REF 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 3Global 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.

Download Print Version | Download XLSX

Figure 7REF 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.