Natural and anthropogenic methane ﬂuxes in Eurasia: a mesoscale quantiﬁcation by generalized atmospheric inversion

. Eight surface observation sites providing quasi-continuous measurements of atmospheric methane mixing ratios have been operated since the mid-2000’s in Siberia. For the ﬁrst time in a single work, we assimilate 1 year of these in situ observations in an atmospheric inversion. Our objective is to quantify methane surface ﬂuxes from anthropogenic and wetland sources at the mesoscale in the Siberian lowlands for the year 2010. To do so, we ﬁrst inquire about the way the inversion uses the observations and the way the ﬂuxes are constrained by the observation sites. As atmospheric inversions at the mesoscale suffer from mis-quantiﬁed sources of uncertainties, we follow recent innovations in inversion techniques and use a new inversion approach which quantiﬁes the uncertainties more objectively than the previous inversion systems. We ﬁnd that, due to errors in the representation of the atmospheric transport and redundant pieces of information, only one observation every few days is found valuable by the inversion. The remaining high-resolution quasi-continuous signal is representative of very local emission patterns difﬁcult to analyse with a mesoscale system. An analysis of the use of information by the inversion also re-veals that the observation sites constrain methane emissions within a radius of 500 km. More observation sites than the ones currently in operation are then necessary to constrain the whole Siberian lowlands. Still, the ﬂuxes within the constrained areas are quantiﬁed with objectiﬁed uncertainties. Finally, the tolerance intervals for posterior methane ﬂuxes are of roughly 20 % (resp. 50 %) of the ﬂuxes for anthropogenic (resp. wetland) sources. About 50–70 % of Siberian lowlands emissions are constrained by the inversion on average on an annual basis. Extrapolating the ﬁgures on the constrained areas to the whole Siberian lowlands, we ﬁnd a regional methane budget of 5–28 TgCH 4 for the year 2010, i.e. 1–5 % of the global methane emissions. As very few in

Abstract. Eight surface observation sites providing quasicontinuous measurements of atmospheric methane mixing ratios have been operated since the mid-2000's in Siberia. For the first time in a single work, we assimilate 1 year of these in situ observations in an atmospheric inversion. Our objective is to quantify methane surface fluxes from anthropogenic and wetland sources at the mesoscale in the Siberian lowlands for the year 2010. To do so, we first inquire about the way the inversion uses the observations and the way the fluxes are constrained by the observation sites. As atmospheric inversions at the mesoscale suffer from mis-quantified sources of uncertainties, we follow recent innovations in inversion techniques and use a new inversion approach which quantifies the uncertainties more objectively than the previous inversion systems. We find that, due to errors in the representation of the atmospheric transport and redundant pieces of information, only one observation every few days is found valuable by the inversion. The remaining high-resolution quasicontinuous signal is representative of very local emission patterns difficult to analyse with a mesoscale system. An analysis of the use of information by the inversion also reveals that the observation sites constrain methane emissions within a radius of 500 km. More observation sites than the ones currently in operation are then necessary to constrain the whole Siberian lowlands. Still, the fluxes within the constrained areas are quantified with objectified uncertainties. Finally, the tolerance intervals for posterior methane fluxes are of roughly 20 % (resp. 50 %) of the fluxes for anthropogenic (resp. wetland) sources. About 50-70 % of Siberian lowlands emissions are constrained by the inversion on average on an annual basis. Extrapolating the figures on the constrained areas to the whole Siberian lowlands, we find a regional methane budget of 5-28 TgCH 4 for the year 2010, i.e. 1-5 % of the global methane emissions. As very few in situ observations are available in the region of interest, observations of methane total columns from the Greenhouse Gas Observing SATellite (GOSAT) are tentatively used for the evaluation of the inversion results, but they exhibit only a marginal signal from the fluxes within the region of interest.

