Interactive comment on “ A multi-year methane inversion using SCIAMACHY , accounting for systematic errors using TCCON measurements ” by S .

Abstract. This study investigates the use of total column CH4 (XCH4) retrievals from the SCIAMACHY satellite instrument for quantifying large-scale emissions of methane. A unique data set from SCIAMACHY is available spanning almost a decade of measurements, covering a period when the global CH4 growth rate showed a marked transition from stable to increasing mixing ratios. The TM5 4DVAR inverse modelling system has been used to infer CH4 emissions from a combination of satellite and surface measurements for the period 2003–2010. In contrast to earlier inverse modelling studies, the SCIAMACHY retrievals have been corrected for systematic errors using the TCCON network of ground-based Fourier transform spectrometers. The aim is to further investigate the role of bias correction of satellite data in inversions. Methods for bias correction are discussed, and the sensitivity of the optimized emissions to alternative bias correction functions is quantified. It is found that the use of SCIAMACHY retrievals in TM5 4DVAR increases the estimated inter-annual variability of large-scale fluxes by 22% compared with the use of only surface observations. The difference in global methane emissions between 2-year periods before and after July 2006 is estimated at 27–35 Tg yr−1. The use of SCIAMACHY retrievals causes a shift in the emissions from the extra-tropics to the tropics of 50 ± 25 Tg yr−1. The large uncertainty in this value arises from the uncertainty in the bias correction functions. Using measurements from the HIPPO and BARCA aircraft campaigns, we show that systematic errors in the SCIAMACHY measurements are a main factor limiting the performance of the inversions. To further constrain tropical emissions of methane using current and future satellite missions, extended validation capabilities in the tropics are of critical importance.


