the Creative Commons Attribution 3.0 License. Atmospheric Chemistry

A four-dimensional variational (4D-Var) data assimilation system for inverse modelling of atmospheric methane emissions is presented. The system is based on the TM5 atmospheric transport model. It can be used for assimilating large volumes of measurements, in particular satellite observations and quasi-continuous in-situ observations, and at the same time it enables the optimization of a large number of model parameters, specifically grid-scale emission rates. Furthermore, the variational method allows to estimate uncertainties in posterior emissions. Here, the system is applied to optimize monthly methane emissions over a 1-year time window on the basis of surface observations from the NOAA-ESRL network. The results are rigorously compared with an analogous inversion by Bergamaschi et al. (2007), which was based on the traditional synthesis approach. The posterior emissions as well as their uncertainties obtained in both inversions show a high degree of consistency. At the same time we illustrate the advantage of 4D-Var in reducing aggregation errors by optimizing emissions at the grid scale of the transport model. The full potential of the assimilation system is exploited in Meirink et al. (2008), who use satellite observations of column-averaged methane mixing ratios to optimize emissions at high spatial resolution, taking advantage of the zooming capability of the TM5 model.


Introduction
Inverse modelling has been widely used as a tool to improve our knowledge on sources and sinks of atmospheric trace Correspondence to: J. F. Meirink (meirink@knmi.nl)gases based on measurements of concentrations in the atmosphere (Enting, 2002).Since such inversions of atmospheric transport are often ill-conditioned, a priori information on the spatial and temporal distribution of sources and sinks derived from emission inventories is normally used to regularize the problem.The goal of inverse modelling is then to optimize these inventories and reduce their uncertainties using observations.
Regarding atmospheric methane, mainly surface observations have been used for this purpose so far (Hein et al., 1997;Houweling et al., 1999;Mikaloff Fletcher et al., 2004;Bergamaschi et al., 2005;Chen and Prinn, 2006;Bousquet et al., 2006).While the existing network of surface measurements constrains global emissions relatively well, large continental regions (e.g. in the tropics) are poorly monitored.Quasi-continuous in-situ observations close to source regions can provide important information on regional emissions (Bergamaschi et al., 2005), but the number of such measurements remains very limited.Recently, satellite observations of column-averaged methane mixing ratio became available from the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) instrument on board ESA's environmental satellite ENVISAT (Buchwitz et al., 2005;Frankenberg et al., 2005Frankenberg et al., , 2006;;Buchwitz et al., 2006).Bergamaschi et al. (2007) (hereafter B07) used these observations -along with the conventional surface measurements -for the first time to optimize continental-scale emission rates.
Traditionally, most inverse modelling studies have been based on the synthesis approach (Enting, 2002).This approach is mainly applicable when emissions are optimized on a coarse resolution (e.g., for a limited number of pre-defined regions).Grid-based optimization has been performed using Published by Copernicus Publications on behalf of the European Geosciences Union.
the so-called adjoint technique (Houweling et al., 1999;Kaminski et al., 1999), but was in these studies restricted to relatively small sets of observational data.To take full benefit from the existing and future satellite measurements, approaches that can handle large amounts of observations together with a large control vector are required.Two branches of promising and computationally feasible techniques (of which the latter also uses the adjoint model) are ensemble data assimilation (Evensen, 1994) and fourdimensional variational (4D-Var) data assimilation (Talagrand and Courtier, 1987), both building upon developments in Numerical Weather Prediction (NWP).
In recent years, variational data assimilation has been increasingly applied for the estimation of emission rates of atmospheric constituents.Examples include the inversion of emission rates based on TIROS Operational Vertical Sounder (TOVS) CO 2 observations (Chevallier et al., 2005), Measurements of Pollution in the Troposphere (MOPITT) CO measurements (Stavrakou and Müller, 2006), surface observations of various short-lived trace gases (Elbern et al., 2007), and LIDAR measurements of dust aerosols (Yumimoto et al., 2007).In synthetic frameworks, 4D-Var has been used to study the utility of CO 2 observations from the planned Orbiting Carbon Observatory (OCO) mission (Baker et al., 2006;Chevallier et al., 2007) and CH 4 observations from SCIA-MACHY (Meirink et al., 2006) for reducing uncertainties in emission estimates.An ensemble data assimilation system for inverse modelling of surface fluxes was, for example, presented by Peters et al. (2005).
The purpose of this paper is to present and evaluate a new 4D-Var system for inverse modelling of methane emissions.This new system is a further development of the work by Meirink et al. (2006).The main changes include the use of TM5 instead of TM4 as underlying atmospheric transport model, the application of a different minimization algorithm, and an improved implementation of the background (also termed prior) error covariance matrix and preconditioning.A major advantage of the new set-up is that it allows to estimate uncertainties of the posterior emissions.
To evaluate the new system, a 1-year inversion of surface observations will be compared in detail with an analogous synthesis inversion by B07.The comparison shows a high degree of consistency, but also illustrates the advantage of the variational method.In Meirink et al. (2008) the 4D-Var system is used to analyse SCIAMACHY observations with focus over South America, exploiting the zooming capability of the TM5 model.
The paper is structured as follows.In Sect.2, the transport model and its adjoint as well as the 4D-Var implementation are presented.Section 3 contains the results of the comparison with B07's synthesis inversion.Some specific issues related to convergence of the variational algorithm, diagnostics of the assimilation system and aggregation errors are discussed in Sect. 4. Finally, the main conclusions are summarized in Sect. 5.