Introduction
Methane (CH 4 ) in the atmosphere contributes to climate forcing as a greenhouse gas and is involved in the atmospheric oxidizing capacity (Forster et al., 2007). Characterizing the variability of the atmospheric CH 4 composition requires accurate understanding of the methane biogeochemical cycle, in particular of the surface-atmosphere fluxes, of their spatial distribution and of their temporal variability. The quantification of these contributions to the methane cycle still experiences high uncertainties (Kirschke et al., 2013). The global surface to atmosphere CH 4 fluxes range between 500 and 600 TgCH 4 yr −1 (1 Tg = 10 12 g). Two of the main contributors to the global CH 4 budget are natural emissions from inundated areas and anthropogenic sources from coal, oil and gas extraction and distribution. Inundated areas are responsible for 145-260 TgCH 4 yr −1 , i.e. 25-50 % of total emissions with a very high heterogeneous spatial distribution and yearto-year variability (e.g. Bousquet et al., 2006;Dlugokencky et al., 2009;Bergamaschi et al., 2009). The anthropogenic sources from fossil fuel burning and leaks account for 100-150 TgCH 4 yr −1 , i.e. 20-30 % of total emissions, according to the EDGAR inventory (depending on the year and the version of the inventory; http://edgar.jrc.ec.europa.eu).
The West Siberian Plain concentrates significant sources of CH 4 of both wetland and anthropogenic types (Lechtenböhmer et al., 2005;Spahni et al., 2011). On the one side, with 50-70 % of its area covered by peatlands (Peregon et al., 2009), about 13 % of global wetlands are located in the West Siberian Plain. On the other side, Russia produces 20 % of the natural gas in the world, mostly extracted in Siberia; and 0.1 to 10 % (i.e. 0.5-40 TgCH 4 yr −1 ) of this gas is estimated to leak into the atmosphere (e.g. Hayhoe et al., 2002). Large amounts of methane are also released during the oil welling, of which Russia is also a major producer (∼10-15 % of the global production), with 1-2 % of the oil production leaked into the atmosphere as methane (e.g. U.S. Environment Protection Agency, 2011). Documenting the emissions of methane in the West Siberian Plain is thus critical to reduce the uncertainties on the global methane budget.
However, accurately quantifying the wetland and anthropogenic emissions in the West Siberian Plain is challenging. On the one hand, wetland emissions at high latitudes like in Siberia exhibit a clear year-to-year variability  and a distinct seasonality (e.g. Pickett-Heaps et al., 2011) due to high sensitivity to the soil temperature and humidity, to the water table depth and to the total inundated surface which can vary up to 25 % from year to year (Ringeval et al., 2010). This high sensitivity to the local meteorological parameters could cause still unobserved drastic increases of CH 4 emissions from boreal wetlands with climate change (Bohn et al., 2007). On the other hand, anthropogenic sources are mainly related to uncontrolled leakage which is difficult to estimate. Quantifying these leaks raises many issues: pipelines with tiny leaks from chemical perme-ability span over thousands of kilometres, single leaks range on a spectrum of several orders of magnitude (from the lower with pipeline permeability leaks, to the higher during the welling), and leaks can appear (and disappear when detected and repaired) very quickly.
Despite the importance of quantifying and understanding the contribution of the West Siberian Plain to the global CH 4 budget, few studies have been dedicated to this region. Glagolev et al. (2011) carried out extensive field measurements of local wetland CH 4 emissions in order to characterize the emission patterns of each different environment. They upscaled their results to the whole region using wetland distribution maps. However, considering the discrepancies on the wetland distribution (Peregon et al., 2009;Frey and Smith, 2007) and their extension variability (Ringeval et al., 2010), numerous hard-to-quantify uncertainty sources are expected from this approach. For the quantification of the anthropogenic sources, Dedikov et al. (1999) for instance measured mixing ratios of methane close to gas lines and gas facilities to deduce emission factors. They upscaled their figures to the Russian territory and got an emission factor for CH 4 lower than 1 % of the total production of natural gas. Reshetnikov et al. (2000) reviewed the existing literature about Siberia and found emission factors from 0.4 to 14 % of the total production. Another approach is the analysis at different places and dates of the local variations of the atmospheric composition (mixing ratios and isotopic fraction) in CH 4 (and related species). The variations of the atmospheric composition provides information on the relative contribution of the different local processes in the Siberian budget. Such an analysis has been carried out using observations from mobile platforms, such as aircraft profiles (Yamada et al., 2005;Umezawa et al., 2012) or train and ship measurements (Tarasova et al., 2006), and fixed stations (Sasakawa et al., 2012). Nevertheless, these studies are not systematic and comprehensive: they give local instantaneous information and exclusively knowledge about the relative contributions of particular processes to the total regional fluxes. High spatial and temporal resolutions, absolute and consistent quantification and separation of the different main sources is essential to better characterizing and estimating the contribution of Siberian emissions to the global methane budget.
A first step towards this goal was done by Winderlich (2012). In this work, a systematic analysis of the variability of the atmospheric composition at the ZOTTO (Zotino Tall Tower Observatory) tall tower (described in Winderlich et al., 2010) and at a small set of auxiliary sites was carried out. Their approach relied on atmospheric inversion techniques based on the data assimilation Bayesian theory (Enting et al., 2012;Tarantola, 1987) and in principle allows an objective use of the information (here the variability of the atmospheric composition). This provided insights on, e.g. the likely under-estimation of the local anthropogenic sources from gas leak in some inventories (EDGAR v4.1 in Winder-A. Berchet et al.: CH 4 flux in Eurasia 5395 lich, 2012). However, regional atmospheric inversions with only a few sites give only local constraints on the surface fluxes. More critically, a regional inversion with a small number of observations experiences difficulties in identifying and separating the different contributions to CH 4 emissions. This mis-separation can be related to mis-quantified sources of errors in the atmospheric inversion in addition to the lack of information. This issue is dominant in the West Siberian Plain because of the co-located wetland and anthropogenic emissions. The objective of quantifying and separating the regional Siberian sources requires an inversion based on numerous observation sites, with a comprehensive approach of quantification of the uncertainties.
Recent literature highlights the need for precisely and objectively quantifying all errors in the inversion (transport, representation, flux distribution, etc.) in order to evaluate their impact on the inversion results (Lauvaux et al., 2009;Winiarek et al., 2012;Berchet et al., 2013b;Ganesan et al., 2014). In Berchet et al. (2015), we proposed a general method in order to objectively quantify most of the critical sources of errors in the inversion. This improved algorithm is based on a Monte Carlo approach superimposed to maximum likelihood estimators (Chapnik et al., 2004;Michalak and Kitanidis, 2005).
For the first time in Eurasia, we use this improved algorithm on a network of eight surface sites (Sasakawa et al., 2010;Winderlich et al., 2010) covering a large part of the Siberian lowlands and of five remote sites that constrain the air masses coming into the domain and getting out of it. These sites, which have been operated since the mid-2000s, are implemented into the inversion system with objectified uncertainty quantification from Berchet et al. (2015). Here, our goal is to deduce an accurate quantification of the fluxes at the mesoscale with a temporal resolution of a few days from the variability of the atmospheric CH 4 composition at the 8+5 observation sites.
We explain the theoretical background used in our study in Sect. 2. In Sect. 3, the data sets and models used in the inversion are introduced. We then present the results on the fluxes and the limitations of the inversion in Sect. 4 and Sect. 5. Our inversion is then evaluated in Sect. 6 by using in situ measurements and independent satellite observations, as the very few available surface observations do not allow to keep enough of them out of the inversion for evaluation.

Motivations towards marginalizing
As the atmosphere mixes irreversibly air masses from different CH 4 sources, using the atmospheric signal as in an atmospheric inversion cannot lead to a deterministic characterization of the surface fluxes. In the classical Bayesian framework (Tarantola, 1987), the objective of the inversion is to inquire about the probability density function (pdf) of the surface fluxes, or more generally of the state of the system, with some knowledge about the atmospheric composition and on the prior state distribution. The sought pdf can be written: In this formula, the vector y o gathers all the available observations; H is the observation operator converting the information in the state vector to the observation space; the vector x depicts the state of the system (mainly CH 4 surface fluxes in our case) and x b is the background vector including the prior knowledge on the state x of the system. H typically represents the discretization of the problem to be computed numerically and the atmospheric transport from the emission areas to the observation sites. In the following, we consider H linear and associate it to its Jacobian matrix H. In our case, we simulate the atmospheric transport only on a domain of limited area as detailed in Sect. 3.3 and 3.4. Additional information about the atmospheric composition at the boundaries of the domain of interest is then necessary to compare observed with simulated atmospheric composition. Therefore, the state vector x encompasses the surface fluxes, but also the lateral boundary conditions related to the observed baselines at each observation site.
Thus, the inversion computes the pieces of information contained within the observations, the prior state and the representation of the transport with their associated uncertainties (e.g. measurement errors, uncertainties in the flux inventories, etc.). With the usual Gaussian assumption, all the uncertainties are considered as normal pdfs and can be described with modes and uncertainty covariance matrices. The inversion then deduces an optimal posterior state vector x a and posterior uncertainties P a . Within the Gaussian assumption, the posterior state vector and uncertainty matrix can thus be explicitly defined: In Eq. (1), the matrix K = BH T (R+HBH T ) −1 is the Kalman gain matrix. The matrices R and B are the covariance matrices describing the observation and background uncertainties. Observation uncertainties encompass measurement, discretization and transport errors. Background uncertainties include the uncertainties in the spatial distribution of the fluxes, in their temporal variability and in their absolute value. As long as these uncertainty matrices are known, the inversion only faces technical issues (e.g. matrix inverses and products in large dimension problems) for numerical implementation. However, only the uncertainties in the measurements are objectively quantified during the calibration process. The errors in the transport or in the prior fluxes are not perfectly known and, in most cases, they are built relying on some expert knowledge about the system. But this subjective knowledge can lead to ill-specified matrices, which have a dramatic impact on the inversion results (e.g. Cressot  et al., 2014). Recent studies inquired into objectified ways of specifying these matrices (e.g. Michalak and Kitanidis, 2005;Winiarek et al., 2012;Berchet et al., 2013b). The approach in these papers was to find optimal uncertainty matrices R and B along an objective statistical criterion: the maximum likelihood. The implementation of the method gave encouraging results, but the impact of the uncertainties within the maximum likelihood computation were complicated to evaluate. Berchet et al. (2015) used a general marginalization approach in order to quantify objectively all the uncertainties impacting on the posterior fluxes. In the following, we summarize this approach.

Method outline
The marginalization consists in computing the complete pdf p(x|y o − Hx b , x b ) as a weighted sum of the Gaussian pdfs p(x|y o − Hx b , x b , R, B) calculated for each possible uncertainty matrices R and B. This can be written as follows: To compute the marginalized integral in Eq.
(2), a large number (60 000 in our case) of posterior vectors x a and posterior uncertainty matrices P a is computed through individual inversions as in Eq. (1) with different uncertainty matrices R and B. This Monte Carlo sampling is carried out based on an estimate of the pdfs of the uncertainty matrices, p(R) and p(B). This estimate is deduced from the objectified maximum likelihood approach which produces a first guess for uncertainty matrices (R max , B max ).
Posterior uncertainties and correlations are in the end computed from the ensemble of 60 000 computed individual inversion results as follows: with N the number of Monte Carlo draws. As detailed in Sect. 2.1, the main motivation for marginalizing the classical inversion framework is to use the available information (in situ observations, flux inventory, transport model) in a way which is as objective as possible. The marginalized inversion gives an explicit and objectified access to pieces of information (described in details in Sect. 2.3) required: 1. to evaluate the efficiency of the observation network for constraining regional emissions and to give guidelines for improving monitoring deployment; 2. to inquire about atmospheric inversion skills in terms of resolved temporal and spatial resolution and of emission process separation, and to deduce observation and modelling requirements for future better inversions; and 3. to assess the robustness of emission inventories and process-based surface-atmosphere exchange models.
These three key points are further discussed in Sect. 4. One of the main drawbacks in our method is its numerical cost. Computing 60 000 individual inversions is cumbersome and requires extensive amounts of memory. Smartly chosen filtering criteria, observation sampling and flux aggregation patterns must be carried out before the marginalized inversion in order to reduce its complexity without degrading the inversion optimality. As recent efforts have been made to reduce the subjectivity in the inversion, we rely on objective criteria (though not computed in an exact manner due to computational limitations) to complete the observation sampling, flux aggregation and filtering. These criteria are shortly developed in Sect. 2.4 and detailed in Berchet et al. (2015).

Output analysis
As described above, the marginalized regional inversion as we compute it answers three main questions. Below are given details on how these questions are treated.

Network efficiency
The evaluation of the network coverage is carried out through the explicit computation at the maximum likelihood of the influence K max H and sensitivity matrices HK max following Cardinali et al. (2004), with K max the Kalman gain matrix defined as in Eq. (1) with R max and B max . With these two matrices, we can compute the unconstrained emissions ǫ and the weight ω of each observation site i in the inversion as the sum of the weight ω i,j of each observation j carried out at the site i: with φ i the emissions at the pixel i. In Eq. (4a), the higher the score ω i the more the inversion uses the site i. Observation sites downwind emissions will have a strong impact on the inversion, but observations constraining air masses coming into the domain of interest are also key stones for regional inversions. Eq. (4b) depicts emitted CH 4 that is not seen or constrained by the inversion. Regions with the highest ǫ i are unseen areas with strong emissions (according to the inventories), where additional observation sites would be required. The colour bar paints the altitude in metres above sea level (from ETOPO1 database; Amante and Eakins, 2009). Red dots (resp. orange triangle) depict hotspots of CH 4 emissions (based on EDGAR v4.2 inventory; see Sect. 3.3) related to oil welling and refineries (resp. gas extraction and leaks during distribution in population centres). Purple squares highlight the observation site localizations. Blueish shaded areas represent inundated regions, wetlands and peatlands (from the Global Lakes and Wetlands Database; Lehner and Döll, 2004). The Siberian budget in Sect. 5.2 is calculated on the hatched area.

Solved spatial and temporal scales
The spatial and temporal resolutions that the inversion can solve are described by posterior error covariances. Misseparated regions are usually detected through so-called flux dipoles (e.g. Rödenbeck et al., 2003). In our case, as we explicitly and objectively compute the posterior matrix P a from the Monte Carlo ensemble x a , posterior correlations are used for post-processing groups of ill-separated regions (regions are defined through a dynamically chosen aggregation pattern prior to the inversion; see Sect. 2.4.2). Both strongly positive and negative correlations r i, j point to illseparated regions i and j . For instance, some neighbouring emissions or successive ones could be mixed in the atmospheric signal and then be not separated from the inversion point of view. Following the observing system simulation experiments (OSSEs) carried out in Berchet et al. (2015), we group two posterior flux regions i and j when This post-processing makes it possible to group posterior fluxes in order to avoid dipoles and reduce uncertainties on bigger regions. Thus, it is also possible to filter out regional fluxes that are not separated from the boundary conditions. Every group of correlated regions including some contribution from the boundary conditions are excluded from the analysis of inversion results, because errors on the lateral boundary conditions can mislead the inversion about the regional fluxes. Through this post-grouping, we can assess the typical size of aggregated regions that the inversion can constrain with our set-up (observations and transport model resolution). The typical timescale on which the inversion can apply increments to the emission can also be assessed with this postprocessing. As we want to separate contributions from anthropogenic and wetland emissions, the post-grouping may group or not co-located emissions from different processes and then give insights into the separation ability of our regional inversion.

Posterior flux analysis
Ideally, an atmospheric inversion provides insights about emissions. From the correlation grouping applied to the Monte Carlo ensemble of posterior fluxes, we can compute tolerance intervals of posterior fluxes so that 68.27 % of the ensemble is within the interval. The number 68.27 % makes the tolerance intervals equivalent to the ±σ interval in the Gaussian framework, as the Monte Carlo posterior ensemble is not necessarily a Gaussian distribution. The inversion thus indicates that the fluxes we are inquiring into are very likely in the defined posterior tolerance interval. Deviating posterior tolerance intervals compared to prior fluxes point at required updates in the used prior database. Below we present our visualization approach to control posterior fluxes. Figure 5 synthesizes the inverted methane fluxes for Siberian lowlands (hatched domain in Fig. 1). As detailed in Sect. 2.3.2, anthropogenic emissions (inverted at a monthly scale) can be grouped with wetland emissions (considered at the weekly scale). So, the lowest common multiple on which the fluxes can be analysed is the monthly scale. Then, for the Siberian lowlands, for each month, we define the proportion The proportions of anthropogenic and wetland emissions that are constrained, constrained but mixed with another type of emission or unseen by the observation network are represented in the pie charts in Fig. 5. Finally, within the proportion of constrained regions, we analyse the inversion correction on anthropogenic, wetland and mixed emissions (bar diagrams in Fig. 5). For each type of emission (anthropogenic, wetland and mixed), we present in the bar diagrams of Fig. 5 the tolerance intervals of posterior fluxes and the prior uncertainties as calculated by the maximum likelihood algorithm (that is to say B max ; see Sect. 2.2). For each month, the tolerance intervals on the total prior and posterior methane budget in the Siberian lowlands is also given in TgCH 4 (upper left corner of each graph).

Size reduction and filters
As suggested in Sect. 2.2, the marginalized inversion requires some filtering, observation sampling and flux aggregation, so it can be numerically computed. Below, we explain how we carry out these pre-processing in a way that do not dampen the advantages of the marginalization.

Observation sampling
Here we try to reduce the dimension of the observation space. At the regional scale, considering the spatial resolution of our transport model, only the synoptic variability of the observed signal is relevant. We then decide to keep only one piece of information per site and per day as it is commonly done at the global and continental scales. In addition, simulated vertical mixing close to the surface where observations are carried out is known to be flawed when the planetary boundary layer (PBL) is shallow (typically at night and in winter in Siberia). We then sample the observations during the afternoon when the PBL is higher than 500 m as suggested by prior studies (e.g. Berchet et al., 2013b) and we pick the observed and simulated mixing ratios at the time when the observations are minimum.
As surface emissions dominate surface sinks for CH 4 , keeping the minimum observed mixing ratio by afternoon is equivalent to detecting the time when the PBL is at its maximum, hence when the atmospheric model is the more accurate. This criterion filters out outliers generated by local influences which cannot be reproduced by an atmospheric transport model with a resolution larger than 25 km, and which only add noise to the system.
For our case study, out of 127 000 hourly measurements available in 2010, 30 000 pass through the PBL height and night filters (see black dots in Fig. 2). Out of these 30 000 data points, about 2000 daily aggregates are selected. Details by observation site are given in Table 2 (last two columns).

Flux aggregation and constraints
The following procedures are meant to decrease the state space dimension. To define aggregation patterns, we divide our domain into 35 physical regions for each emission process (according to vegetation types, demography, industrial activity, etc.) as a basic pattern. This basic pattern is chosen so that the mesh gets finer closer to the observation network. To further reduce the number of aggregated regions, under-constrained regions are grouped together. This is carried out by analysing the observation network footprints estimated with a Lagrangian model (Sect. 3.2) which offers an efficient way to compute them.
In addition to the footprint aggregation, the influence matrix K max H makes it possible to quantify observational constraints on the fluxes. Below a given threshold of constraint for a flux (related to the flux contribution to the atmospheric signal), the inversion cannot deduce any valid information on the flux. For this reason, we also filter out a region i with very low constraints from the marginalization if

Plume filtering
It is known that atmospheric transport models suffer from temporal and spatial mismatches when simulating air masses with strong mixing ratio gradients (e.g. for plumes welldelimited from the background air masses). When a plume is transported to the wrong place and time, as we sample air masses at a given point location, the differences between simulated and observed mixing ratios can reach very high values. Such strong model-observation differences lead to unrealistic corrections on fluxes.
To dampen such undesirable effects, the resolution of the transport representation is chosen finer close to the observation network. In addition, we introduce a procedure to filter out plume-like air masses from the inversion input. As explained in Sect. 2.2, a maximum likelihood estimation is computed prior to the marginalization. We take advantage of the maximum likelihood computation to detect air masses critically ill-reproduced by the transport model. As the maximum likelihood estimation computes optimal uncertainty matrices R max and B max , we filter out observations with too a high computed uncertainty. That is to say, for each observation i, the data point is excluded if This criterion does not necessarily exclude only plumes generating a big signal, but only the ones that are very poorly reproduced by the transport model. One should notice that this criterion is computed in association with the low constraint criterion of Eq. (5). That is to say, a region which always illuminates the observation network through plume-like air masses will have all its constraining observations filtered out. As a consequence, it will Table 1. Site characteristics. The altitudes of the sites are given as m above sea level (a.s.l.) and the inlet height is in m above ground level (a.g.l.). The frequency column depict the type of measurements in the site: C = quasi-continuous, F = flask sampling.  All the mixing ratios are reported to the NOAA 04 scale before being implemented into the inversion system. JR-STATION and ZOT sites are located in the vicinity of anthropogenic and wetland sources. These local sources strongly influence the Siberian network, which was designed to monitor regional emissions. As a consequence, numerous observations from the Siberian network are ill-reproduced by our transport model and then are filtered out from the inversion as "plume" observations (following criteria in Sect. 2.4; see also Fig. 2 and Table 2 in Sect. 6.1). In particular, measured mixing ratios at BRZ site are not well simulated by our model, possibly due to missing local emissions in the prior. Most BRZ observations are thus filtered out from the inversion. Due to logistical issues or instrument dysfunctions, observation sites do not provide measurements all year round. Fig. 2, described in Sect. 4.1, details the temporal availability of the observations. The sampling bias is known to impact the inversion results (Villani et al., 2010). The issue is discussed in Sect. 4.1, but the general method developed in Sect. 2 consistently takes into account such a bias regarding increased posterior tolerance intervals and decreased constraints on the emissions.

Biogeosciences
The observation vector y o is defined after Sect. 2.4.1 sampling method. The final size of y o implemented in the inversion is 2000. On average, 0.4 observations per station per day are validated for the inversion.

Estimates of the network footprints
As the observations that will be implemented in the system are known, the observation network footprints, necessary to choose the aggregation patterns in order to define the prior state vector x b and the observation operator H (as detailed in Sect. 2.4.2), can be computed. As we do not carry out a quantitative analysis of the footprints, we only need a rough estimation of the footprint patterns. Thus, we compute simulations with the Lagrangian dispersion model FLEXPART version 8.2.3 (Stohl et al., 2005) to get such an estimation. To build the footprints, we compute numerous back trajectories of virtual particles from the observation sites at the times when measurements are available and valid for the inversion.
The model is forced by ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim data at an horizontal resolution of 1 • × 1 • , with 60 vertical levels and 3 h temporal resolution (Uppala et al., 2005). Virtual particles are released in a 3-D box (10 km per side and 1000 m high) centred around each observation site with 10-day lifetime backwards in time. The footprints are computed on a 0.5 • ×0.5 • horizontal grid, following the method of Lin et al. (2003), taking into account the boundary layer height at each particle location. This method considers that only the particles within the boundary layer are influenced by surface emissions and that the boundary layer is well-enough mixed to be considered as uniform.

Prior fluxes and state vector: x
The inversion system optimizes prior fluxes grouped in the regions aggregated by the pre-processing procedure (see Sect. 2.2). The prior spatial distribution and temporal variability of the fluxes are deduced from the following: 1. EDGAR database v4.2 FT2010 (http://edgar.jrc.ec. europa.eu) for year 2010 for anthropogenic emissions, 2. LPX-Bern v1.2 process model (Stocker et al., 2014) at a monthly scale for wetland emissions 3. GFED v4 database at a daily scale for wildfires.
In Fig. 1, the distributions of the anthropogenic hotspots of emissions and of the wetland regions are represented, superimposed over the regional topography. Anthropogenic emissions of methane in the region are mainly hotspots related to the intense oil and gas industry in the Siberian lowlands and to the leaks in the distribution system in population centres in the vicinity of the Trans-Siberian Railway in the southern part of the Siberian plain. Wetland emissions are mainly confined to the lower part of Siberia in the West Siberian Plain, half of which is lower than 100 metres above sea level. Wildfires occur mainly in spring and summer in the Eurasian forest-covered areas; they emit CH 4 as intense hotspots. The EDGAR inventory uses economic activity maps by sectors and convolved with emission factors calculated in laboratories or with statistical studies (Olivier et al., 2005). The Bern-based land surface processes and exchanges (LPX-Bern v1.2) model is an update of the dynamic global vegetation model LPJ (Lund-Potsdam-Jena Dynamic Global Model) (Spahni et al., 2011). It includes a dynamical simulation of inundated wetland areas (Stocker et al., 2014), dynamic nitrogen cycle , and dynamic evolution of peatlands Stocker et al., 2014). The model uses CRU TS (Climate Research Unit time series) 3.21 input data (temperature, precipitation rate, cloud cover, wet days), observed atmospheric CO 2 and prescribed nitrogen deposition (Lamarque et al., 2011) for each year for the simulation of dynamic forest and peatland vegetation growth. The GFED v4 database is built from the 500 m Collection 5.1 MODIS DB burned-area mapping algorithm (Giglio et al., 2009). CH 4 emissions at monthly and daily scales are deduced from the burnt areas using the Carnegie-Ames-Stanford Approach (CASA model;Potter et al., 1993) and emission factors (van der Werf et al., 2010).
We are aiming at a separation of the types of emissions at the mesoscale. We therefore aggregate the emissions along the three different types of sources, with specific spatial patterns and temporal profiles for each type of emissions. Anthropogenic sources are hotspots that emit all year round. Wetlands are responsible for diffuse fluxes on large areas, with a high temporal variability depending on the local weather conditions (typically temperature or water table depth). The emissions of CH 4 from wildfires come from point sources and occur on relatively short periods (Kasischke and Bruhwiler, 2002). Consequently, we do not aggregate the different types of emissions along the same spatial patterns and temporal intervals. Anthropogenic emissions are aggregated by month, while wetlands and wildfires, which have quicker time responses to meteorological changes, are grouped by periods of 10 days. In the following, we discuss the inversion results only in term of anthropogenic and wetland emissions. Indeed, as the wildfire emissions generate plumes relatively well-defined from the ambient air, the marginalized inversion exclude from the system all the emis-  Table 2. Correlations of observed and simulated mixing ratios, prior and posterior to the inversion. Posterior correlation coefficients r are presented for the filtered data points used in the inversion, but also for those filtered out. The number of daily available observations is also reported alongside with the number of data used in the inversion. sion contribution related to fires (according to the procedures described in Sect. 2.2). For the computation of the observation operator H (see Sect. 3.4), we use a regional chemistry-transport model with a domain limited in space and time. Initial and lateral boundary conditions (hereafter ICs and LBCs) are then also to be optimized in the system. Prior lateral mixing ratios are deduced from simulations at the global scale by the general circulation model LMDz with the assimilation of surface observations outside the domain of interest (Bousquet et al., 2011). LBCs are assimilated by periods of 10 days. We arbitrarily aggregate LBCs along four horizontal components (by side of the domain) and two vertical ones (1013-600 and 600-300 hPa). Though we are mainly focused on Siberian lowlands, the domain of model computation has been chosen spanning over all Eurasia. This is expected to attenuate the impact of the rough global resolution in LMDz boundary conditions on the simulated variations of mixing ratios at the observation sites. Indeed, the central region is thousands of kilometres away from the sides of the domain.
To summarize, all the pieces of information in the observations are assimilated to constrain 1700 aggregated regions of flux and boundary conditions: 10 × 12 month regions for anthropogenic emissions, 10 × 36 10-day period for wildfires, 25 × 36 10-day period for wetlands, 9 (4 sides × 2 horizontal levels + roof top) × 36 10-day period for the lateral boundary conditions. After the filtering of Sect. 2.4, the dimension of the state space is reduced from 1700 to 495.

The observation operator: H
We explicitly define the observation operator H by computing the forward atmospheric transport from the regions of aggregated emissions (defined in Sect. 3.3) to the observation sites. As CH 4 is a reactive species, the observation operator should include the oxidation by OH radicals. However, as the residence time of the air masses in the domain of simulation is short (a few days to a few weeks) compared to CH 4 atmospheric life time (8-10 years; Dentener et al., 2003), ignoring OH sinks only generates small differences in the simulated mixing ratios. Additionally, OH sinks are mostly responsible for large-scale gradients while the regional inversion focuses on the synoptic scale. Thus, the regional inversion system attributes OH sinks to global boundary conditions, and not to regional fluxes.
Thus, for each aggregated region, we calculate the socalled response functions using the transport module of the Eulerian mesoscale non-hydrostatic chemistry transport model (CTM) CHIMERE Menut et al., 2013). This model was developed in a framework of air quality simulations (Schmidt et al., 2001;Pison et al., 2007), but is also used for greenhouse gas studies (Broquet et al., 2011;Berchet et al., 2013b). We use a quasi-regular horizontal grid zoomed near the observation sites after Sect. 2 considerations. The domain of interest is of limited area and spans over the mainland of the Eurasian continent (see Fig. 3). As we are interested in mesoscale fluxes, we take a spatial resolution larger than 25 km. The average side length of the grid cells is 25 km close to the western Siberian stations and 150 km away of the centre of the domain. lated meteorological fields from the ECMWF with a horizontal resolution of 0.5 • × 0.5 • every 3 h (Uppala et al., 2005).