Introduction
To reduce the uncertainty of climate change prediction, improved understanding of the global cycles of long-lived greenhouse gases is needed. This is true in particular for methane, the second most important anthropogenic greenhouse gas, which has shown important growth rate variations in the past 2 decades that are only poorly understood Kirschke et al., 2013;Aydin et al., 2011;Pison et al., 2013). To study the drivers of changes in the methane growth rate requires tools for quantifying its sources and sinks at large scales, which remains a scientifically challenging task.
Inverse modelling techniques make use of independent information from atmospheric measurements, to improve the available emission estimates and to reduce their uncertainties. In the past, these techniques have most commonly been applied to ground-based measurements from global monitoring networks (Bousquet et al., 2006Houweling et al., 1999;Hein et al., 1997). In recent years, measurements from satellites have become available, such as from SCIAMACHY (Bovensmann et al., 1999;Frankenberg et al., 2005Frankenberg et al., , 2008aFrankenberg et al., , 2011Schneising et al., 2011) and more recently also from GOSAT (Yokota et al., 2009;Kuze et al., 2009;Butz et al., 2011;Schepers et al., 2012), and significant efforts have been spent to investigate their use in inverse modelling Bergamaschi et al., 2009Bergamaschi et al., , 2013Fraser et al., 2012).
Satellite instruments offer spectacular improvements in measurement coverage, particularly over tropical continents, which have been the Achilles heel of the surface monitoring networks for many years. Using satellite retrievals, several studies have addressed the size of the tropical methane flux, pointing to previously underestimated emissions (Frankenberg et al., 2005Bergamaschi et al., 2009;Beck et al., 2012;Bergamaschi et al., 2013) that were corrected downward in part after revisions of the methane and water vapour spectroscopy (Frankenberg et al., 2008a, b). Despite these efforts, however, inversion-derived tropical methane emissions remain highly uncertain due to the methodological difficulty of accounting for poorly characterized measurement and transport model uncertainties.
The use of SCIAMACHY retrievals in inverse modelling requires a bias correction algorithm (Bergamaschi et al., 2009. The common approach is to account for inconsistencies between surface and satellite measurements using a set of bias correction functions that are optimized along with the methane sources and sinks. Here, surface measurements are used in addition to the satellite data to provide additional constraints on the surface fluxes and act as anchor points to constrain the parameters of the bias correction functions. However, it is not clear how to best define the bias correction functions and to account for their uncertainties. The problem is further complicated by transport model uncertainties and model-data representation errors, which have systematic error components also. The approach taken by Meirink et al. (2008b) and Bergamaschi et al. (2013) is to optimize a quadratic function of latitude for each month of available data, to account for inconsistencies between modelled and observed latitudinal and seasonal variations in the total column-averaged methane mixing ratio (XCH 4 ).
For inverse modelling applications spanning many years of data this approach has two important limitations. First, since the bias functions have no direct relationship with known shortcomings of the measurements, it is unclear to which extent the optimized parameters account for systematic error in the measurements, or systematic model error, or real signals in the satellite data. Secondly, the number of bias parameters increases proportionally with the number of years, introducing additional degrees of freedom capable of absorbing much of the information on inter-annual variation provided by the measurements.
The purpose of this study is to further investigate the role of measurement bias in methane inversions using satellite data, to assess how sensitive studies such as Bergamaschi et al. (2013) are to the bias correction method that is used. For this purpose, an alternative approach to bias correction is introduced, which makes use of measurements from the Total Carbon Column Observing Network (TCCON) (Wunch et al., 2011a) to identify systematic errors in satellite retrievals. Several versions are tested to study the sensitivity of the optimized fluxes to the bias correction method, including a direct comparison with results from Bergamaschi et al. (2013). Further, we test whether comparisons with independent measurements allow evaluation of the bias correction approach.
Inverse modelling calculations are carried out using TM5 4DVAR (Meirink et al., 2008b;Bergamaschi et al., 2013) for the period 2003-2010, with and without the use of SCIA-MACHY retrievals . This period covers the transition from stationary to increasing methane mixing ratios, and therefore allows us to contribute to the ongoing discussion about the underlying causes. Optimized fluxes derived using only surface measurements act as a reference for assessing the contribution from SCIAMACHY in joint inversions of surface and satellite data.
Section 2.1 describes the various elements of our inversion setup, including a priori assumptions, and the bias correction functions that are used. In Sects. 3.1 and 3.2, we analyse the statistical consistency of the inversion-derived fits to the data and investigate the performance of the optimized TM5 model against independent measurements. In Sect. 3.3 we analyse the posterior fluxes, and assess the sensitivity of the inversion-derived flux adjustments to the bias correction algorithm. Section 4 discusses the implications for our understanding of the recent evolution of the global methane cycle.
Here it is used for estimating large-scale emissions of CH 4 , by optimizing the agreement between measured and model simulated mixing ratios. This is achieved by minimizing a Bayesian cost function J(x), defined as with respect to a state vector of surface fluxes x. H is a Jacobian matrix of the sensitivity of the observed mixing ratios to the surface fluxes, calculated using the atmospheric transport model TM5 (Krol et al., 2005). The misfits to the measurements y and a priori fluxes x 0 are weighted by their respective uncertainties (R and B).
Equation (1) is solved using the variational approach described by Meirink et al. (2008b). In short, an iterative procedure is followed of alternating forward and adjoint transport model simulations. The forward model is used to calculate model-data mismatches (y − Hx) for each successive update of the flux vector x. The adjoint model (H T ) uses the modeldata mismatch vector to evaluate the cost function gradient (∇J(x)), represented by The trial state vector x and corresponding cost function gradient are input to a Lanczos minimization algorithm (Lanczos, 1950), yielding an updated state vector. The procedure is terminated after a fixed number of iterations or after the remaining gradient passes a preset convergence criterion. The inversion results presented in this article are obtained after 55 iterations, which corresponds to a gradient norm reduction of about a factor 5000. TM5 simulations were performed for the period 2003-2010 using meteorological fields from the ECMWF ERA Interim reanalysis project (Dee et al., 2011), preprocessed for use in TM5 at 4 • × 6 • (latitude × longitude) and 25 hybrid sigma pressure levels from the surface to the top of the atmosphere. The inversion makes use of the "cy1" version of TM5-4DVAR, extended with a parametrization of horizontal diffusion using the modified version of Prather et al. (1987) described in Monteil et al. (2013). This scheme has been introduced to account for the underestimated inter-hemispheric mixing rate in TM5 reported in Patra et al. (2011). The tuning parameter for the horizontal diffusion strength (D in Prather et al., 1987) is adjusted such that the TM5 simulated north-south gradient for the year 2003 is brought in agreement with available measurements of SF 6 (see Monteil et al., 2013).
The inversion solves for the monthly mean CH 4 flux into each surface grid box of the model. Note that, in contrast to Bergamaschi et al. (2013), no distinction is made between contributions from different processes, because the measurements provide only limited process-specific information at the resolution of the transport model. Also, they provide limited information to distinguish between surface fluxes and atmospheric sinks. For this reason we choose to solve for surface fluxes only, which are thus optimized given the assumed atmospheric sinks.

A priori fluxes
Anthropogenic emissions from fossil fuel use, agricultural practices and waste treatment are taken from the EDGAR4.1 emission inventory (European Commission, Joint Research Centre (JRC)/Netherlands Environmental Assessment Agency (PBL), 2010), which provides annual emissions at 0. by the fractional increase of global fuel production and consumption since 2005. Our extrapolation procedure differentiates between different fuel types (oil, coal and gas), and sources classified as "energy production". For the remaining fossil sources, such as the EDGAR category "residential commercial and others", we choose the BP reported production/consumption categories that best reproduce the EDGAR emissions for the period until 2005.
For agricultural sources, we apply the same procedure using FAO reported indices for gross production, rice production and cattle number. Remaining categories, for example waste treatment, are extrapolated using a linear regression of the global EDGAR4.1 emissions for the period 2000-2005. This procedure yields an increase of anthropogenic CH 4 emissions between the years 2003 and 2010 of 37 Tg CH 4 yr −1 . EDGAR4.1 does not provide information on intra-annual emission variations. In TM5 we neglect seasonal variations in fossil and agricultural sources, except emissions from rice agriculture, which are distributed seasonally following the rice cropping calendar of Matthews and Fung (1987).
For emissions from natural wetlands we use the output from a transient simulation of the process-based model LPJ-WhyMe (Spahni et al., 2011), without the correction derived from atmospheric measurements, covering our simulated period until the end of 2008. For the remaining years we use LPJ-WhyMe emissions averaged over the period 2003-2008. These emissions include emissions from wet soils, and also account for microbial soil oxidation under dry conditions. In regions with rice agriculture, LPJ-WhyMe emissions are not used, to avoid double counting. Biomass burning emissions are taken from GFED3.1 , except for emissions from agricultural waste burning, which are included in EDGAR4.1.
Geologic sources of methane have received much attention in recent years (Etiope and Klusman, 2002;Etiope and Milkov, 2004;Kvenvolden and Rogers, 2005), with reported global emissions ranging between 13 and 80 Tg CH 4 yr −1 (Judd, 2004;Etiope et al., 2008). Their representation in global models is complicated by the limited number of mud volcanoes for which the emissions have been quantified. The available estimates vary widely and represent only a small fraction of the extrapolated global totals. Digital global maps that would be needed to construct gridded global emissions for use in models do not exist or are not made available for commercial reasons. This is true in particular for micro seepage from sedimentary basins on land, with an estimated global contribution exceeding 10 Tg CH 4 yr −1 and covering 10 % of the global drylands (Etiope and Klusman, 2010). For this reason, we account only for 1300 mud volcanoes and seeps that are reported worldwide and are archived in the GLOGOS database (http://www.searchanddiscovery. com/documents/2009/090806etiope/). To avoid biasing our emissions to extensively studied sites, we distributed them evenly over all reported coordinates.
For marine seepage, the situation is similar. Despite coordinated activities to archive all available information (Bange et al., 2009), no global inventory is available yet. Our approach is to distribute a conservative estimate of global emissions from marine seeps and hydrates of 17 Tg yr −1 uniformly over the continental shelves (see Table 1). This simplified approach does not account for what is known from ongoing regional efforts to quantify, for example, the hot spot of geologic emissions in Azerbaijan (Etiope and Klusman, 2002) or the reported marine seepages in the Laptev Sea (Shakhova et al., 2010). A priori flux uncertainties are set to 50 % of the combined monthly flux in each grid box of the model, with uniform and exponentially decaying spatial and temporal correlation scales amounting to, respectively, 1000 km and 3 months.
Tropospheric methane oxidation is represented using methyl chloroform optimized OH distributions (Huijnen et al., 2010). A single scaling factor of 0.92 is applied to the seasonal varying global OH field of Spivakovski et al. (2000), which is recycled each year of simulation, hence ignoring inter-annual variability. We do not account for the estimated contribution of marine Cl radicals to the tropospheric CH 4 sink (Allan et al., 2001(Allan et al., , 2005Lassey et al., 2011). Stratospheric CH 4 lifetimes are derived from the Cambridge 2-D model (Law and Pyle, 1993), and account for CH 4 oxidation by OH, O( 1 D), and Cl radicals.
Like most off-line transport models, TM5 tends to underestimate the stratospheric age of air (Jones et al., 2001). To account for the generally overestimated contribution of upper stratospheric CH 4 to XCH 4 in TM5, we correct model sampled vertical profiles above 50 hPa using a CH 4 climatology based on HALOE/CLAES observations (Randel et al., 1998). This procedure includes a linear correction for the CH 4 increase since the period of observation. Table 1 summarizes the prescribed a priori fluxes, including their corresponding global emissions for the year 2005.

Measurements
Inverse model calculations were carried out using either surface measurements (referred to as NOAA-only) or a combination of surface measurements and SCIAMACHY satellite retrievals (referred to as SCIAMACHY). The surface measurements are taken from selected sampling sites of NOAA-ESRL's global cooperative air sampling network . Sites were selected (see Fig. 1) that had been operational for the full period of our simulations, with limited data gaps. The measurements, as used in the inversion, represent averages of samples that fall in the same model grid box in 3-hourly time intervals.   The data covariance matrix accounts for measurement uncertainties and model representation errors. We use the local concentration variability as a proxy for model representation error (also referred to as model data mismatch error). The assigned value is the maximum of: (i) the standard deviation of the samples that are averaged, and (ii) the standard deviation of the simulated CH 4 mixing ratios in the surrounding grid boxes. The measurement uncertainty is set to 2.4 ppb for surface measurements, which represents the average pair difference of the NOAA flask measurements. All errors are assumed independent.
We use SCIAMACHY retrievals from the IMAPv5.5 proxy retrieval reported by Frankenberg et al. (2011). The retrieved CH 4 / CO 2 ratios have been converted into XCH 4 using CarbonTracker-derived XCO 2 (Peters et al., 2007). Optimized CarbonTracker fields were unavailable for the last 3 years of our simulation. The CarbonTracker estimates for 2007 were extrapolated for those years using the global growth rate of CO 2 , estimated by averaging the surface measurements from South Pole, Mauna Loa, and Alert. The use of SCIAMACHY retrievals is limited to 50 • N to 50 • S.
Representation errors for SCIAMACHY are calculated from the modelled local variability in XCH 4 , using the simulated columns in the surrounding grid boxes. For the measurement uncertainty we take the maximum of the retrieval uncertainty and the standard deviation of the averaged retrievals. The spatio-temporal interval for averaging is the same as for the surface measurements. The uncertainties of the averaged SCIAMACHY retrievals are assumed to be uncorrelated. To account for correlated retrieval error of nearby retrievals, we treat the errors of the retrievals within the 6 • × 4 • grid boxes and 3-hourly intervals as fully correlated. This approach yields 1σ uncertainties in the range of 15 to 45 ppb. In addition, we account for measurement bias as described in the next section.
The SCIAMACHY bias correction, described in detail in the next section, makes use of TCCON data from the GGG2012 release. Figure 1 shows the sites that were used. The estimated 1σ single measurement uncertainty for CH 4 , after correction for systematic errors using aircraft profile measurements, amounts to 3.5 ppb (Wunch et al., , 2011a.

Bias correction
A common approach to account for systematic measurement errors in inversions is to define a set of functions describing the spatio-temporal dependence of these errors, with uncertain parameters that are optimized in the inversion along with the other elements of the state vector. Following this approach, the comparison between modelled and observed mixing ratios can be written as where and f b (α i ) represent, respectively, the random and systematic component of the data uncertainty. The systematic error components f b (α i ) are functions of the added state vector elements α i . As for the random error , the systematic error has contributions from both the measurements and the model. In practice, well-quantified biases are usually directly corrected in the model or the measurements. Therefore, the bias correction we are dealing with here targets ill-defined residual biases. In our experience, the optimization of bias terms in the inversion can lead to solutions that are out of their a priori assigned uncertainty ranges. This is a sign that a bias function ends up accounting for other uncertain elements in the inversion, such as unaccounted model errors or even the sources and sinks, which are meant to be estimated by the inversion. In that case, it may be preferable to prescribe the initial bias correction, without further optimization. In this study, we perform inversions with and without optimization of bias corrections to investigate the impact. Different strategies can be followed for defining f b (α i ), using 1. empirical functions of space and time, 2. information about known retrieval uncertainties, 3. differences between retrievals and independent measurements.
An example of strategy 1 is the use of seasonal and latitudinally varying functions in Bergamaschi et al. (2010Bergamaschi et al. ( , 2013, based on the difference between SCIAMACHY and TM5 after optimization using surface measurements. The limitation of this approach is that some unknown portion of the correction may reflect real signal, which is accounted for as bias and therefore cannot improve the model anymore. In our view, a preferable approach is to combine strategies 2 and 3, by making use of independent measurements to support the quantification of known terms in the retrieval error budget. A good example is the study by Wunch et al. (2011b) on the use of TCCON data for evaluating systematic errors in ACOS-GOSAT XCO 2 retrievals. For SCIAMACHY XCH 4 , we constructed a model of known retrieval uncertainties on the basis of Butz et al. (2010), accounting for spectroscopic errors varying with the sampled air mass and residual aerosol errors varying with the difference in surface albedo between the weak short-wave infra-red (SWIR) CO 2 and CH 4 absorption bands. Optimization of the bias coefficients in the inversion, however, turned out to be rather inefficient in correcting apparent inconsistencies in the optimized fits to surface and total column measurements (not shown). Analysis of the difference between co-located SCIAMACHY and TCCON measurements pointed to a seasonal pattern that correlates well with climatic variables such as temperature and humidity. These parameters lag variations in solar zenith angle (SZA), due to the inertia in the earth's climate system, which may explain why the use of air mass for seasonal bias correction did not work well. Figure 2 demonstrates the relation between the SCIA-MACHY seasonal bias with respect to TCCON and ERAinterim derived specific humidity. Vertical profiles of specific humidity are used that have also been used in the IMAP v5.5 retrieval. The values represent averages over the lowest 3 km of the vertical profile (up to ∼ 700 hPa over sea). Influences of water vapour on the SCIAMACHY CH 4 retrieval were discussed by Frankenberg et al. (2008a), who corrected the water vapour spectroscopy in the 1.6 micron CH 4 absorption band. It is not clear whether water vapour is the cause of the seasonal bias discussed here, or that it only happens to covary with a different underlying cause. In either case, the correlation of the SCIAMACHY-TCCON residuals with water vapour offers a useful proxy for SCIAMACHY bias correction.
This finding motivated us to define a bias correction using a limited set of functions accounting for the main differences between SCIAMACHY and TCCON (i.e. a combination of strategies 1 and 3). Besides seasonality, our bias correction method accounts for a global offset and a latitudinal correction. Note that the humidity proxy already introduces a latitudinal correction (as well as corrections on other timescales, including that of inter-annual variability). However, differences with TCCON remain, which are corrected separately. Unfortunately, the TCCON network is not dense enough to derive a latitudinal correction function. To extend coverage, we make use of GOSAT proxy XCH 4 retrievals reported by Schepers et al. (2012), which are significantly more accurate than SCIAMACHY retrievals (see Schepers et al., 2012;Monteil et al., 2013). This is confirmed by Fig. 3, which shows that residuals between TCCON and seasonality corrected SCIAMACHY show a similar latitudinal variation as residuals between SCIAMACHY and GOSAT. To isolate the large-scale component of this relationship, the SCIA- MACHY minus GOSAT residuals have been smoothed using a boxcar filter (width = 20 • latitude).
In summary, the bias correction functions that we will investigate are the following: -S: a global uniform scaling factor only.
-SQ: as "S", extended with the seasonal correction using humidity.
The corresponding bias coefficients α are determined by linear regression to the SCIAMACHY-TCCON residuals, using TCCON measurements from the sites shown in Fig. 1 for the years 2009 and 2010. To obtain a single set of coefficients for use in S, SQ and SQG, the regressions were carried out in this order. Each extension was regressed to the TCCON-SCIAMACHY residuals up to that extension. To avoid interference with S, functions Q and G were normalized by taking out their means. Table 2 summarizes general statistics of the difference between TCCON and co-located SCIAMACHY retrievals for the different stages of bias correction. As can be seen in Table 2, the uniform scaling factor decreases the rms difference between TCCON and SCIAMACHY by about a factor of two. Smaller additional improvements are obtained after introducing the seasonal and latitudinal corrections. A large fraction of the remaining rms is explained by scatter in the retrievals on relatively small scales.

Experimental setup and scenarios
Inverse modelling calculations have been performed for the period 2003-2010. For efficient calculation of inversions covering 8 years, the target period is split into three blocks of 4 years each, overlapping the adjacent inversion(s) by 2 years. The calculations for these blocks are carried out in parallel. The initial concentration field was optimized following the approach of Meirink et al. (2008b).      discontinuities at the inversion boundaries. For the SCIA-MACHY inversions larger differences were found initially due to trade-offs between the optimization of bias parameters and the initial condition. This problem has been circumvented by using the NOAA-only optimized initial con-ditions as prior for the SCIAMACHY inversions, with the corresponding prior uncertainties reduced by a factor of 10. Table 3 summarizes the inverse modelling calculations that were performed. The SCIAMACHY inversions consist of two sets: one with and the other without optimization of bias parameters (referred to as "flex" and "fix" respectively). In the "flex" inversions, the bias coefficients of the corresponding "fix" inversions are taken as first guess. The a priori uncertainty of the bias coefficients has somewhat arbitrarily been set to 100 %. This is a conservative estimate given the statistics of the TCCON corrections. However, an unquantified additional uncertainty comes from the assumption that the TCCON comparisons are representative for the full global domain. The bias coefficients for the seasonal and latitudinal bias functions represent fractional adjustments of the TCCON-derived a priori corrections (referred to as "normalized corrections" in Table 3). For the humidity correction, the TCCON-derived regression coefficient amounts to 3.42 ppb (g kg −1 ) −1 (which corresponds to α 2 = 1 in the table). As can be seen in Table 3, the bias optimization leads to significant corrections of the global scaling factor (α 1 ) and the latitudinal correction (α 3 ), whereas the seasonal correction (α 2 ) remains close to the TCCON-derived estimate.
Since the 4-year blocks are computed independently, the optimized bias parameters vary between the blocks. It is possible to pre-optimize the bias parameters to enforce consistency between the blocks. This is not done in this study, because the differences between the blocks may provide valuable information, for example about instrumental drift (see Sect. 3.3 for further discussion).

Posterior concentrations and bias optimization
Before we analyse the statistics of inversion fit residuals and comparisons with independent data, we first verify that the inversion results support the bias correction functions derived in the previous section.

Seasonal cycle
To test the seasonal cycle correction, the inversion-derived posterior CH 4 mixing ratios are sampled at the TCCON sites, using averaging kernels that were made available for the FTS at Lamont and site-specific seasonally varying a priori profiles. Figure 4 shows average seasonal cycles for northern and southern hemispheric TCCON sites. Averages were taken by combining samples for the period 2009-2010, and smoothing the data by boxcar filtering. For the Southern Hemisphere only two sites are available (Lauder and Wollongong), which do not provide sufficient data. This problem was solved by extending the averaging period to 2007-2010 and combining the Bruker 120 and 125 data from Lauder. For the Northern Hemisphere six sites were selected with good coverage during 2009-2010, representing North America (Park Falls, Lamont), Europe (Bialystok, Garmisch, Orleans), and Japan (Tsukuba). Uneven data coverage between years and sites, in combination with a positive trend in this period, may influence the shape of the derived seasonal cycles. However, since all simulations are sampled in the same way as the data, the comparison is still valid. The results in Fig. 4 for the northern hemispheric TC-CON sites clearly show a phase-shifted seasonal cycle in the uncorrected SCIAMACHY inversion (S fix ). In contrast, the seasonality of the NOAA-only inversion agrees closely with TCCON. This is explained in part by the constraints of the surface measurements on the seasonal cycle in the troposphere. However, the seasonal dynamics of tropopause height contributes significantly to the seasonality of total column CH 4 also (see e.g. Warneke et al., 2006;Washenfelder et al., 2003). The close agreement between the NOAA-only simulated and TCCON observed seasonal cycle therefore supports the TM5 simulated seasonality of the stratospheric partial column.
The results for the seasonally corrected SCIAMACHY inversion (SQ fix ) confirm that the humidity proxy largely accounts for the seasonal misfit of the uncorrected IMAPv5.5 retrievals. In the Southern Hemisphere the seasonal misfit between the uncorrected SCIAMACHY data and TCCON is much smaller, in agreement with a smaller seasonal variation of specific humidity. This means that the humidity proxy accounts for a north-south asymmetry in the seasonal bias of the SCIAMACHY retrieval.

SCIAMACHY vs. NOAA-only
Next we investigate whether the SCIAMACHY bias corrections improve the consistency of the NOAA-only optimized model. Figure 5 shows maps of the residual differences between the two, averaged seasonally. These differences can be interpreted as the signal that is brought in by SCIAMACHY, in addition to what is provided by the surface network. At temperate latitudes the residuals are on average −30 ppb, roughly consistent with the 2 % underestimation of the SCIAMACHY retrievals at the TCCON sites. Over tropical Africa and tropical South America, however, the residuals change sign with regional averages up to 18 ppb. As can be seen in the other panels of Fig. 5, this difference between the tropics and extra-tropics is corrected to a large extent, albeit not fully, by the seasonal bias correction. This is because the humidity proxy varies not only temporally, but also spatially with a significant latitudinal component. The SQG fix bias correction leads to a systematic overestimation of XCH 4 compared to the NOAA-only inversion (13 ppb on average), which is largely reduced if the scaling factors are further optimized in the SQG flex inversion (−5 ppb on average). This result can be explained by the general tendency of the NOAA-only inversion to slightly underestimate TCCON. Note that the differences that remain after bias correction in Fig. 5 include real signals of methane sources that are not represented by the NOAA-only inversion.

Extended comparison to TCCON
The SCIAMACHY calibration using TCCON mainly represents the period after the launch of GOSAT (2009GOSAT ( -2010, when most sites became operational. It raises the question of whether the inferred bias corrections are valid also for other years. With only a limited number of available TCCON stations for years prior to 2009, this is difficult to assess. Nevertheless, an attempt is made as shown in Fig. 6. Reasonable agreements are found, except that the "flex" inversions show systematic underestimation of the TCCON measurements, which is seen most prominently at Darwin. The NOAAonly inversion underestimates TCCON also, with systematic differences increasing towards southern latitudes (3 ppb for Park Falls, 10 ppb for Lauder). This finding is in line with Monteil et al. (2013), who point to inaccuracies in the TM5 simulated atmospheric transport or the chemical oxidation of methane, which are not optimized in the current inversion setup. It should be realized that the aircraft measurements for calibrating the TCCON network cover only a small fraction of the presented time series, and therefore the calibration of the TCCON-observed inter-annual variability has its limitations also. The comparison at Park Falls shows clear examples of the phase-shifted seasonality of the S fix corrected SCIAMACHY data, in particular during the first months of 2005 and 2006. Overall, the SQ fix inversion shows the best agreement with the three long-term TCCON time series (rms = 9.0 ppb), followed by SQG fix (rms = 9.6 ppb).

Statistics of fit residuals
Next, we analyse the distributions of residuals between the optimized models and the NOAA measurements that were used in the inversions (see Fig. 7). Deviations from a Gaussian distribution centred around the origin may point to remaining inconsistencies between the bias-corrected SCIA-MACHY data and the surface measurements. As can be seen, the a priori model systematically overestimates the measurements. The NOAA-only inversion shows on average the smallest residuals. This result is expected, since this inversion fits the surface measurements, without constraints to also fit satellite retrievals. SCIAMACHY inversions with fixed biases lead to systematic overestimation of the NOAA flask measurements, which is effectively corrected in the "flex" inversions by reducing the uniform bias component (see Table 3). Of the alternative SCIAMACHY bias corrections that were tested, SQG flex and SQ flex show the smallest residuals, with little difference in performance between these two.

Performance evaluation using independent data
Aircraft measurements are particularly useful for evaluating inversion optimized concentration fields, since they address the modelled vertical profile of CH 4 , which connects the observational constraints on the surface and the total column that are used in the inversions. In all comparisons between model and measurements the model was interpolated to the coordinates and times of single measurements. Figure 8 shows comparisons between inversion optimized mixing ratios and measurements from the HIPPO aircraft campaigns (Wofsy, 2011;Kort et al., 2012), representing north-south cross sections across the Pacific Ocean covering altitudes between the planetary boundary layer and the tropopause. The HIPPO measurements were taken from campaigns 1-3, which took place between January 2009 and April 2010. The comparisons between model and measurements for these campaigns are combined in one figure. The HIPPO measurements have been linked to the NOAA scale using coincident measurements from discrete samples collected during the flights and analysed at NOAA. The remote Southern Hemisphere offers a good opportunity to verify this approach, since the concentration gradients in the marine boundary layer are small. Therefore a good agreement is expected between the NOAA optimized model and the HIPPO measurements. Indeed, at latitudes southward of about 30 • S differences are within 3 ppb. The differences between S fix and HIPPO show the expected pattern of over-and underestimation of respectively tropical and extra-tropical CH 4 in absence of a latitudinal varying SCIAMACHY bias correction. When such corrections are applied the latitudinal differences reduce, although minor differences remain (within ∼ 10 ppb). The inversions with the fixed, TCCON-derived, bias corrections tend to overestimate the HIPPO measurements by 6-12 ppb. Further optimization of the bias coefficients in the inversion im-proves the agreement with HIPPO, with a low bias of close to 5 ppb. When comparing the results in Figs. 6 and 8 it becomes clear that surface and aircraft data support lower CH 4 mixing ratios than TCCON. This happens despite the fact that TCCON has been validated using aircraft measurements also (see e.g. Deutscher et al., 2010;Wunch et al., 2010;Geibel et al., 2012). The remaining differences could be explained by the stratospheric contribution to the total column, which is the subject of ongoing investigations (see Monteil et al., 2013). The residuals between HIPPO and TM5 point to overestimated mixing ratios in the lower stratosphere at higher latitudes (> 50 • ). The SCIAMACHY retrievals at those latitudes, however, are not used in the inversion.
The question remains, whether independent data support the residual differences between bias-corrected SCIA-MACHY and NOAA-optimized TM5 over tropical continents. To address this question, comparisons are presented between optimized TM5 and aircraft measurements from the BARCA campaigns (Beck et al., 2012) over the Amazon rain forest (see Fig. 9). Results are shown for campaigns in November 2008 (end of the dry season) and May 2009 (end of the wet season). All measurements are combined in single profiles, since the resolution of the model is not considered Atmos. Chem. Phys., 14, 3991-4012 adequate to resolve regional gradients within the Amazon basin. Overall, the results confirm the conclusions from Beck et al. (2012) that using of only surface data in the inversion leads to underestimated posterior mixing ratios. Including SCIAMACHY increases the mixing ratios over South America, to an extent which varies considerably between the different bias correction approaches. Using only a uniform scaling factor leads to overestimated mixing ratios. The SQ and SQG inversions are in closer agreement with the measurements. The relative performance of the "flex" and "fix" inversions differs between the two campaigns. Whereas the "flex" inversions are in close agreement with measurements from the BARCA-A campaign, they tend to underestimate the observed mixing ratios during BARCA-B. None of the inversions are capable of capturing the seasonal dynamics of the BARCA measurements. The seasonal bias correction hardly influences the difference between BARCA 1 and 2, because the ECMWF-derived regional humidities for the end of the dry and wet season are rather similar.

Inversion-estimated fluxes
In this section we analyse the inversion-derived surface fluxes of CH 4 , and their sensitivity to the SCIAMACHY bias correction. We concentrate on fluxes that are integrated annually and over large enough regions to avoid interpretation of results that are compromised by the generally limited spatiotemporal resolution of the global inversions. Figure 10 shows flux time series for the period of July 2003 to July 2010 for six large regions. Monthly time series of annual fluxes are derived using a boxcar filter with a width of a year. The uncertainties represent 2σ intervals of the annual uncertainty for each region, computed by integrating the posterior covariance matrix in space and time, in 6-month intervals. The first and final half years of the inversion time window have been removed to avoid influences of inversion spin-up and spindown, and end effects of the box car filter. Note that each line is a composite of three 4-year inversions since the full period is split up in three overlapping blocks (as described in Sect. 2). The 2-year overlap allows us to reduce end effects of each of these blocks. Full time series are constructed by connecting the blocks in the middle of the overlapping periods. The NOAA-only and "fix" inversions show insignificant jumps between the blocks. For the "flex" inversions, differences are larger because the optimized bias correction parameters in subsequent blocks may differ. Nevertheless, these differences remain small compared with the inversionderived inter-annual variability. Compared with the a priori fluxes, the a posteriori fluxes show increased inter-annual variability (IAV), by 40 % for NOAA-only and 70 % for the SCIAMACHY inversions for the regions and data shown in Fig. 10. This increase may in part be real, since the a priori fluxes do not account for some sources of inter-annual variation, such as the year-toyear variation in wetland area. In addition, we do not account for inter-annual variability in the OH sink. As the sink is not optimized in our inversion, the influence of its IAV on the observed mixing ratios is projected on the surface fluxes. Studies of OH using methyl chloroform point to inter-annual variations within ±5 % (see e.g. Holmes et al., 2013), which corresponds to variations in methane oxidation of about 25 TgCH 4 yr −1 . This is significant, but also difficult to account for in global inversions since the available a priori estimates vary substantially and the CH 4 measurements used in this study provide limited independent constraints on surface emissions and atmospheric sinks.        The increased IAV of the SCIAMACHY inversions can be explained by stronger signals of source variability in the SCIAMACHY data than in the surface measurements in regions such as the tropics, where there is a paucity of surface measurements. As discussed earlier, however, we have limited means to verify the performance of SCIAMACHY on inter-annual time scales. Therefore, we cannot exclude the possibility that the enhanced variability is caused in part by changes in the performance of the instrument over time. Some regions, such as "Globe" and "Tropics", show a sim-ilar phasing of the SCIAMACHY and NOAA-only inferred IAVs, but with a different amplitude. Those results are more likely explained by the strengthening of a common signal by joined observational constraints, rather than temporal variations in instrument performance.
The   (Marengo et al., 2008). Our inversion-derived flux anomaly is not easily explained by this climatological anomaly, since reduced emissions would have been expected during the first half of 2005 rather than the second half extending into 2006. An alternative explanation could be the loss of detector pixels in the 1.6 µm absorption band of CH 4 in SCIAMACHY's channel 6. Artificial jumps in the retrieval due to variations in the dead and bad pixel mask are avoided in the retrieval by using only those pixels that were functioning properly throughout the whole analysed period. The instrumental damage nevertheless had an impact on the retrieval, as indicated by a marked increase in the variability of the data after November 2005 reported by Frankenberg et al. (2011).
The difference between the red and green lines in Fig. 10 highlights the impact of SCIAMACHY bias correction. Overall, it leads to a reduction of tropical emissions compensated by increases in the extra-tropics, consistent with the latitudinally varying adjustments of the total columns discussed earlier. The size of this difference varies between about 10 and 50 Tg yr −1 between the regions that were analysed. For emissions integrated over large zonal bands, as shown for "Tropics" and "ExtraTrop NH", SQG flex is fairly consistent with NOAA-only, although towards smaller regions the differences become larger. Despite the shift towards higher latitudes in the bias-corrected inversions, tropical emissions remain higher compared with the a priori by 40 Tg yr −1 , and 10 Tg yr −1 compared with the NOAA-only inversion.
Average emissions over tropical America vary between 75 and 88 Tg yr −1 for our SCIAMACHY inversions, in close agreement with the 79 Tg yr −1 reported in Frankenberg et al. (2008a) using SCIAMACHY IMAPv50 retrievals. Our NOAA-only emissions, however, are lower by 27 Tg yr −1 compared with that study, which can be explained by the SF 6 -tuned inter-hemispheric exchange in TM5 (Monteil et al., 2013). The most significant increase in tropical emissions is found for Africa, where the SQx inversions show on average a 60 % increase compared with the prior. In contrast to this, the NOAA-only inversion remains close to the prior.
The differences in temporal variability between S fix and SQG flex are relatively small compared with the latitudinal differences described above. This is true also for the other SCIAMACHY inversions, as can be seen in Fig. 11, and more clearly after removing the mean (see Fig. 12). The remaining differences are on average well within the posterior flux uncertainty, whereas the differences between the absolute fluxes grossly exceed those uncertainty ranges. It can be concluded that the biases of SCIAMACHY introduce important uncertainty in the spatial distribution of the global CH 4 emissions, but affect the deseasonalized temporal variability to a lesser extent. Note that the mean-subtracted time series show a striking difference in trend between the "flex" and "fix" inversions for the tropics. This arises from the independent optimization of bias parameters in the three inversion blocks. It leads to a systematic reduction in the long-term trend, with important implications for the inversion-derived trend over Southeast Asia.

Discussion
Overall, the results presented in the previous section confirm that the SCIAMACHY-inferred CH 4 emissions are sensitive to the bias correction method. This is true in particular for the emission difference between the tropics and the extratropics. The range of results obtained using different bias correction formulations exceeds the inversion-derived posterior flux uncertainty by up to a factor of five. Differences in inter-annual variations are much smaller, but difficult to evaluate because of the low density of the TCCON network during a large part of the SCIAMACHY operation. Further analysis of the fit residuals at background sites of the NOAA network (not shown) showed no clear relation with the differences between the NOAA and SCIAMACHY derived fluxes, which would have hinted to remaining inter-annually varying biases. Because of this, no further attempt was made to investigate the use of inter-annually varying bias functions. The time dependence of the inversion-optimized bias parameters points to a longer-term drift, which influences the inferred emission trends most clearly in the tropics. Regression of the residuals between the long-record TCCON time series and inversion-optimized TM5 confirms a drift using "fix" bias corrections of 1.6 ppb yr −1 . The largest drift is found for Darwin, although the regression is influenced by relatively large residuals at the end of the record. Without Darwin the slope of the regression drops to 0.6 ppb yr −1 , and is adjusted to −0.2 ppb yr −1 when bias parameters are optimized in the inversion.
The comparisons with independent measurements show residuals with changing sign dependent on the bias correction method. It indicates that the systematic errors we are trying to correct are an important factor limiting the performance of the inversion-optimized model. Extended bias Atmos. Chem. Phys., 14, 3991-4012, 2014 www.atmos-chem-phys.net/14/3991/2014/  correction can further reduce the difference between SCIA-MACHY and NOAA-only inversions as demonstrated by Bergamaschi et al. (2013). However, the options for deriving such corrections from independent data such as TCCON are limited. This highlights the importance of measurements from monitoring networks for surface and total column CH 4 , as well as aircraft measurement programmes. The recent expansion of the TCCON network greatly improved its usefulness, from which current and future satellite missions will benefit much more than SCIAMACHY. Still, the coverage is not yet sufficient, as indicated by the results of this study, which would have benefited from extended coverage over tropical continents.
To investigate common signals of SCIAMACHY in our inversions and those reported by Bergamaschi et al. (2013), the time series of inversion-derived fluxes are compared in Fig. 13. Although Bergamaschi et al. (2013) (referred to as B13 hereafter) make use of TM5 4DVAR also, the inversion setups are quite different. Besides the use of different bias correction functions, important differences include the use of either Gaussian or log normal probability density functions for the surface fluxes, optimization of net fluxes or separate flux components, and the selection of measurement sites. The differences between the NOAA only and SCIA-MACHY+NOAA inversions are smaller in B13, which can be explained by the larger number of bias correction parameters that are accounted for (288 vs. 9 in our "flex" inversions). Some common differences between inversions with and without satellite data are found, such as the minimum in emissions over tropical America during 2006.
Besides the differences in inter-annual variability highlighted in Fig. 13, multi-year mean fluxes show important differences also. A summary of mean fluxes is presented in Fig. 14, including all inversions of this study combined  with results from B13. Integrated over the globe, the average of the posterior emissions of our inversions amount to 527 ± 4 TgCH 4 yr −1 , in close agreement with B13. Note that, despite the significant differences in inter-annual variation between the inversions, the 8-year integrated fluxes span only a small range. The Northern Hemisphere accounts for 78 % of global emissions, without a significant shift from a priori to a posteriori. The results of B13 show a shift towards the Southern Hemisphere of 42 Tg yr −1 compared with our inversions, consistent with the SF 6 calibrated interhemispheric exchange used in this study but not in B13. A sizable shift of 50 Tg yr −1 is found from the extra-tropics to the tropics, although this number varies substantially between the inversions (σ = 25 Tg yr −1 ). S fix and S flex show a net sink in the SH extra-tropics, which is brought in a more realistic range when the seasonal bias correction is used. Since we assume Gaussian flux uncertainties, nothing prevents the inversion from changing the sign of the fluxes. Note that the values found here are small compared with the corresponding uncertainties.
To address the question of which region contributed most to the CH 4 increase since 2007, the difference is plotted between a 2-year period before and after July 2006, when according to our inversion results the emission increase started (see Fig. 15). As can be seen, the global difference varies between 27 and 35 Tg yr −1 across our inversions, most of which is attributed to the tropics with the northern hemispheric part of this zone contributing most. Splitting the tropics up by continent the resolution of the inversion becomes limited, as indicated by the size of the error bars exceeding the flux difference. The largest portion is attributed to Southeast Asia (9 ± 13 Tg yr −1 ), consistent with the growing demand for energy and food in rapidly developing economies. However, the low significance of this value does not allow us to draw strong conclusions.
Atmos. Chem. Phys., 14, 3991-4012, 2014 www.atmos-chem-phys.net/14/3991/2014/ Fig. 13. As Figure 12, with comparison to the results from the NOAA-only (S1-NOAA) and NOAA+SCIA (S1-SCIA) inversion from Bergamaschi et al. (2013). Compared to our estimates, B13 report a smaller emission increase across the 2007 growth rate transition (16-20 vs. 27-35 Tg yr −1 ). The B13 estimates in Fig. 15 are slightly higher than reported in the original publication, due to a difference in the time intervals used to quantify the emission increase. The same holds for the comparison with numbers reported by Bousquet et al. (2011). Their INV2 shows a difference between 2006 and 2007-2008 of 23 Tg yr −1 , compared to 19 Tg yr −1 as the average of our inversions for the same period. Our inversions show smaller increases using these time intervals, because emissions start increasing already in the course of 2006. Nevertheless, evaluated over a 4-year time period our SCIAMACHY inversions show larger increases compared with B13 and Bousquet et al. (2011), which we attribute to the reduction in the IMAPv5.5 retrieved XCH 4 in the tropics in 2006. Bousquet et al. (2011) conclude that an important fraction of the increase can be attributed to increasing natural wetland emissions. Although our inversion setup does not separate between different processes, the emission increase over Southeast Asia agrees well with B13 who point to increasing anthropogenic emissions in that region. Our uncertainty analysis, however, indicates that the inversion is incapable of resolving the difference between these explanations. The three studies do agree on an important contribution from the tropics. The answer to the question of which processes are primarily responsible may come from isotopic analysis. The recent downward trend (see http://www.esrl.noaa.gov/gmd/  Fig. 14. Comparison of inversion derived fluxes integrated over regions and over the period 2003-2010, including results from the NOAAonly (S1-NOAA) and NOAA+SCIA (S1-SCIA) inversion from Bergamaschi et al. (2013). Error bars represent 2σ uncertainties. For the region definition see Figure 1.      dv/iadv/) in global δ 13 C-CH 4 can be interpreted as a shift in the global CH 4 source mixture towards microbial methane production. However, even in the absence of an isotopic shift in the source mixture, a change in the CH 4 growth rate itself is expected to deplete atmospheric δ 13 C-CH 4 by increasing the isotopic disequilibrium. Whether this can account for the observed trend will have to be investigated in further detail in the future.

Conclusions
We have investigated the use of SCIAMACHY retrievals spanning multiple years, for estimating temporal variations in the global sources of CH 4 . Inverse modelling calculations were performed using TM5 4DVAR, for the period 2003-2010. The challenge of using CH 4 retrievals from SCIAMACHY is to account for systematic retrieval uncer-tainty, which complicates the interpretation of inverse modelling results. We made use of the available TCCON measurements to define and calibrate bias correction functions. In comparison with TCCON, the SCIAMACHY IMAPv5.5 XCH 4 retrievals show an important seasonally varying bias, which can be effectively corrected using specific humidity in the lower troposphere as a proxy. Inverse modelling calculations were carried out using alternative formulations of the bias correction functions, with the coefficients either calibrated to TCCON or optimized in the inversion. Evaluation of the results against aircraft measurements from the HIPPO and BARCA campaigns shows important sensitivity to the bias correction approach, confirming that systematic uncertainties are a critical factor limiting the performance of the inversion. The TCCON total column measurements are of crucial importance for identifying and quantifying those uncertainties. To take full advantage of satellite Atmos. Chem. Phys., 14, 3991-4012, 2014 www.atmos-chem-phys.net/14/3991/2014/ retrievals, extended validation capabilities are needed, in particular, over tropical continents. Our inversions shift emissions from the extra-tropics to the tropics by 50 Tg yr −1 compared with the prior, but this number varies by 25 Tg yr −1 among the different setups that were tested. The inversion-derived inter-annual variability is less sensitive to the bias correction method, in part because limited information is available to base such corrections on. Comparisons with inversion results from Bergamaschi et al. (2013) show that the SCIAMACHY-derived interannual emission variations become less robust when the bias correction is extended with additional degrees of freedom on inter-annual time scales. Integrated over large scales, the use of SCIAMACHY data leads to a 22 % increase in the inversion-estimated inter-annual variability compared with the use of only surface measurements. Increased variability is expected since the use of SCIAMACHY improves the coverage over tropical continents, which host important sources of IAV such as biomass burning and tropical wetlands. The fact that SCIAMACHY amplifies tropical variability that is already seen in the NOAA-only inversion, increases our confidence that part of the difference is real. In our inversions, the observed transition to increasing CH 4 mixing ratios in 2007 is attributed mostly to the tropics. The difference in global emissions between a 2-year period before and after July 2006 amounts to 27-35 Tg yr −1 . Within the tropical band an important contribution is found from Southeast Asia, although the associated posterior flux uncertainties are too large to identify the emissions from growing Asian economies as the main cause. Therefore our results are also consistent with a scenario of coincident increases in emissions from tropical wetlands.