Transport model and adjoint
The TM5 model is a global chemistry-transport model with a two-way nested zooming capability (Krol et al., 2005).TM5 is driven by meteorological fields (6-h forecasts) from the European Centre for Medium Range Weather Forecasts (ECMWF) operational model.We use the methane tracer version as described in B07, in which chemical oxidation of CH 4 is calculated from prescribed OH fields.These OH fields were obtained from a full-chemistry model simulation and calibrated with methylchloroform observations, resulting in a tropospheric CH 4 lifetime of 9.4 years (Bergamaschi et al., 2005).The model is operated on a basic horizontal resolution of 6 • ×4 • globally.The zooming option via 3 • ×2 • to 1 • ×1 • nested grids in specific regions is not used here, but in Meirink et al. (2008).In the vertical direction, 25 layers have been defined as a subset of the 60 layers used operationally in the ECMWF model in 2003.Emissions are described in Sect.2.3.
For 4D-Var assimilation, an adjoint model is required.The coding of the adjoint was done manually by matrix transposition, except for the slopes advection scheme (Russell and Lerner, 1981), for which an adjoint was originally generated by the Tangent and Adjoint Model Compiler (TAMC) (Giering and Kaminski, 1998).For the different processes (vertical diffusion, convection, oxidation by OH, and emissions) the construction of the adjoint was a rather straightforward task.Most complications encountered were related to the merging and division of grid cells, which occurs in the communication between parent and child regions in the zoom grid and in the reduced grid that is applied near the poles to ensure numerical stability at reasonably large time steps.The exact reconstruction of so-called inactive variables, such as the air masses, was also a difficult task.More details on the construction of the adjoint model can be found in the supplementary material (http://www.atmos-chem-phys.net/8/6341/2008/acp-8-6341-2008-supplement.pdf).The adjoint was used recently in a study on the usefulness of 14 CO measurements for determining OH concentrations (Krol et al., 2008).
The correctness of the adjoint was verified by checking the equality where M and M T denote the forward and adjoint TM5 model operators applied over a certain time window, x and y are arbitrary forward and adjoint model states, and denotes the inner product.In the standard TM5 model the tracer slopes in the advection scheme are limited to avoid negative concentrations at the edges of grid boxes.These limiters represent a non-linearity, but in the methane tracer version they are not needed, so that the forward model operator is linear and can thus be written as a matrix M. It was verified that the relative difference between the left-hand and right-hand terms in Eq. ( 1) is ∼ 10 −14 (i.e. machine precision), for arbitrary time windows up to ∼ 1 year and arbitrary model states x and y.

4D-Var
Compared to Meirink et al. (2006), the set-up of 4D-Var has been further updated in a number of aspects, most importantly the minimization algorithm.Therefore, the major components of the previous system are summarized and the new developments are described in detail.It should be noted that 4D-Var as applied in the context of source and sink optimization differs from the conventional use in NWP.In the latter, measurements are used to optimize three-dimensional meteorological fields over time windows on the order of a day, serving as initial conditions for forecasting.Here, measurements are used to optimize the two-dimensional distribution of surface emissions (although the three-dimensional initial concentration field is also optimized), and the time window is on the order of months.However, the mathematical framework is identical.
The optimization problem has the following ingredients: 1.A set of observations y with a corresponding error covariance matrix R. The vector y contains the available observations during the time window of assimilation.
For convenience the time dimension is not explicitly denoted here.
2. A set of model parameters x (the control vector) with a corresponding background error covariance matrix B. In our case, the control vector can be written as x = (s T 1 , . . ., s T m , c T , p T ) T , where s i are monthlymean surface emissions for source category i and m is the number of source categories that are distinguished, c is the three-dimensional concentration field at the start of the assimilation window, and p contains any additional parameters, for example model parameters in parameterizations of physical processes.In the present paper, p is not used, but in Meirink et al. (2008) it contains parameters that model a bias in the satellite observations that are assimilated.

A model operator H, with
Hx providing the equivalent of the observations y.This operator consists of application of the forward model M to the time of the measurements followed by application of an observation operator, which interpolates and/or vertically integrates methane concentrations on the model grid to produce equivalents of the observations.Since this observation operator is linear in our case (also for the SCIAMACHY observations assimilated in Meirink et al. (2008)), the operator H is linear.
The goal is to find the optimal control vector given the background estimate and observations taking into account their respective uncertainties.In 4D-Var this analysis x a (also termed posterior solution) is obtained by iteratively minimizing the following cost function J with respect to x: The minimization algorithm requires calculations of the gradient of the cost function: After further differentiation, it can be seen that the Hessian of the cost function is independent of x: Starting from the property that Eq. ( 3) evaluated at x a is zero, it can be shown (e.g., Fisher and Courtier, 1995) that the covariance matrix of analysis errors (hereafter mostly referred to as posterior errors) A equals the inverse of the Hessian: This equation illustrates that the assimilation of observations leads to posterior errors that are smaller than the prior errors.
In other words: the inversion yields an error reduction (also termed uncertainty reduction in this paper), which is defined as 1 − σ a /σ b , where σ b and σ a denote the (potentially aggregated) prior and posterior errors, respectively.We use the same minimization method as employed in the ECMWF 4D-Var (Fisher and Courtier, 1995).A more detailed description can be found there and also in Chevallier et al. (2005).In short, the method is based on the Lanczos algorithm (Lanczos, 1950) and allows to simultaneously minimize the cost function and derive the leading eigenvalues and eigenvectors of the posterior error covariance matrix.With preconditioning (Courtier et al., 1994), which is applied to reach faster convergence to the minimum of the cost function, the Hessian with respect to the preconditioned control variable χ becomes After N iterations, the minimization algorithm has estimated the N leading eigenvalues λ i and eigenvectors ν i of this matrix (with λ 1 >λ 2 > . . .>λ N ).Writing the Hessian in terms of this eigen decomposition, inverting, and transforming back from χ to x, gives the approximation of the posterior error covariance matrix: This approximation converges to the true posterior error covariance matrix as λ i goes to one with increasing number of  iterations.From Eq. ( 8) it follows that the diagonal elements of the approximation of A strictly decrease with every iteration added.Therefore, the method always yields an overestimate of the diagonal elements of A, which represent the posterior variances.The a priori error covariance matrix B is much too large (typically ∼ 10 11 numbers) to be stored in memory.However, by making some reasonable assumptions, the required storage can be sharply reduced.First, the covariance matrix is split up in a correlation matrix C and a diagonal matrix D containing the standard deviations: Furthermore, the following assumptions are made: (i) the emissions in different source categories, the initial concentration field, and the additional parameters have mutually uncorrelated errors; (ii) the error correlation matrices for emissions, C s i , can be split into independent spatial, C h s i , and temporal, C t s i , parts; (iii) the error correlation matrix for initial concentrations, C c , can be split into independent horizontal, C h c , and vertical, C v c , parts.With these assumptions, we can write: where ⊗ is the Kronecker matrix product and C p is the error correlation matrix for the additional parameters.Herewith, the storage requirements are reduced to a number of C h matrices, which contain the square of the number of grid points in a horizontal plane.This is typically on the order of 10 7 , small enough to be stored in computer memory.The C h matrices are modelled by a Gaussian function of the distance between grid cells, with pre-defined correlation length scales L s i and L c .Similarly, the C t s i are modelled by an exponential function of the time difference, with pre-defined correlation time scales τ s i .As outlined in Meirink et al. (2006), C v c is determined with the National Meteorological Center (NMC) method (Parrish and Derber, 1992).Finally, C p is usually an identity matrix, stating that the errors in the parameters are uncorrelated.The square roots of the various correlation submatrices, needed for the preconditioning (Eq.6), are determined by eigenvalue decomposition.

Inversion set-up
The CH 4 inversions presented in B07 provide an excellent test case for the newly developed 4D-Var assimilation system.Briefly, scenario S1 in B07 comprised a 1-year synthesis inversion of monthly-mean methane emissions from 11 source categories and for a number of large geographical regions, based on surface observations from the NOAA ESRL network.We have repeated this inversion with our 4D-Var approach, copying the set-up as well as possible.Specifically, we use the same forward model, including meteorological input, the same a priori emissions, and the same set of surface observations and corresponding errors.
Total a priori emissions for the year 2003 are listed in Table 1 for the 11 source categories and for -depending on the category -1 to 7 large regions as shown in Fig. 1.
Methane surface observations are taken from the NOAA ESRL global cooperative air sampling network (Dlugokencky et al., 1994(Dlugokencky et al., , 2003)).Only flask measurements from marine and continental background sites are used.The 32 selected sites are listed in B07 (Table 1) and depicted in Fig. 2. The y i in the cost function are individual surface  near large CH 4 gradients.Measurement errors are assumed to be uncorrelated.The inversion is carried out in two cycles.
In the second cycle only those observations are assimilated that differ less than three times the observation error from the posterior model simulation of the first cycle.This approach was used by B07 to get rid of the possible detrimental impact of observations that cannot be reproduced by the model simulation even after modifying the surface emissions.
The only important difference between the present inversion and B07 is the control vector and therefore also the prior error covariance matrix.In B07, the control vector consisted of monthly emissions for the different source categories and regions plus a scaling factor for the initial concentration field.In the 4D-Var, the control vector consists of monthly emissions for the same categories but for all individual grid cells plus the grid-scale initial concentration field.For technical reasons, we did not perform a two-months spin-up as in B07, but we did include January 2004 in the assimilation window, since the measurements in this month provide a constraint on emissions in the last months of 2003.The dimension of the control vector thus becomes: (13 months×11 source categories+25 vertical levels of the initial concentration field)×(60×45 grid cells)=453 600, compared to 541 in B07.The actual number of degrees of freedom in the 4D-Var control vector is only a small fraction of its dimension due to the assumed spatial and temporal error correlations.
The a priori errors of grid-cell emissions are set to the same relative levels as the a priori errors of big-region emissions in B07, varying from 30% for enteric fermentation to 80% for, e.g., wetlands.Spatial correlations are modelled by normal distributions where the correlation length scales have been chosen such that the a priori errors aggregated to big regions are close to those in B07.This yields very large length scales of 20,000 km for wild animals, termites, soil oxidation and ocean, 4000 km for wetlands, and 5000 km for the other categories.Temporal correlations are specified exactly the same as in B07, i.e.: no correlations for wetlands, rice paddies and biomass burning, and exponential correlations with a 9.5-months correlation time scale (month-to-month correlation of 0.9) for all other sources.Table 1 shows the resulting a priori errors aggregated to a whole year and to big regions.They are indeed close to the values in the synthesis inversion.
Apart from the above reference inversion (termed scenario A), an alternative inversion was conducted (scenario B), in which the prior error spatial correlation lengths were set to 1000 km for all source categories.To arrive at the same globally aggregated prior uncertainty as in the reference inversion, it was found that the prior standard deviations had to be multiplied by a factor ≈2.5.

Convergence
The convergence of the minimization algorithm is analysed in Fig. 3.The norm of the gradient of the cost function decreases rapidly and steadily.The iteration is stopped when a reduction of 10 10 is reached.This point is reached after 56 iterations for the present inversion.A good proxy for convergence are the eigenvalues of the Hessian, estimated by the Lanczos algorithm.Fig. 3b shows that the eigenvalues decrease rapidly and have reached values close to one in the last iteration.A more direct way of evaluating convergence, is to inspect the estimate of the Hessian itself, and changes therein.Figs.3c and d illustrate this in the form of uncertainty reduction, defined in Sect.2.2.At the end of the inversion, the mean grid-scale uncertainty reduction reaches a value of around 5%.As will be discussed in Sect.4.1, full convergence has not yet been reached at this stage.When aggregated over months, categories and large regions, the uncertainty reduction converges much faster (Fig. 3d), and stable estimates are obtained after about 45 iterations.Note that the uncertainty reduction on the grid scale is very low (and it is even lower when the spatial error correlation length is decreased).However, it becomes much larger when aggregated over months, categories and large regions, showing that the measurements are indeed useful to constrain emissions on larger spatial and temporal scales.
In contrast to the errors in posterior emissions, the posterior emissions themselves have converged in about 20 iterations (not shown).Thus, if one is interested in emissions only, relatively few iterations are sufficient, but if reliable posterior error estimates are needed, many more iterations are required.This implies that large gradient norm reductions must be achieved, for which an exact adjoint model is a prerequisite.

Emission increments
In Fig. 2  Both inversions propose enhanced emissions in the tropics and decreased emissions in the northern-hemispheric (NH) high latitudes.This large-scale increment is mainly driven by the need for the system to correct the too large interhemispheric (IH) concentration gradient in the a priori TM5 simulation and reduce the emissions of NH wetlands during summer.When a closer look at the details is taken, some differences emerge.Most notably, the emission increments in South-East Asia are different.The 4D-Var suggests a quite strong decrease in emissions, while the synthesis inversion  provides a more mixed pattern.Further analysis shows that these differences can be largely attributed to opposite increments in emissions from rice cultivation (see below).In principle, 4D-Var is more flexible by optimizing (spatially correlated) grid-cell fluxes versus fixed flux patterns in a few big regions in the synthesis inversion.In order to compare to the synthesis inversion, very large correlation length scales have been used in the present 4D-Var inversion, so that this difference in flexibility should be relatively small.It is thus not completely clear what causes the specific differences between both inversions on the smaller scales.Certainly, the problem at hand has many near-optimal solutions.Therefore, it is probably to be expected that both inversion systems come up with slightly different answers.
Table 1 gives a more quantitative comparison of the emissions estimated by the synthesis and 4D-Var inversions.The overall agreement is good, although for some categories there are substantial differences.For example, the global total rice emission increment is +3.7 Tg and −4.0 Tg for the synthesis and 4D-Var, respectively.The difference is, however, still only slightly larger than the posterior uncertainty estimates of 6.7 and 6.5 Tg, respectively.Interestingly, both inversions yield nearly vanishing posterior oceanic emissions.In Sect.4.3 we will show that this is likely an artefact related to the large regions / large correlation lengths employed.

Comparison with measurements
Figure 4 shows a comparison of prior and posterior model simulations with observations at the same 8 stations displayed in B07's Fig. 4. In general, the observed seasonal cycles and synoptic-scale events are captured very well already in the prior model simulation.There are two clear exceptions.First, at the NH high-latitude stations, here represented by Barrow, the prior simulation is much too high, particularly in summer.This is explained by the high boreal wetland emissions in the applied a priori inventory, which are substantially reduced by the inversion, leading to a good correspondence between the posterior model simulation and the observations.Second, at the southern-hemispheric (SH) stations, the prior simulation is drifting away from the observations, leading to an increasing underestimation in the course of the year.This feature reflects the general overestimation of the IH gradient in the free-running model.The cause of this overestimation may be an error in the prior emissions or in the seasonality and IH gradient of OH, but could also be related to errors in model transport.Detailed comparison of SF 6 simulations with observations also showed a slight overestimation of the IH SF 6 gradient, although the derived IH mixing time of 10.4 months was well within the typical range of atmospheric transport models (Bergamaschi et al., 2006).Since model transport and OH distribution are assumed to be perfect in the inversion, the overestimated IH gradient is corrected by enhancing emissions in the tropics and the SH.
The green curves in Fig. 4 correspond with a simulation based on optimized emissions from B07's synthesis inversion.The time series from both inversions are hardly discernible at most locations most of the time, highlighting again the high degree of consistency between the 4D-Var and the synthesis inversions.

Uncertainty reduction
The posterior errors derived from the 4D-Var inversion using Eq. ( 8) have been aggregated to the big regions of the synthesis inversion in Table 1.The resulting uncertainty reductions in the 4D-Var and synthesis inversions are very similar.For the global total emission both approaches yield an error reduction estimate of almost 90%.Differences of more than 10% occur for only four category-region combinations.
From Eq. ( 5) it is evident that the uncertainty reduction of a surface flux is determined by (i) the extent to which this flux is seen by the observations and (ii) the ratio of the a priori error (which is assumed to be proportional to the a priori flux) and the observation errors.Indeed, the uncertainty reductions are generally largest for the NH regions in which most measurement stations are located and over which substantial emissions take place.At first glance, the ocean category seems to be a notable exception: it shows a large error reduction of around 50% despite its modest prior uncertainty.The explanation for this large error reduction is given in Sect.4.3.

Convergence
There are a few issues evolving from Fig. 3 that deserve further discussion.First of all, the mean grid-scale uncertainty reduction (panel c) may appear to have converged after 55 iterations, but when more iterations are performed, it still shows an additional increase (not shown).From a series of inversions we could accurately approximate the convergence behaviour by an exponential function.From this approximation it is estimated that for the present inversion the actual (converged) mean uncertainty reduction is about 7.5% instead of 5%.Correspondingly, the Hessian eigenvalues are close to 1 after 55 iterations, but remain large enough in further iterations to yield non-negligible contributions to the posterior error in Eq. 8.
Experiences with similar inversions of CO 2 emissions show a slower convergence of Hessian eigenvalues than in our case (F.Chevallier, personal communication).Obviously, the speed of convergence is related to prior and observational errors and their correlations.One reason why there is slower convergence for CO 2 may be that the prior error distribution is more spatially homogeneous than used here for CH 4 , giving effectively more degrees of freedom to modify the emissions.Indeed, a test inversion with homogeneous CH 4 a priori error distributions took about 40 iterations more to converge than the reference inversion.
A second feature is that posterior errors aggregated over space and time converge much faster than on the grid scale (compare Fig. 3c and d).This can be explained as follows (see also Fisher and Courtier, 1995;Chevallier et al., 2005): in the first iterations the large-scale emission patterns are determined, while in later iterations the smaller-scale patterns are finetuned.The increments in later iterations are mainly located in regions where the prior error is relatively small.These increments do influence the mean grid-scale uncertainty reduction, but they have very little impact on the aggregated uncertainty reduction.Most often we are interested in uncertainty reductions over somewhat larger regions than the grid scale.For this purpose the Lanczos algorithm appears to give accurate estimates within a reasonable number of iterations.This remains true when large volumes of satellite data are assimilated, as is shown in Meirink et al. (2008).

Diagnostics
A useful diagnostic, indicating whether the assimilation may not be optimal, i.e. measurement and prior errors have been improperly set, is the χ 2 (e.g., Tarantola, 2005): When x b and y are interpreted as random vectors, being randomly distributed around the truth x t and y t , respectively, the random variable in Eq. ( 11) should follow a χ 2 distribution with n degrees of freedom, where n is the number of observations.In 4D-Var, H and H T are not available in matrix representation.Therefore, we use an equivalent representation in terms of x a (Tarantola, 2005): For the reference inversion, a satisfying value of χ 2 /n=1.30 is obtained.The alternative inversion scenario B with shorter spatial error correlation lengths of 1000 km yields χ 2 /n=1.08,thus somewhat closer to 1, reflecting the larger number of degrees of freedom in the control vector.Whereas the χ 2 -criterion from Eq. ( 11) can be naturally applied to a subset of the observations (e. or within a specific time window), this is not the case for Eq. ( 12) as it contains a background term which cannot be coupled to specific observations.Still, the analysis-minusobservation residuals Hx a − y are often studied for selected observations.In case of uncorrelated observation errors, one then writes: where σ y i is the uncertainty of measurement y i in an observation subset with size n s .Given the omission of the background term, the calculated χ 2 should be smaller than n s .However, in our case the background term of the cost function contributes typically 10% to the total cost, so that χ 2 /n s from Eq. ( 13) should still be close to 1. B07 obtained χ 2 /n s = 1.2 for the whole set of observations.Our reference inversion yields the same value, while scenario B yields χ 2 /n s = 0.9.It thus appears that 4D-Var, having a higher flexibility in changing emission patterns, is indeed able to reach a closer fit to the observations than the synthesis inversion, but when it is run with large correlation lengths the fit to the observations is very similar.It should be noted that χ 2 /n s is considerably influenced by the 2-cycle inversion approach (see Sect. 2.3).In the first cycle χ 2 /n s is about 1.8, much higher than in the second cycle, because posterior outlier observations (typically 2.5% of the total number of measurements) have not yet been removed.
Between the measurement stations there is considerable variability in χ 2 /n s .The lowest values (χ 2 /n s ≈0.2) are obtained for the Antarctic stations, indicating that the measurement precision estimate of 3 ppb may be too conservative.On the other hand, χ 2 /n s goes up to 2 at some stations (e.g.AZR).These high values may be related to unresolved smallscale variability in the model or to an underestimate of the representation error for some sites.The wide range of chisquared for different stations demonstrates that one should be careful when using a single global chi-squared statistic for diagnosing the overall performance of the assimilation system.

Uncertainty reduction and aggregation error
As shown in Table 1, the posterior ocean emissions in both the synthesis and the 4D-Var inversions are nearly zero, while at the same time a large uncertainty reduction of ≈50% is calculated.This feature is related to the so-called aggregation error (Kaminski et al., 2001;Engelen et al., 2002).Due to the large spatial error correlation lengths assumed in the 4D-Var inversion and the big regions applied in the synthesis inversion, in fact all marine background stations contribute to an uncertainty reduction of the global ocean emission.
The alternative inversion scenario B with shorter error correlation lengths (see rightmost columns of Table 1) further illustrates the impact of the aggregation error.The uncertainty in the global ocean emissions is now reduced by only 5%.At the same time, the ocean emissions are hardly reduced.Also for other source categories large changes in posterior emissions are obtained.In general, the posterior totals are much closer to the prior, except for wetland emissions, which increase by more than 20 Tg yr −1 .The global total emissions over all categories appear to be relatively tightly constrained at around 525 Tg yr −1 , given the OH distribution currently applied in our inversions.
The spatial patterns of the uncertainty reduction in inversion scenarios A and B are shown in Fig. 5.With increasing correlation length, the uncertainty reduction on the grid scale increases, and the spatial pattern shows broader structures.The reason for this is that more measurements contribute to uncertainty reduction in a certain grid box.In other words: the region of influence of a measurement increases artificially.The plot for the reference inversion looks a little noisy.In contrast, the spatial patterns of uncertainty reductions for individual categories are smooth (not shown).The noisiness is thus mainly caused by the aggregation over different source categories.
The possibility to estimate emissions at high resolution is a major advantage of 4D-Var compared to the traditional synthesis approach, since it reduces the risk of aggregation errors due to measurements having artifically large regions of influence.Another advantage of grid-based inversions is that the information present in observations with high spatial and temporal resolution, such as satellite and high-frequency insitu measurements, can be better exploited.It has to be noted that these advantages weaken considerably if the prior error correlations of emissions and observations cannot be properly determined, which is unfortunately often the case.In any case, the 4D-Var offers a flexible and computationally efficient means to perform inversions with both large control vectors and many assimilated observations.

Conclusions
We have presented a 4D-Var system for optimizing methane emissions based on observations of atmospheric methane concentrations.The main advantage of the 4D-Var system compared to the classical synthesis inversion is that emissions can be optimized at high spatial resolution and that large volumes of observational data can be taken into account.Furthermore, a specifically useful feature of our inversion system is that it allows to estimate uncertainties of the optimized emissions.A 4D-Var inversion with large spatial correlation lengths of prior emission errors was rigorously compared with a previously performed synthesis inversion, showing a high degree of consistency in terms of both posterior emission estimates as well as their uncertainties.At the same time, an inversion with smaller spatial error correlation lengths demonstrated the advantage of 4D-Var in reducing aggregation errors.The full benefit of the new system becomes apparent when larger numbers of observations are assimilated.This is illustrated in Meirink et al. (2008), who use SCIAMACHY satellite observations to infer global CH 4 emissions with a focus on South America.

Fig. 1 .
Fig. 1.Definition of geographical regions used in the synthesis inversion by B07.

Fig. 3 .
Fig. 3. Convergence of the 4D-Var minimization.As a function of iteration i are shown: (a) the log of the norm of the cost function gradient relative to the prior simulation, (b) the log of the eigenvalue λ i of the Hessian matrix, (c) the average uncertainty reduction of 1-month, 1-category, 1-grid-cell emissions, and (d) the aggregated emission uncertainty reduction over 1 year, all categories, and 5 different large regions (as represented in Fig. 1).

Fig. 4 .
Fig. 4. Comparison of model simulations based on prior emissions (blue), posterior 4D-Var emissions (red), and posterior synthesis inversion emissions (B07, scenario S1) (green) with observations (black symbols) at eight NOAA stations during 2003.Black vertical bars indicate the total observation error (including representativeness error).

Fig. 5 .
Fig. 5. Emission uncertainty reduction aggregated over the complete year 2003 and over all source categories for (a) 4D-Var reference scenario A and (b) scenario B using shorter prior error correlation lengths.Note the different color scales in both panels.

Table 1 .
Prior and posterior emissions and uncertainties for synthesis inversion of B07 (scenario S1) and two 4D-Var inversions aggregated over the complete year 2003 and over the regions defined in Fig.1.Units are Tg CH 4 yr −1 .The fractional uncertainty reduction is given in the columns 'U.r.'.Two 4D-Var inversions with different error correlation length scales (see text) are shown.4D-Var a posteriori emissions deviating more than 20% from the synthesis inversion estimate have been printed in bold face.The same holds for uncertainty reductions deviating more than 0.1 from the synthesis inversion estimate.
(Bergamaschi et al., 2005)mpared with modelled 3-hourly mean concentrations.The measurement error is assumed to be 3 ppb, to which an estimate of the representativeness error is added, based on the modelled concentration differences with neighbouring grid cells(Bergamaschi et al., 2005).This yields values smaller than 1 ppb for the Antarctic stations, and typical values of 10 to 15 ppb for most other stations, with occasional increases to several tens of ppb's at stations www.atmos-chem-phys.net/8/6341/2008/Atmos.Chem.Phys., 8, 6341-6353, 2008