Independent observations for evaluation
Any inversion has to be confronted to independent data in order to evaluate its results. Few in situ measurements of CH 4 mixing ratios are available in Siberia. We choose to assimilate all surface observations described in Sect. 3.1 for the optimization of CH 4 fluxes. NIES (National Institute for Environmental Studies, Japan) and LSCE (Climate and Environment Sciences Laboratory, France) carry out aircraft measurements in the region (Paris et al., 2010;Umezawa et al., 2012;Berchet et al., 2013a), but these measurements are still difficult to compare to mesoscale models. In addition, their spatial and temporal coverage is poor for the year 2010, and they are not numerous enough to get significant validation insights.
For year 2010, the only remaining observations with sufficient spatial coverage and temporal availability are the total columns retrieved by the Greenhouse Gas Observing SATellite (GOSAT). In Sect. 6.2, we evaluate the results of the inversion against GOSAT data. The Japanese satellite GOSAT was launched by the Japan Aerospace Exploration Agency (JAXA), NIES and the Japanese Ministry of the Environment (MOE) in January 2009. It has a polar sun-synchronous orbit at 667 km and provides a full coverage of the Earth every 3 days with a swath of 750 km and a ground pixel resolution of 10.5 km at nadir. The TANSO-FTS instrument observes the solar radiation reflected at the surface and the top of the atmosphere in the short wave infrared (SWIR) domain that allows deducing total columns of methane (XCH 4 ) in cloud-free and sunlight conditions. We use version 3.2 of the TANSO-FTS bias-corrected XCH 4 proxy retrievals performed at the University of Leicester (Parker et al., 2011). The XCH 4 retrieval algorithm uses an iterative retrieval scheme based on Bayesian optimal estimation and associated to averaging kernels and a priori profiles. The retrieval accuracy is estimated to be about 0.6 % (i.e. ∼ 10 ppb). The retrieval algorithm needs CO 2 mixing ratios as a proxy for the light path. We use the 4-D CO 2 analysis from the surface air-sample inversion by Chevallier et al. (2010) (MACC v10.2). We obtain ∼ 25000 GOSAT observations in 2010 over the domain of interest.
In order to compare the observations of the total columns to the model, we use the averaging kernels to compute prior and posterior model equivalents. The regional CTM CHIMERE has a top layer at 300 hPa in our set-up. The stratospheric contribution to the total columns is deduced from the global model LMDz used for the initial and lateral boundary conditions (described in Sect. 3.3). The average observed XCH 4 is ∼ 1775 ppb over the domain throughout the year. The prior average in LMDz XCH 4 is ∼ 1820 ppb. This bias is attributable to the excessive injection of tropospheric air into the stratosphere in our version of LMDz. In Sect. 6.2, the bias on XCH 4 of 45 ppb has been corrected to allow the observation-model comparison.

