Variational inverse modelling within the Community Inversion Framework to assimilate δC(CH4) and CH4 : a case study with model LMDz-SACS

Framework to assimilate δC(CH4) and CH4 : a case study with model LMDz-SACS Joël Thanwerdas1,*, Marielle Saunois1, Antoine Berchet1, Isabelle Pison1, Bruce H. Vaughn2, Sylvia Englund Michel2, and Philippe Bousquet1 1Laboratoire des Sciences du Climat et de l’Environnement, CEA-CNRS-UVSQ, IPSL, Gif-sur-Yvette, France. 2INSTAAR University of Colorado, Boulder, CO, United States Correspondence: J. Thanwerdas (joel.thanwerdas@lsce.ipsl.fr)

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 mixing ratios resumed their increase at a large rate, exceeding 10 ppb/yr in 2014 and 2015. Trends in atmospheric CH 4 are caused by a small imbalance between large sources and sinks. Assessing their spatio-temporal characteristics is particularly challenging considering the variety of methane emissions. Yet, identifying and quantifying the processes contributing to these changes 5 is mandatory to formulate relevant CH 4 mitigation policies that would contribute to meet 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 methane measurements provided to the scientific community increased over the last decades. Consequently, top-down estimates using inversion methods emerged and became relevant, along with bottom-up estimates, to explain and quantify the 10 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 2 measurements (Enting and Newsam, 1990;Newsam and Enting, 1988). The inverse problem is considered as "ill-posed" (non-uniqueness of the solution, no continuity with the data) and therefore necessitates as many constraints as possible to be regularized. Several methods have been designed over the years, among which analytical (e.g., Bousquet et al., 2006;Gurney et al., 2002), ensemble (e.g., Zupanski et al., 2007;Peters 15 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 20 the pixel scale using large volumes of observational data, although its main limitation is the numerical cost to access posterior uncertainties.
Inversion systems generally assimilate measurements from ground-based stations and/or satellites to constrain the global sources and sinks of CH 4 , starting from a prior knowledge of these. These systems are very effective to provide total emission estimates (e.g., Saunois et al., 2020;Bergamaschi et al., 2018Bergamaschi et al., , 2013Saunois et al., 2017;Houweling et al., 2017, and references 25 therein). However, differentiating the contributions of multiple co-located CH 4 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., 2013Bergamaschi et al., , 2010. The atmospheric isotopic signal contains additional information on methane emissions that can help to separate emission categories based on their source origin. The atmospheric isotopic signal δ 13 C(CH 4 ) is defined as: 30 where R and R std denote the sample and standard 13 CH 4 : 12 CH 4 ratios. We use the VPDB scale with R std = 0.00112372 (Craig, 1957) throughout this paper. CH 4 source isotopic signatures δ 13 C(CH 4 ) source notably differ between emission categories ranging from 13 C-depleted biogenic sources (approx. -62 ‰) and thermogenic sources (approx. -44 ‰) to 13 C-enriched thermogenic sources (approx. -22 ‰) (Sherwood et al., 2017;Schwietzke et al., 2016). Consequently, δ 13 C(CH 4 ) depends on both the CH 4 emissions and their isotopic signatures. Saunois et al. (2017) pointed out that many emission scenarios inferred from atmospheric inversions are not consistent with δ 13 C(CH 4 ) observations and that this constrain must be integrated into the inversion systems to avoid such inconsistencies. Since the 1990s, δ 13 C(CH 4 ) have been monitored at multiple sites, although less than for total CH 4 , providing opportunities to use this constraint within an inversion framework. In addition, these values have been shifting towards smaller values since 2006 (Nisbet et al., 2019) when CH 4 trends resumed their increase, suggest- 5 ing that this isotopic data can help to understand the processes that contributed to the regrowth. However, implementing the assimilation of such measurements into an inversion system is not straightforward and introduces additional complexity.
Hereinfafter, the assimilation of δ 13 C(CH 4 ) 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 have already been attempted in previous studies focusing on CH 4 (e.g., Thompson et al., 2018;McNorton et al., 2018;Rigby et al., 2017;Rice et al., 2016;Schaefer 10 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. 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.

15
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 two decades : future studies will address the estimation of CH 4 sources 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. Sect. 4 presents our conclusions and recommendations on using such a multi-constraint variational 20 system.

Theory of variational inversion
The notations introduced here follow the convention defined by Ide et al. (1997). The observation vector is called y o . It includes here all available observations, namely CH 4 and δ 13 C(CH 4 ) measurements retrieved by surface stations, over the full simulation 25 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 4 surface fluxes, initial CH 4 mixing ratios, source signatures δ 13 C(CH 4 ) source and initial δ 13 C(CH 4 ) 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 30 also assumed to be unbiased and Gaussian and are described within the error covariance matrix B. H is the observation operator that projects the control vector x into the observation space. This operator mainly consists of the 3-D Chemistry-Transport Model (CTM) (here LMDz-SACS 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. 12 CH 4 and 13 CH 4 , and observations, e.g. δ 13 C(CH 4 ) (see In a variational formulation of the inference problem that allows for H non-linearity, the cost function J is defined as : The cost function is therefore a sum of two parts : -The first part is induced by the differences between the posterior and prior variables (J b ).
-The second is induced by the differences between simulations and observations (J o )

10
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: H * denotes the adjoint operator of H. Although the variational method is a powerful approach for dealing with large numbers of observations and control variables (several hundred thousands), it implies the inversion of both error matrices R and B. In 15 most applications, R is considered diagonal as point observations are distant in time and space, allowing inversing it easily, although that assumption may change 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 to compute its inverse by blocks.  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 4 , formaldehyde (CH 2 O), carbon monoxide (CO) and molecular hydrogen (H 2 ) together with methyl chloroform (MCF). For the purpose of this study, this system has been 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 5 previous SACS module to any possible set of reactions. The adjoint code has also been implemented to allow variational inverse Here, KIE is defined by KIE = k 12 /k 13 where k 12 is the constant rate of the reaction involving 12 CH 4 and k 13 is the constant rate of the same reaction involving 13 CH 4 . Additional information is provided in the supplement (Text S2).

The Chemistry-Transport Model
The chemistry-transport LMDz-SACS is used to test the new variational inverse modelling system that is described in the next 15 section.

Technical implementation of the isotopic constraint
The isotopic multi-constrain system was implemented in the Community Inversion Framework (CIF), supported by the European Union H2020 project VERIFY (http://www.community-inversion.eu). The CIF has been designed to allow comparison of different approaches, models and inversion systems used in the inversion community (Berchet et al., 2020). Different atmo-20 spheric transport models, regional and global, Eulerian and Lagrangian are implemented within the CIF. The system presented in this paper has been 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 : -Assimilate δ 13 C(CH 4 ) and CH 4 observations together.
-Independently optimize fluxes and isotopic signatures for multiple emission categories.

25
-Optimize δ 13 C(CH 4 ) and CH 4 initial conditions. 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 four main steps presented below. For clarity, we only present here the optimization of CH 4 fluxes and associated source signatures but CH 4 and δ 13 C(CH 4 ) 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 12 CH 4 and 13 CH 4 mass fluxes for each category following the Eq. (5)-(7) below.
T OT , F i 12 and F i 13 are the CH 4 , 12 CH 4 and 13 CH 4 fluxes in kg m −2 s −1 of a specific category i, respectively. M T OT , M 12 and M 13 are the CH 4 , 12 CH 4 and 13 CH 4 molar masses, respectively. δ 13 C(CH 4 ) i source is the source isotopic signature of the category i.
The 12 CH 4 and 13 CH 4 fluxes are then provided by summing all categories and used by the model LMDz-SACS to simulate the 12 CH 4 and 13 CH 4 atmospheric mixing ratios over the time-window considered. Finally, the simulated 10 values are converted back to CH 4 and δ 13 C(CH 4 ) simulated equivalent of the assimilated observations using Eq. (8) and (9) below : 2. These simulated values are then compared to the available observations in order to compute H(x) − y o which is further used to infer the cost function and generate CH 4 and δ 13 C(CH 4 ) adjoint forcings (indicated by the "*" star superscript symbol) that compose the vector δy * : 20 This vector is normally used directly as input to the adjoint model (see Eq. 4) but in the new system, the adjoint forcings CH 4 and 13 C(CH 4 ) must first be converted into the adjoint forcings 12 CH 4 and 13 CH 4 .
3. The newly designed adjoint code that converts CH 4 and δ 13 C(CH 4 ) adjoint forcings into 12 CH 4 and 13 CH 4 adjoint forcings is based on the Eq. (11)-(13) depending on the type of the initial observation. are adjoint forcings associated with δ 13 C(CH 4 ) 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 12 CH 4 and 13 CH 4 fluxes of a specific category i denoted F * ,i 12 and F * ,i 13 . Equations (14) and (15) convert them back to sensitivities to the initial control 5 variables, denoted F * ,i T OT and δ 13 C(CH 4 ) * ,i source .
4. The minimization algorithm uses these sensitivities to compute the gradient of the cost function. It then finds an optimized control vector that reduces the cost function and that is used for the next iteration.

Setup of the reference simulation
The reference configuration (REF) is a variational inversion that optimizes the CH 4 emission fluxes and δ 13 C(CH 4 ) source isotopic signatures of five different categories (biofuels-biomass burning, microbial, fossil fuels, natural and wetlands) and the CH 4 /δ 13 C(CH 4 ) initial conditions. The assimilation time-window is the period 2012-2017. The five categories originate from an aggregation of ten sub-categories (Table 1) and chosen to be as isotopically consistent as possible. Sinks are not optimized 15 here.

Control vector x and B matrix
We adopt the CH 4 emissions compiled for inversions performed as part of the Global Methane Budget (Saunois et al., 2020).

25
Source isotopic signatures are provided either at the pixel scale (for wetlands), at the regional scale based on TransCom regions (Patra et al., 2011)   H(x) δy* = R -1 (H(x) -y 0 ) Figure 1. The minimization iteration process in the newly designed system. The "step" black circles with gold border indicates 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 both in the minimization process (red ellipse) and as inputs for Step 3. The minimization iteration process followed by the previous system is also illustrated in the supplement (Fig. S1)  infer each category diagonal term (see Table 1). No temporal correlations are considered here. Finally, prior uncertainties on initial conditions are set to 10 % for CH 4 (∼ 180 ppb) and 3 % for δ 13 C(CH 4 ) (∼ 1.4 ‰). CH4 station 13 C-CH4 station  (Table S3 and S4).
in-situ measurements of CH 4 mixing ratios. The stations are displayed in Fig. 2. Table S3 in the supplement provides a list of these 66 stations and specific information.  Table S4 in the supplement provides a list of these 18 stations and specific information. The observed high-frequency temporal variability cannot be adequately reproduced by the LMDz-SACS model. Therefore, instead of assimilating the real observations, we used a smooth curve fitting the real observations. The fitting curve is a function including 3 polynomial parameters (quadratic) and 8 harmonic parameters. One sensitivity inversion aims at estimating the error introduced by this simplification (simulation S2 10 in Table 2).
The R matrix introduced in Sect. 2.1 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 : measurement and model error variance. Measurement errors account for instrumental errors while model errors encompass transport and representativity errors induced by the model : Here, we use the provided observation errors to fill the R measurement diagonal matrix. Globalview-CH 4 (Globalview-CH 4 , 2009) values are used to represent model errors and prescribe variances at each station for CH 4 mixing ratio measurements in order to fill the R model diagonal matrix. This simple approach has been used previously in atmospheric inversions (Locatelli 5 et al., 2015(Locatelli 5 et al., , 2013Yver et al., 2011;Bousquet et al., 2006;Rodenbeck et al., 2003). Errors in Globalview-CH 4 are computed at each site as the Root-Mean-Square-Error (RMSE) of the measurements on a smooth curve fitting them. As Globalview-CH 4 does not provide errors for δ 13 C(CH 4 ) measurements, the same method has been applied here. RMSE of the measurements on a smooth curve fitting them over the period 2012-2017 is prescribed as the standard deviation for each site providing δ 13 C(CH 4 ) measurements. 10

Spin-up
The model has been spun-up during 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, δ 13 C(CH 4 ) values have been offset to fit the δ 13 C(CH 4 ) global-mean in January 2012 and CH 4 mixing ratios have been scaled to fit the CH 4 global-mean mixing ratios in January 2012. Due to the non-linearity of transport and mixing, offsetting δ 13 C(CH 4 ) initial values in a forward run 15 can generate errors. This impact is discussed later using a configuration where δ 13 C(CH 4 ) initial conditions have not been offset (S1). Table 2. Nomenclature and characteristics of the configurations. Details are provided in Sect. 2.5. ** Prior uncertainties on initial δ 13 C(CH4)

Sensitivity tests
conditions have been set to 10 %.
Configuration Including REF, a set of 9 different configurations has been designed to assess the impact of assimilating δ 13 C(CH 4 ) observations in addition to CH 4 observations and also to evaluate the sensitivity of the inversion results to the system's setup.
Multiple parameters have been tested throughout the various configurations : NOISO has no isotopic constraint. Therefore, this configuration only simulates CH 4 and assimilates CH 4 observations. δ 13 C(CH 4 ) initial conditions in S1 are not offset and are directly taken from the spin-up. S2 assimilates the real δ 13 C(CH 4 ) observations instead of the fitting curve data. In S3, 5 the δ 13 C(CH 4 ) model uncertainties are divided by a factor 2. T1 uses 10 sub-categories instead of 5 aggregated categories, increasing the degrees of freedom. 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 10 for each category. In T3, the δ 13 C(CH 4 ) source signatures 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 4 emissions. Finally, T4 applies both changes from T2 and T3. Table 2  3 Results

Minimization of the cost function
The minimization process is performed using the M1QN3 algorithm (Gilbert and Lemaréchal, 1989 burden is increased by a factor 2 in comparison to an inversion without the isotopic constraint due to the doubling of simulated tracers ( 12 CH 4 and 13 CH 4 ). 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  bution to Jo. c) δ 13 C(CH4) contribution to Jo. d) RMSE associated to observed-simulated CH4. e) RMSE associated to observed-simulated δ 13 C(CH4). For clarity reasons, S1 and S3 initial values are not displayed because they are too large compared to REF.
that increasing the degrees of freedom also increases the computational burden. For S1, it highlights the benefits of offsetting δ 13 C(CH 4 ) initial conditions. As we assume no correlation of errors in R, J o (see Eq. 3) can be divided into CH 4 and δ 13 C(CH 4 ) contributions. Figure 3 shows that all configurations lead to a fast reduction of the δ 13 C(CH 4 ) contribution. During the first ten iterations, it decreased from 50-90% (depending on the configuration) to 10-20 %. Conversely, the CH 4 contribution increased from 10-50 % to 80-5 90%. By adjusting the source isotopic signatures (all configurations besides T3-T4), the system was able to efficiently and rapidly reduce the discrepancies between simulated and observed δ 13 C(CH 4 ). As a result, the δ 13 C(CH 4 ) RMSE decreased very rapidly during the first ten iterations while the CH 4 RMSE due to CH 4 discrepancies decreased at a roughly constant rate. Consequently, the system is preferentially adjusting δ 13 C(CH 4 ) over CH 4 values to reduce the cost function. The decrease rate associated with δ 13 C(CH 4 ) RMSE can be increased by reducing the model uncertainties prescribed to the δ 13 C(CH 4 ) 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 less iterations than REF to reach similar δ 13 C(CH 4 ) RMSE reduction but 7 additional iterations to reach similar CH 4 RMSE reduction. T3 and T4 configurations constrains the isotopic signatures, thus the reduction of the δ 13 C(CH 4 ) contribution necessitates 25 more iterations than REF to reach similar RMSE reduction. 5 To summarize, the decrease rate associated with δ 13 C(CH 4 ) RMSE is highly dependent on the prescribed uncertainties in δ 13 C(CH 4 ) observations and the ability of the system to adjust source signatures.
3.2 CH 4 and δ 13 C(CH 4 ) fitting As expected, the assimilation process greatly improves the agreement between simulated and observed values for both CH 4 and δ 13 C(CH 4 ). Figure 4 shows the globally-averaged time-series of CH 4 and δ 13 C(CH 4 ).

10
CH 4 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 lead to simulated CH 4 mixing ratios in good agreement with observations and the improvement is therefore relatively small. In addition, all configuration results regarding CH 4 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 4 observations.

15
Prior δ 13 C(CH 4 ) 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 is likely due to an underestimation (too negative values) of some source isotopic signatures or a poor prior estimation of the source partitioning, i.e. an underestimation of 13 C-enriched sources (fossil fuels or biomass burning) or an overestimation of 13 C-depleted sources (biogenic). The data assimilation process reconciles simulated and observed δ 13 C(CH 4 ) (RMSE is 0.086 ± 0.008 ‰) for all configurations, albeit small differences depending on the setup 20 emerge.
The S-group provides a better match than the T-group (0.081 ± 0.003 ‰ versus 0.091 ± 0.007 ‰). Furthermore, the fit is very similar within the S-group. In contrast, the spread in the T-group is larger with δ 13 C(CH 4 ) 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 regional-specific estimates of prior source signatures instead of global values may   Moreover, T3 provides the poorest δ 13 C(CH 4 ) fitting at AMY (0.24 ‰). Therefore, setting 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 where these stations are located.

Global and regional emission increments
We are primarily interested in the additional information provided by the assimilation of δ 13 C(CH 4 ) data. Rather than discussing the regional and global CH 4 emissions and comparing these results to previous estimates, we investigate the differences   Global increment differences in MC (-6.4 Tg yr −1 ) and FF emissions (8.6 Tg yr −1 ) are mainly due to regional increment 20 difference in China and Temperate Asia. MC regional increment difference is equal to -2.1 Tg yr −1 in Temperate Asia and -2.4 Tg yr −1 in China. Similarly, FF increment difference is equal to 1.5 Tg yr −1 in Temperate Asia and 5.0 Tg yr −1 in China.
WT global increment difference (-5.7 Tg yr −1 ) is mainly due to differences in Canada (-2.0 Tg yr −1 ) and South America (-2.3 Tg yr −1 ) 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 yr −1 principally owing to 25 increment differences in South-East Asia (1.7 Tg yr −1 ), Canada (0.4 Tg yr −1 ) and Africa (0.4 Tg yr −1 ). The Natural (NAT) category exhibit very little changes (less than 1 Tg yr −1 ), even in relative values (see Fig. S3 in the supplement).
S-group configurations infer results remaining 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.   55 -56] estimates. In other words, the quality of isotopic signature (values and distributions) appears to be critical for the robustness of the system's source estimates.

Global and regional source signature increments
Source isotopic signatures are also optimized by the system. Figure 7 provides the flux-weighted source isotopic signatures for different regions. It shows the difference between REF posterior and prior estimates.

5
All source signature are shifted upwards by the inversions, in order to correct the too strong negative trend in δ 13 C(CH 4 ). At the global scale, flux-weighted source signatures of WT, FF, MC and BB are increased by 1.7, 0.5, 0.9 and 0.5 ‰, respectively.
The global source signature is increased from -53.9 ‰ (prior) to -52.6 ± 0.2 ‰ (posterior) depending on the configuration (see Table S6). The posterior global signature is strongly dependent on the total fractionation effect. This effect tends to deplete air in 13 CH 4 , shifting the δ 13 C(CH 4 ) to more positive values as the CH 4 molecules emitted by the sources are removed from the The WT source signature exhibits the larger upward shift, from a global value of -60.8 ‰ to -59.1 ‰. This large difference for an average signature is due to upward shifts in Boreal regions (North America, Russia) but also in South America and 15 Temperate Asia. The MC source signature is increased by 0.9 ‰ mainly due to changes in Asia. The FF source signature is increased by 0.5 ‰ globally due to a large increment in China (+ 1.2 ‰). Finally, the BB source signature is reevaluated in South-East Asia (+ 1.4 ‰) and Canada (+ 0.8 ‰).
These changes are consistent within the S-group (see blue errorbars in Fig. 7), although small variations are visible (e.g. ± 0.3 ‰ for WT in Canada). The source signature is therefore modified nearly to the same extent in all regions, no matter which  (Table S6). T1 (yellow dot in Fig. 7), with more optimized categories than the others, shows small differences at global scale (less than 0.3 ‰ for all categories), although differences of more than 1 ‰ are visible in China. Therefore, increasing the number of degrees of freedom lead to similar flux estimates but can affect the signatures at regional scale.
T2 estimates are shifted upward to reach a less negative global source isotopic signature without getting closer to the regional 5 distribution of the S-group. This is likely caused by the scarcity of δ 13 C(CH 4 ) stations and correcting this behavior seems challenging without additional observations. The problem might be circumvented by using the region scale rather than the pixel scale to optimize isotopic signature values. Future inversions will test this assumption.

Posterior uncertainties
Formally, posterior uncertainties are given by the Hessian of the cost function. This matrix can hardly be computed at an  (Table 3) exhibit an uncertainty of 10 %, 7 %, 19 % and 38 %, respectively. BB is the most uncertain 15 estimate relatively to its intensity, although FF show the largest absolute uncertainty (23 Tg yr −1 ).

Conclusions
We present here a new variational inversion system designed to assimilate observations of both a specific trace gas and its isotopic data. This system allows to optimize both the tracer emissions and the associated isotopic signatures for multiple source categories. To test this system we have assimilated CH 4 and δ 13 C(CH 4 ) data retrieved at different measurement sites 20 over the globe.
Different configurations have been tested in order to assess the sensibility of the system to the setup. We have shown that offsetting the δ 13 C(CH 4 ) initial conditions before the inversion (S1), using δ 13 C(CH 4 ) curve fitting data instead of the original observations (S2) and reducing the prescribed uncertainties in the δ 13 C(CH 4 ) observations (S3) have very little effect on the inferred fluxes (less than 2 Tg yr −1 for each category at global scale). However, offsetting the δ 13 C(CH 4 ) initial conditions Tg yr −1 ), the posterior isotopic signatures can be modified in some regions (more than 1 ‰ in China). Also, starting from 30 mean global values for the source signatures (T2) makes the system unable to retrieve the regional-specific isotopic signatures from REF. Increasing the number of δ 13 C(CH 4 ) observations could help to cope with this issue. Finally, configurations which constrain the source signatures (T3-T4) show differences in global flux estimates of more than 10 Tg yr −1 , compared to REF.
This emphasizes the need for good prior source signature estimates.
The major caveat 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 8 CPUs can be run in parallel, resulting in an elapsed time of 5-6 weeks to run one of the inversions of 5 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. In addition, variational inversions as implemented in the CIF are not enabled to provide a quantification (even approximated) of the posterior uncertainties. Dedicated efforts need to be done to address this issue in the future, at an achievable numerical cost.
This system is implemented within the CIF framework and can therefore be used for inversions with the various CTMs 10 embedded in the CIF, provided the adjoint codes of the models exist. 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 here the sensitivity of the system to these parameters as only 15 technical aspects of the system were tested. This will be part of a future analysis.
δ 13 C(CH 4 ) is not the only isotopic data that can be assimilated in such a system. Many δD(CH 4 ) observations have also been retrieved during the period 2004-2010 at many different locations. These isotopic values can provide additional information that can further help to discriminate the co-emitted CH 4 fluxes (Rigby et al., 2012). Moreover, ethane (C 2 H 6 ) is co-emitted with CH 4 by fossil fuel extraction and distribution (Kort et al., 2016;Smith et al., 2015) and observations are available at a 20 multitude of sites since the early 1980s. Therefore, assimilating this data can provide additional constraint. The system will therefore be improved in the future in order to assimilate δ 13 C(CH 4 ), δD(CH 4 ) and C 2 H 6 observations together.
Data availability. The CIF codes and documentation pages are available here: community-inversion.eu.
Author contributions. JT implemented the variational inversion system within the CIF with the precious help of AB. JT designed, run and analyzed the tested configurations. AB, MS, IP and PB provided scientific and technical expertise. They also contributed to the analysis of 25 this work. BV and SEM provided scientific expertise regarding δ 13 C(CH4) observations. JT prepared the manuscript with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.