Diagnostics of the marginalized inversion
The marginalized inversion described in Sect. 2 provides us with tolerance intervals of posterior fluxes, posterior correlations of errors and influence indicators. As the marginalized inversion filters out some data and regions, we present and analyse here the overall performance of our inversion, the effects of the data selection in the inversion and the implied limitations.

Temporal observational constraints
As explained in Sect. 2, the method developed by Berchet et al. (2015) filters out numerous observations and emission regions. Some observations available in the domain in 2010 are set aside before the inversion because of known flaws in CTMs. But the marginalized inversion also flags out additional observations when they are measured within plumes difficult to inverse. The remaining pieces of data do not all have the same weight in the inversion. Contrary to most classical inversion methods which cannot afford the computation of the explicit sensitivity matrix (see Sect. 2.3.1) informing in the weight of individual observations, the marginalized inversion allows us to explicitly analyse the use of the observations in the system. In Fig. 2, we represent the observations filtered out along the PBL height criterion before the inversion, the ones flagged out during the inversion and the relative weight (in degrees of freedom signal; dfs) of the remaining used observations. Many observations cannot be assimilated (black dots), especially in winter, when the very cold conditions (temperatures lower than −20 • C in average) related to the Siberian High generate very stable atmospheric conditions. In these conditions, the local emissions, which cannot be well assimilated in our inversion system because of the performances of the mesoscale transport model, significantly influence the observations. In addition to the numerous not-assimilated observations, all daily observations that are not filtered do not necessarily convey the same amount of information: all the blue circles in Fig. 2 depict pieces of data with a limited influence on the inversion (ω i,j < 0.6). The two main explanations for this inability to assimilate all the available pieces of data is the chosen scale of interest and the integrating character of the atmosphere. First, as we are interested in mesoscale fluxes, the system has been chosen with a spatial resolution of 25 to 100 km. All the variability in CH 4 mixing ratios driven by single local plumes cannot be reproduced in the system. Second, the limitation of the atmospheric inversion comes from the fact that the atmosphere behaves as an integrator, hence attenuating some information in the atmo-5404 A. Berchet et al.: CH 4 flux in Eurasia spheric signal. Tracking back the atmospheric signal to the fluxes then has intrinsic limitations. This limited capability of the system drastically reduces the number of usable pieces of information. Out of the 127 000 hourly measurements available in 2010, the pre-processing (as defined in Sect. 3.1) only retains 2000 daily aggregates into the inversion system. The system then excludes some observations and, at the end, only 800 data points remains, with 460 pieces of information (i.e. the trace of the sensitivity matrix; see details in Cardinali et al., 2004) carried by the atmospheric signal. Many observations give redundant information in our specific inversion framework at the mesoscale. The observations that cannot be processed by the mesoscale marginalized inversion carry information about local emissions.

Network range of constraints
The temporal use of the available observations matters in the inversion, as much as the relative use of the different observation sites. Fig. 3 displays the location and the average weight of each observation site (the coloured background is described and discussed further in Sect. 4.2). It is divided into four panels in order to separate the use of information related to anthropogenic (top panels) and wetland (bottom panels) emissions. As the marginalized inversion raw results are also processed in order to detect the regions that are mis-separated from the boundary conditions (see Sect. 2.3.2), data from edge observation sites is noticeably less used than for central sites.
Comparing anthropogenic (top panel) and wetland (bottom panel) maps, we notice that the weights of the observations are smaller for anthropogenic hotspot emissions. As expected, the inversion experiences difficulties in constraining emission hotspots, compared to diffuse fluxes. Concerning the spatial distribution of observation weights, wetlandrelated constraints follow the heuristics that the closer the observation site is to the fluxes, the higher the constraints to the inversion system is. Anthropogenic-related constraints do not exhibit such a pattern. For instance, NOY, close to the main oil extraction fields, has a lower observational influence than BRZ, remote from hotspots. Looking at wetland-related influences, BRZ has a bigger influence than NOY, while, as for anthropogenic emissions, wetlands emit more intensely in the vicinity of NOY. Then, in our mesoscale system, a surface observation site must be not too close, but not too far, from an emission hotspot to optimally constrain it. There is no generic criterion for this optimal distance to the observation sites as it depends on the atmospheric transport and on the intensity of the hotspots.
Looking at the differences in the relative weights of the observation sites between the raw inversion results and the LBC-separated ones (not shown here), one can notice that the sites at the edge of the domain of interest are logically dedicated to constraining the LBCs. Even the relative weight of the observation sites surrounding the Siberian lowlands is significantly reduced. Additional observations away from the region of interest would be necessary to better constrain the LBCs.

Constrained regions
The spatial distribution of the observational constraints on the fluxes is calculated from the sensitivity matrix (see Sect. 2.3). The information in this matrix is convolved with the prior distribution of the fluxes to deduce the maps in the right column of Fig. 3. The spatial distribution of the constraints on the fluxes depends on the intensity of the emissions and the distance to the observation sites. In the left column of Fig. 3, the closer to the Siberian network, the higher the constrains, independently of the intensity of emissions. For example, the western part of the Russian Federation contains most of the anthropogenic emissions of the country (roughly 20 TgCH 4 yr −1 according to EDGAR FT2010). But the constraints are lower for this region than for the Siberian lowlands with smaller emissions (8 TgCH 4 yr −1 according to EDGAR FT2010). The post-processing excluding the regions that are mis-separated from the LBC by the inversion highlights this pattern. The observation sites within the denser part of the network seem to constrain emissions within a radius of roughly 500 km.
In the right column of Fig. 3, we display the average fluxes which are not considered as constrained by the inversion (as detailed in Sect. 2.3.1). Despite the limited range of the observation sites and the high number of filtered-out data points, wetlands in the Siberian lowlands are constrained in a way such that the remaining unconstrained fluxes are of the same order of magnitude as minor wetland emitting region (e.g. in the far eastern parts of Russian shores on the Arctic ocean and the Pacific ocean). As a consequence, with the existing network constraining major wetland areas, minor wetland regions now contribute equally to the uncertainties on Siberian CH 4 budget. This points at a needed extension of the monitoring network toward these minor wetland areas. For anthropogenic emissions, the constraints on oil and gas related emissions are still too low. Anthropogenic emissions in the Siberian lowlands are still to be inquired into to reduced uncertainties on Russian CH 4 emissions.
In Fig. 5 described in Sect. 2.3.3, we also explicitly compute the portion of constrained emissions per month. On average, the major part of the emissions is not constrained by the inversion. The maximum proportion of constrained emissions is reached in summer with 50-70 % of the emissions constrained. In contrast, in winter, only a small part of the emissions are constrained. The proportion is critical in January and November when only 0-9 % of the emissions are constrained. In Sect. 4.1.1, we noticed that most of the observations are flagged out in winter because of the very low boundary layer. Consequently, emissions in the winter months are not well constrained. Overall, the configuration of the network gives valuable insights on Siberian lowland emissions, but it is not entirely adapted to our objective of constraining the Siberian lowland CH 4 budget, even in summer when supply issues do not prevent the acquisition of observations and when the atmosphere is mixed enough for the CTM to accurately reproduce the transport patterns. Additional observation sites would be needed for a complete resolution of the regional fluxes.

Solved time and space resolution
In order to inquire about possible improvements in the regional inversion, we compute the typical temporal and spatial resolutions the inversion can solve. These scales are plotted in Fig. 4.
For each pixel of the domain of interest, we consider the groups the pixel belongs to at the different periods of the year (months for anthropogenic emissions, 10-day periods for wetland emissions). We then average the size of the selected groups along the year. We do the same for the duration of each period. The lowest common time step for anthropogenic and wetland emissions is month. Anthropogenic fluxes are thus consistently solved at a monthly scale. As wetland emissions can be grouped with anthropogenic ones, their solved time resolution can be increased from the 10day basis. The temporal resolution is then computed only for wetland emissions in Fig. 4.
As expected, the resolved spatial resolution is better close to the network. Thus, most of Siberian lowlands are constrained with a typical resolution below 2 × 10 6 km 2 . The best resolution in our system configuration is roughly 700 000 km 2 . As we chose aggregation patterns with a mesh of about 300 000 km 2 in Sect. 2.4.2, this confirms that our aggregation procedures are not too coarse.
However, the resolved spatial resolution suggests that numerous additional monitoring sites would be required to identify emission patterns in relation to hydrological and meteorological parameters as wetlands can react quickly and with high gradients to changes in the weather or in the water table depth. A high temporal resolution would be required in the inversion to link wetland emissions to these physical parameters. We see in Fig. 4 that most wetland fluxes are constrained with a temporal resolution of typically 2 weeks, which is too long to resolve quick changes in emissions. Wetlands along the Yenisei River, far from anthropogenic emissions are solved with a better temporal resolution than wetlands in central lowlands. Unfortunately, this encouraging resolution is compensated by very high posterior uncertainties in the Yenisei fluxes (not shown). Overall, wetland emissions are resolved at a temporal resolution that allows the detection of the seasonal cycle, but not sufficient for linking wetland emissions to physical parameters varying at the synoptic scale.
To summarize, the inversion approach that we developed allows a precise quantification of the use of the observations. We also can deduce where the inversion results are the most reliable from the spatial influence of the network. At the mesoscale and in the Siberian framework, it appears that (i) hourly and even daily measurements are difficult to assimilate, (ii) anthropogenic emission hotspots require observation sites remote from them to be inverted, (iii) diffuse wetland emissions can be constrained with sites located close to them, (iv) the observation sites constrain fluxes within a radius of ∼ 500 km in our mesoscale inversion. This knowledge could find applications in network design and in the choice of the type of measurements.

Results of the marginalized inversion
In Sect. 5.1, we describe the methane fluxes resulting from the inversion system that are explicitly constrained and not mis-separated from the lateral boundary conditions. We then extrapolate the fluxes on the constrained regions to the Siberian lowlands and discuss the total budget of CH 4 of the region (Sect. 5.2).

A. Berchet et al.: CH 4 flux in Eurasia
Overall, on a yearly basis, posterior anthropogenic and wetland emissions are roughly equal to the prior. The tolerance interval is 1-13 TgCH 4 yr −1 for wetland emissions, and 6-16 TgCH 4 yr −1 for anthropogenic emissions. The ranges of uncertainties are reduced by 40 % for wetland emissions and 57 % for anthropogenic emissions from prior to posterior fluxes. Siberian oil and gas extracting activity and population centres are then responsible of 1.5-4.2 % of the global anthropogenic emissions. As a large portion of the oil and gas extracted in Siberia is exported to Europe, the inclusion of the emissions from the Siberian extraction process would have a significant impact on the total European emissions of ∼ 27 TgCH 4 (as computed from EDGAR v4.2 for all Europe apart from Russia, Ukraine and Belarus). The uncertainties of the wetland emissions are still too high to provide valuable insights for the modelling of these emissions, such as the start of the emitting season, or correlation to the precipitation rates or temperatures.

Wildfire influence
A peak of anthropogenic emissions with large uncertainties occurs in August. The large tolerance interval is still compatible with the prior scenario. Though, analysing into details posterior fluxes in August aggregation groups, we can propose two plausible explanation to this increase un-captured by prior data sets.
First, an increase of 0.7 TgCH 4 in August anthropogenic emissions is suggested in oil and gas extraction regions. This could be explained by the activities of prospecting companies. They make use of the decrease in the household demand in summer to carry out maintenance and welling operations on the oil and gas welling sites. These operations are related to punctual leaks and purging releases of gas.
Second, August is also a month with numerous forest fires. In particular, very large wildfires occurred in western Russia in August 2010, upwind our observation network. As said in Sect. 3.3, wildfires are systematically excluded from the marginalized inversion because they generate intense plumes difficult to take into account in our model. But, as western Russian wildfires took place far away of the observing sites, emitted plumes could have mixed before reaching these sites. Thus, the inversion system would have failed in eliminating observations very influenced by wildfires. The increase in the mixing ratios could then have been wrongly attributed by the inversion to an increase in anthropogenic emissions. As a confirmation, an aggregation group in August embraces anthropogenic emissions from western Russia, grouped with anthropogenic and wetland emissions in the Siberian lowlands. The inversion suggests that the emissions are increased by +100-200 % for this group, meaning roughly 0.5-1 TgCH 4 attributable to these intense wildfire episodes.

Performance on filtered-out data
For any inversion system, the inversion results are to be evaluated with independent data sets. A usual way to evaluate inversion results is to carry out a leave-one-out experiment. It usually consists in (1) flagging out all the observations from a site, (2) computing the inversion with the reduced data set, and (3) comparing the prior and posterior simulated mixing ratios to the observed ones at the left-out site. In principle, the leave-one-out inversion should improve the simulated mixing ratios at the left-out sites. That is to say, the posterior mixing ratios should be closer to and more consistent with the observations at the left-out site than the prior ones. In our case, the observation sites are far from each other and not numerous in regard to the size of the domain. Moreover, as developed in Sect. 4.2, the marginalized inversion explicitly informs about the constrained regions. Flagging out one site modifies the inversion results mainly in the surrounding regions where the uncertainty reduction becomes negligible compared to the complete data set. Therefore, the optimized fluxes in our leave-one-out experiments (not shown) remain within the range of uncertainty of another one for constrained regions. This confirms that the method we use consistently accounts for the uncertainties, but is not sufficient to quantitatively evaluate the optimized fluxes.
We can also use the data points which have been filtered out by our system in order to evaluate the inversion results. As the number of filtered-out observations is high, sampling biases may be expected from the filtering procedures. In addition, as only a few data points are assimilated, unrealistic fluxes could have been inferred by the inversion to fit the assimilated data leading to a flawed reproduction of the remaining observations. As shown in Table 2, the marginalized inversion significantly improves the simulated mixing ratios at the sites where data are used as expected. The model results are also well improved for unused data, meaning realistic flux prescription despite the applied filters.
This confirms that our method does not create sampling biases despite the high number of filtered-out data points. It also confirms that the increments on the fluxes are realistic from the point of view of our network.

GOSAT evaluation
Genuinely independent observations are required for a better evaluation of the inversion results. Long-term monitoring surface sites are scarce in Siberia, and airborne (e.g. Paris et al., 2010;Berchet et al., 2013a) or train (e.g. Tarasova et al., 2006) measurement campaigns only provide snapshots of the Siberian atmospheric composition. Therefore, we try to evaluate the marginalized inversion results with satellite data. We choose GOSAT total column biased-corrected retrievals (see Sect. 3.5) and compare them to their simulated into the background. This fast dilution can explain at least part of the factor of difference of 2-3 between observed and simulated CH 4 columns.

Toward using satellite measurements in regional frameworks
Satellite data, as a tool for compensating for the lack of observations in Siberia to evaluate the inversion results, do not seem to be suited for our regional configuration. As satellite observations have very good potential in term of spatial coverage, in situ observation sites, which are difficult to maintain and to extend into the Siberian framework, should be complemented by new satellite data sets.
The inversion system could be developed in order to use GOSAT data as a proxy of the large scale gradients given here by a global model. With such additional observations, the LBCs, representative of the emissions outside the domain of interest, would be better constrained in the inversion. As a consequence, fewer regions would be expected to be misseparated from the LBCs, and the assimilated surface sites at the edge of the network would then provide more information about the emissions within the domain of interest. In addition, high resolution transport simulations close to the hotspots would better represent the dilution of the plumes. The model-observation comparison for quantifying anthropogenic sources would then be more suitable. However, the plumes generated by the hotspots are not necessarily simulated at the right location and date. Indeed, the location of emission hotspots is not always well known a priori and temporal or spatial transport mismatches can occur. Direct comparisons between the observed and simulated CH 4 columns would face the same issue as for surface sites when assimilating hotspots plumes. Recent developments (e.g. Krings et al., 2013) point to the possibility of using 1-D or 2-D high resolution snapshots of the hotspot plumes to infer information about very local emissions. Integrated comparison of the observed and simulated plumes in the CH 4 columns could then be implemented in a mesoscale inversion system.
With such techniques, future satellite missions with active remote sensing (e.g. the joint French-German cooperation Methane Remote Sensing LIDAR Mission, MERLIN) providing high resolution accurate 1-D or 2-D products could be used in regional inversions; the spatial resolution of the products to be used in such an inversion system should of the same magnitude as the mesoscale transport model, i.e. at least < 50 km.

Conclusions
We assimilated the data collected in 2010 at eight surface observation sites measuring atmospheric CH 4 mixing ratios in the West Siberian Plain into a regional atmospheric inversion. It was the first time all these observations were used in a single study. As regional inversions suffer from mis-specified uncertainties, we implemented an enhanced Bayesian method developed by Berchet et al. (2015) in order to get reliable results at the regional scale with an objectified quantification of the uncertainties in the system. This new method allows us to consistently evaluate the local spatial distribution of the sensitivity of the emission areas to the inversion and the usefulness of each available observation. The inversion seems to be able to primarily constrain the emissions close to the observation sites (within a radius of roughly 500 km). The inversion system assimilates daily observation aggregates to constrain the emissions. Amongst all the observation aggregates, and despite the efforts to provide precise and quasi-continuous measurements, our mesoscale inversion system properly uses only one piece of information every few days. This is mainly caused by atmospheric limitations related to transport and mixing; the fewer the assimilated data, the higher the uncertainties after the inversion. Even in the regions close to the stations, the posterior uncertainties (objectively quantified) thus remain larger than 20 % of the prior fluxes for anthropogenic emissions and 50 % for wetlands, although an important error reduction is achieved. Nevertheless, objectified uncertainties allow a robust evaluation of the wide range of proposed wetland and anthropogenic emissions in Siberia. On average, the posterior tolerance interval (defined so that 68.27 % of the Monte Carlo marginalized ensemble is within the interval) on the West Siberian Plain methane budget is 5-28 TgCH 4 for the year 2010.
The year 2010 was the first year when most of the used observation sites were functional. Reproducing our set-up to subsequent years would provide a more robust estimation of the regional fluxes and possibly valuable information about the year-to-year variability of Siberian methane fluxes. Finally, satellite platforms provide an extensive spatial coverage of observational constraints. Implementing such rather uniform observation coverage in a regional framework with few surface sites is tempting. However, with the inversion framework used here, satellite data would be useful only for constraining large scale gradients, hence the lateral boundary conditions. Further work on inversion systems is required so that satellite observations can be used to quantify local emissions in a regional framework like this one.