Atmospheric Chemistry and Physics Discussions Interactive comment on “ Towards better error statistics for atmospheric inversions of methane surface fluxes ”

Abstract. We adapt general statistical methods to estimate the optimal error covariance matrices in a regional inversion system inferring methane surface emissions from atmospheric concentrations. Using a minimal set of physical hypotheses on the patterns of errors, we compute a guess of the error statistics that is optimal in regard to objective statistical criteria for the specific inversion system. With this very general approach applied to a real-data case, we recover sources of errors in the observations and in the prior state of the system that are consistent with expert knowledge while inferred from objective criteria and with affordable computation costs. By not assuming any specific error patterns, our results depict the variability and the inter-dependency of errors induced by complex factors such as the misrepresentation of the observations in the transport model or the inability of the model to reproduce well the situations of steep gradients of concentrations. Situations with probable significant biases (e.g., during the night when vertical mixing is ill-represented by the transport model) can also be diagnosed by our methods in order to point at necessary improvement in a model. By additionally analysing the sensitivity of the inversion to each observation, guidelines to enhance data selection in regional inversions are also proposed. We applied our method to a recent significant accidental methane release from an offshore platform in the North Sea and found methane fluxes of the same magnitude than what was officially declared.


Introduction
Quantifying the methane (CH 4 ) fluxes between the surface and the atmosphere, establishing their temporal variability and spatial distribution, and estimating the anthropogenic and natural contributions to these fluxes is critical for closing the present-day methane budget. One of the approaches used for this purpose, called the atmospheric inversion, assimilates information about atmospheric composition to infer surface fluxes. This type of top-down estimation relies on the assimilation of in-situ observations of atmospheric concentrations (Houweling et al., 1999(Houweling et al., , 2006Hein et al., 1997;Pison et al., 2009;Bousquet et al., 2011Bousquet et al., , 2006Bergamaschi et al., 2005) and/or of remote-sensing data from satellite-based instruments (e.g., Bergamaschi et al., 2009). Using observations for inversions at the global scale reduces the uncertainties on the mean CH 4 flux balances on large regions (typically a few millions of km 2 large). At the regional and mesoscales, high-resolution inversions potentially provide the spatial distribution of the fluxes, so that the characterisation of the processes involved can be improved (Bergamaschi et al., 2010).

A. Berchet et al.: Error statistics for atmospheric inversion
Inversions at any scale depend on simulations of the atmospheric mixing and advection by Chemistry-Transport Models (CTMs) to estimate the influence of emissions and sinks on the atmospheric concentrations where they are measured. Whether they are based on coarse (Chen and Prinn, 2006;Hein et al., 1997), varying  or high (Lauvaux et al., 2008;Sarrat et al., 2007) resolutions, all the CTMs suffer to a certain extent from uncertainties in reproducing the atmospheric concentrations. The uncertainties are due to the transport errors (Baker et al., 2006;Geels et al., 2007;Peylin et al., 2002;Ahmadov et al., 2007;Prather et al., 2008), to the assumption that a point observation can be compared to the mean simulated concentration on the corresponding grid box, i.e., the representation errors (Gerbig et al., 2003;Tolk et al., 2008), or to the errors from aggregating the fluxes on large regions (Kaminski et al., 2001).
In the framework of Bayesian atmospheric inversion (Enting et al., 1993;Tarantola, 1987), the implementation of a system requires obtaining an advanced understanding of the statistics of the observational and instrumental errors, the transport errors, the representation errors, and the errors of the prior distribution and magnitude of the fluxes prescribed in the system. Most of the cited works empirically assigned these error statistics. Objective methods of tuning the errors in the system also exist (Wahba et al., 1994;Dee, 1995;Desroziers and Ivanov, 2001) and have been applied to get the general structure of the errors (Michalak et al., 2005;Winiarek et al., 2012). But these methods rely on subjective prior knowledge on the error structure (e.g., isotropic spatial correlation or temporal decay in the correlations), which can limit the generality of the results.
In this study, we apply three different methods based on the statistical and algebraic properties of the errors, but with a minimum of additional physical assumptions on the error patterns. In inversion systems typically solving fluxes at the model resolution (e.g., ∼ 0.5 • ×0.5 • each week during a season or a year in regional scale studies), this approach would require the handling of matrices of error covariances the total size of which exceeds billions or even trillions of components. To embrace memory limitations and reduce the computation costs, we have chosen a short window of inversion and have aggregated the surface fluxes on synoptic-scale regions. This simplification allows applying powerful generic methods, but induces limitations (Kaminski et al., 2001;Bocquet et al., 2011) that have to be taken into account when moving to a full-resolution inversion system.
Our study exploits a recent unexpected release of CH 4 in the North Sea in spring 2012 to apply this statistical approach and test the ability of a European network of atmospheric observations to detect the leak. On 25 March 2012, an offshore oil platform on the Elgin field, located 200 km east of Scotland shores (57 • N, 1.53 • E), was evacuated due to a gas leak. The company operating the platform gave a rough evaluation of the flux reaching 200 000 m 3 d −1 or 140 metric tons per day (t d −1 ) for CH 4 , which accounts for less than 1 % of the daily regional emissions (within a radius of ∼ 750 km around the leak point) according to the Emission Database for Global Atmospheric Research (EDGAR v4.2;http://edgar.jrc.ec.europa.eu) for the year 2008. The leak was stopped two months after. The methane plume emitted by this point source is difficult to extract from the observation noise and from the variability of the other sources, which makes the assignment of error statistics particularly critical (Winiarek et al., 2012). We develop and apply a regional inversion framework based on CHIMERE CTM simulations ) on a domain covering the European continent (Fig. 1). Relying on objective statistical criteria, we optimise the covariance matrices of the errors of the observations and of the prior state vector (surface fluxes, initial and lateral boundary conditions) for two independent time windows of inversion: the 2 weeks before the beginning of the leak and the 2 weeks after. The three methods of optimisation are implemented with acceptable computation times. They managed to produce error covariance matrices, which are specifically suited to the system and the inversion window, and that are a best guess of the optimum in regard to the objective criteria of each method. Complex error structures are then retrieved. And within the framework of the underlying assumptions, every piece of information provided by the observations and the prescribed fluxes is entirely used. We then use the computed matrices to invert the European fluxes before and after the leak start and test whether or not the atmospheric network detected this methane plume.
In Sect. 2, we describe the inversion methods and the dataset used in the study. We also present the algorithms that we implemented following the literature to specify the inversion system configuration. In Sect. 3, the results of these tuning methods are presented. The inversion results from these sets are analysed in Sect. 4 and their limitations and possible adaptation to larger systems are discussed in Sect. 5.

Theory: analytical framework
We apply classical data assimilation methods based on the Bayesian formalism (Courtier et al., 1994;Enting et al., 1993Enting et al., , 1995Tarantola, 1987). In the following we use the unified notation by Ide et al. (1997). Assuming a Gaussian nature for all the errors, the method basically relies on the minimisation of the cost function:  quet et al., 2011) and EDGAR emissions during 10 days before the period of inversion; and (4) an offset, constant and uniform along the whole domain; the prior offset was calculated from the available observations as an estimation of the background concentrations during the window of inversion; the initial and boundary concentrations are expressed as perturbations from this offset.
This simplification implicitly implies the hypothesis of pure correlation of the information within each aggregated region of the state vector (see Sect. 5.1). We chose an extended domain compared to the network coverage in order to cope with the spatial and temporal ill-representation of the LBC (Lauvaux et al., 2012).

Atmospheric transport model
We use the Eulerian mesoscale non-hydrostatic chemistry transport model CHIMERE for this study . This model was developed in a framework of pollution simulations (Schmidt et al., 2001;Pison et al., 2007), but is also used for greenhouse gas studies (Broquet et al., 2011). We use here a regular horizontal grid of 50 km-side cells with 25 layers geometrically spaced from the surface to 450 hPa (∼ 6000 m). The model time step varies dynamically from 4 to 6 min depending on the maximum wind speed in the domain. The model is an off-line model which needs meteorological fields as forcing. The forcing fields are deduced from interpolated meteorological fields from the European Centre for Medium-range Weather Forecast (ECMWF) with a horizontal resolution of 0.5 • × 0.5 • every 3 h. The model is operated in a domain of limited area spanning over the whole continental Europe (roughly 24 × 10 6 km 2 ; see Fig. 1).

Observations
The study is based on the assimilation of measurements of the atmospheric composition. Concentrations of CH 4 are measured in-situ in 13 European sites (see Fig. 1; details in Table 1) at altitudes from sea level up to 3580 m a.s.l. and with different instrumentation and time resolution. Some stations are equipped with CRDS analysers (frequency of up to 1 Hz magnitude), whereas others are Gas Chromatographs. Hourly aggregates were used as input in y 0 for the inversion. The observation sites can be split into three categories: (1) mountain sites which monitor almost all the time the free troposphere; they are scarcely influenced by the local emissions and are representative of the continental and global budget; (2) coastal sites with primary influence from the ocean when the air flows towards lands; as mountain sites, they are representative of global patterns; and (3) rural sites, inland but remote from anthropogenic emissions hot spots. All instruments are calibrated by tanks traceable to the NOAA 2004 CH 4 scale  with a calibration precision of ±2 ppb.

Error configuration: description of the algorithms
In order to apply the Bayesian inversion framework, a perfect knowledge of the background and observation error statistics is needed. The tuple of covariance matrices (R, B) must then be established. Tuning and calculating optimal covariance matrices has long been of interest in data assimilation (e.g., Talagrand, 1998;Desroziers and Ivanov, 2001;Chapnik et al., 2004). Statistical studies on large sets of data are required to reach a sufficient threshold of information to get a reliable approximation of R and B. In most cases, the sets of data are not available and the covariance matrices are built relying on physical considerations and an expertise on the observation and model behaviours (Bergamaschi et al., 2010). In this section, we describe different objective methods to infer the best tuple of R and B matrices: first, the Desroziers' scheme, second, the maximisation of the likelihood, third, observation space diagnostics. The Desroziers' scheme and the maximisation of the likelihood are computed on the subspace of the diagonal matrices for both R and B, while the observation space diagnostics allow the recovery of full matrices.

Validation test: χ 2 distribution
It can be shown, within Gaussian assumptions, that for the state vectorx a minimising the cost function J , J (x a ) = J o (x a ) + J b (x a ) has the statistics of a χ 2 distribution with a mean equal to d/2, d being the total number of available observations. We then define a χ 2 index 2J (x a ) d that shall be close to 1.
The index can be written: Though insufficient, this test provides a low-cost insurance that we got a well-defined tuple of covariance matrices. The three methods presented below are tested and validated in regard to this test. A deeper analysis of the algorithms' results is presented in Sect. 3.

Desroziers' scheme: subsets application
We describe here a method to roughly infer the shapes of R and B covariance matrices with a very low computation cost. It has been shown by Talagrand (1998) and Desroziers and Ivanov (2001) that for a given (x b , y 0 ), for any subspace j Table 1. Site characteristics. The altitudes of the sites are given as m above sea level (a.s.l.) and the inlet height is in m above ground level (a.g.l.). The sites are grouped into three categories relatively to the topography and their close environment: rural (R), mountain (M), coastal (C). * These sites are recent and still do not have related publications. 1 Observatoire Pérenne de l'Environnement.  Vermeulen et al. (2005) independent of the complementary space in the observation or state space, the optimal tuple of matrices (R, B) follows the expression below: where tr(.) stands for the trace operator, is the expectation of the contribution at the maximumx a to the cost function J of the independent subspace j of the observation (resp. state) space, P j is a projector from the whole observation (resp. state) space to the subspace j and p j (resp. n j ) the dimension of the subspace j . I n stands for the identity matrix in the background space.
One cannot ascertain statistical independence between subspaces of the observation or of the state space before running the algorithm. To apply the method, we then make the following assumptions to divide the observation and the state spaces into 41 independent subspaces: (1) each measurement site is independent from all others, (2) at each site, all the observations during the afternoon (12:00-17:00) are gathered in one subspace (named "day" period in the following) and all the remaining observations during the morning and the night in a second subspace (named "night"); the planetary boundary layer (PBL) during the night is ill-represented by models; ensuing erroneous vertical mixing is expected to deteriorate the ability of the CTM to simulate realistic concentrations during the night; the shared cause of the enhanced errors jus-tifies the assumed dependency of night-time observations; for similar reasons, we group the remaining observations during the afternoon with well-mixed PBL. (3) the LBC are independent from other state dimensions, (4) same with the IC, (5) same with the offset, (6) every aggregated region of emissions is independent from the others.
Doing so, we have 26 independent subspaces within the observation space (13 sites x day/night) and 15 for the state space (12 regions + LBC + IC + offset). Desroziers and Ivanov (2001) proposed an iterative tuning procedure that converges to a tuple that satisfies Eq. (5); we refer to this procedure as Desroziers' scheme (DS). Let us rewrite the cost function: where s b j,k and s o j,k denote the adapted weights for the subspace j at step k of the iterative procedure to balance the observations and the background in the cost function.
Desroziers' scheme is described by the following system of equations for every step k: www.atmos-chem-phys.net/13/7115/2013/ Atmos. Chem. Phys., 13, 7115-7132, 2013 7120 A. Berchet et al.: Error statistics for atmospheric inversion B j and R j stand for the diagonal sub-matrices associated to the subset j . K k is the Kalman gain matrix calculated with (R k , B k ) (see Eq. 2). The method implicitly relies on the χ 2 distribution of the cost function. The computed tuples then converges to a tuple filling the χ 2 criterion. We start this algorithm from a tuple (R 0 , B 0 ) with no physical assumption, i.e., R 0 = I d and B 0 = I n . We then stop the iterative scheme when every subset contribution to the cost is less than 1 % away from its theoretical expected value (right member of Eq. 5); this is equivalent to get a χ 2 test greater than 99 %. The algorithm converges in no more than 15 steps (i.e., a couple of minutes on a standard office computer).

Maximum of likelihood
Desroziers' scheme (DS) relies on coarse approximations and, for example, cannot extract the variability of the observational errors day-by-day and hour-by-hour. The following method allows the computation of a tuple (R, B) which is tuned not by block but by component individually. This improvement implies computational cost drastically higher than for DS method, but still affordable (less than a day) in our case study. In Gaussian assumptions, the likelihood of the observations y 0 for given R and B can be written as follows (Michalak et al., 2005): The function diverges to infinity when R+HBH T ∼ (y 0 − Hx b )(y 0 − H x b ) T . But when supposing that R and B are diagonal definite positive, S cannot be written as a matrix R + HBH T . Hence, in these assumptions, one can prove that the function is bounded and admits a computable maximum (Burg et al., 1982). A proper (R, B) tuple for the inversion system is necessarily a maximum of the function (Dee, 1995).
For memory limitation reasons, we do not maximise the function itself, but equivalently its logarithm: and C a constant not relevant for computing the maximum of the function. |.| stands for the determinant operator. The function maximum cannot be easily computed analytically in general. Hence, we use an ascending pseudo-Newtonian method based on the calculation of the gradient of the log-likelihood. The algorithm converges to a local maximum (Chapnik et al., 2004), but we have no insurance of converging to the global maximum. The result of the algorithm can be very dependent from the (R, B) tuple chosen as a starting point. To ensure the robustness of the result, we test this method with 2 different starting tuples: (1) one constructed relying on expert knowledge (e.g., diagonal background matrix with variances consistent with inventories specification), (2) the other with uniform errors of 50 ppb for observations, the LBC and the offset and 10 % for the emissions. The results are very similar: the difference between the two errors related to an observation j does not exceed 5 % and is less than 1 % in average.
Additionally, to accelerate the convergence of the algorithm, one can notice (details in Burg et al., 1982) that the log-likelihood is maximum only if (R, B) satisfies: Consequently, we force α to stay close to 1 at each step of the algorithm by rescaling (R, B) by α. This normalisation is equivalent to the χ 2 test; hence the constrained maximum of likelihood algorithm necessarily fulfills the χ 2 test. This method makes the observation and background error variances statistically consistent with the prior difference between y 0 and Hx b . When assuming that R and B are diagonal, a unique set of variances matches this criterion on y 0 − Hx b . Hence, a maximum total amount of d pieces of information is used in the algorithm, while more are available in the background and the observation operator for the subsequent computation ofx a . A more precise and complete quantification of the balance between the information used for the optimisation of the matrices and for the inversion itself is difficult in a real case study. A dedicated OSSE (Observing System Simulation Experiment) could improve our knowledge in this direction.

Observation space diagnostics
With the two previous methods, R and B are confined to the sub-space of diagonal matrices. But the errors on the observations are known to be correlated through the H operator errors amongst others. Errors on the background are also correlated, mainly because of shared errors in the inventory methods and in flux process modelling. We carry out in this section an algorithm to produce a setup of non-diagonal matrices R and B.
Inquiring into error correlations requires a huge amount of information. The available information is not sufficient to characterise deterministically the full non-diagonal matrices R and B. This section must then be seen as a way to infer guidelines for covariance building. Desroziers et al. (2005) showed that the innovation vectorŝ Equation (8) is valid if and only if (1) the observation operator is linear, (2) the errors are Gaussian and unbiased, and (3) the matrix HBH T (R + HBH T ) −1 is consistent with the statistics of the background and observation errors.

Design on R
In order to build a matrix R consistent with Eq. (8), we test an iterative algorithm similar to Desroziers' scheme using the following instruction for every k: The tuple is normalised by the associated χ 2 test in order to constrain better the algorithm.
The expectation could be calculated explicitly as a combination of R k , B k and H. However, as suggested by Desroziers et al. (2005), the exact iterative scheme is not expected to converge. In our case, we indeed found no convergence. Therefore, we use a Monte Carlo evaluation of the expectation with 50 000 perturbations of y 0 and x b with Gaussian distribution of covariances R k and B k . The tuples of covariance matrices generated with the Monte-Carlo algorithm have likelihood values bigger (hence closer to the maximum of likelihood) that the one inferred from the explicit algorithm. Their likelihood (−480) is also significantly higher than with the diagonal tuple calculated in Sect. 2.2.3 (∼ −15 000). The non-diagonal tuple is then more 'likely' in the sense of Sect. 2.2.3 compared with the diagonal one.
Any inconsistency in the requirements of Eq. (8) will make the expectation non-symmetric positive definite. In particular, biases in the observations carry information that is considered by the algorithm as random errors. The built matrices are then partly influenced by the intrinsic biases in the system. However, a simple diagnostic on the spectrum shows that after 50 000 perturbations, only less than 0.1 % of the Eigen values are negative. The biases that create asymmetry in the matrix then do not seem to significantly impact our case. This could not always apply to other frameworks. The limitation of the algorithm is to build up statistics from a single occurrence ofŷ 0 andx b and therefore is not fully generic.
But since asymmetry does not seem to be critical in our case, we rebuild a symmetric semi-definite positive matrix by correcting R ′ k+1 spectrum. We start the algorithm from two different tuples of matrices: the one which maximises the log-likelihood in the diagonal assumption, and another with the same B but with R = I d . The convergence is slower with the second starting tuple but the two optimum tuples are similar.
Calculating a matrix expectation of dimension 5000×5000 with only 50 000 Monte-Carlo perturbations is very unstable. Additionally, the perturbations on the observations are generated from y 0 , whereas it should be computed from the unknown perfect unperturbed observation vectorỹ. Hence, the closer R k gets to the optimal value, the weaker is the approx-imation based on y 0 . For this reason, after a few steps of improvement of the log-likelihood value, the algorithm quickly diverges. We then keep the last step before the beginning of the divergence as a guess for the true covariance matrices tuple. Desroziers et al. (2005) also reported an equation constraining B in the observation space:

Design on B
with the innovation vectord a b = H(x a −x b ). The expectation is again based on a Monte-Carlo estimation. We compute an iterative scheme similar to that developed in Eq. (9) with also a normalisation by the χ 2 test. We use the tuple calculated above as a starting point.
Equation (10) relies on the same assumptions as Eq. (8). In particular, bias and other mismatches between HBH T (R+ HBH T ) −1 and the correct statistics can induce inconsistent asymmetries that must be corrected. Another critical point in building B from constraints in the observation space is that the computed expectation could be outside the ensemble {HBH T / B symmetric definite positive}. We then project the expectation onto this structure to recover a compatible B matrix.
Because of the Monte Carlo and the indirect estimation of B, the instability is even sharper than for the computation of R. We then do a unique iteration of the algorithm to get an evaluation of the potential optimal correlation coefficients in B.

Results
We run the three algorithms described in Sect. 2.2 to infer a best guess for the optimum tuple (R, B) of the error covariance matrices. In the following section, we describe the shape of the calculated matrices in regard to known physical patterns of errors. We will refer to the diagonal tuple computed from the Desroziers' Scheme (resp. the Maximum of log-likelihood) as (R DS , B DS ) (resp. (R ML , B ML )); the non-Diagonal tuple from the observation space algorithm will be referred to as (R ND , B ND ).

Patterns in the error variances for the 3 methods
In Fig. 2, the variances of the observation errors in the 3 methods are compared for the period after the leak start. For the ML and ND algorithms, the variances of the observation errors are averaged along the same subspaces as in DS method (see Sect. 2.2.2) to be comparable. Most of the observation errors remain within the same interval (5-20 ppb) with the 3 methods. Their magnitudes are comparable with other But the two parts of the North Sea are expected to be illdistinguished by the system, as confirmed by a posterior correlation coefficient r of −0.78 between the errors in the two regions after the leak start. Then the flux that can be attributed to the leak is defined as the difference of the emissions after and before the leak start over the whole North Sea area. Our inversion computes a flux of +406 ± 33 t d −1 . The inversion detects an increase of emissions from the direction of the leak, but fails to unambiguously affect the increased flux to the proper region. The figure we compute is of the same magnitude (3 times higher) than the estimation given by the operator and does not exceed 18 % of the background emissions related to oil and gas extraction in the North Sea.
Our results also reveal a high dependency to the meteorological situation during an inversion window. The ratio between the increment and the posterior errors on the emission budget in the inner regions is very different for the two periods (∼ 1 before and ≫ 1 after the leak start). The reconstructed error on the total budget largely depends on the correlation coefficients in P a . For the period before the leak start, most posterior correlations are large and positive (r > 0.9). The gradients are then well constrained while the total budget stays uncertain. On the opposite, after the leak start, neighbouring regions exhibit negative correlations by pair (e.g., "NSN" -"NSS", "MGP" -"GBT"). The assimilation of the observations cannot separate the contribution from these close regions, but it leads to a good reduction of the error on the total balance. These two different behaviours may be related to different synoptic regimes during each inversion window: before the leak start, an anticyclone was laying on central and western Europe; after the leak start, air masses coming from North Europe vented the domain.

Limitations and hypothesis probation
All the results depend on strong statistical and physical hypotheses, which may not all be robust. First, we show in Sect. 3.1 that the assumed Gaussian errors of the background can produce physically inconsistent inverted fluxes. Adding Lagrangian correcting factors (e.g., Göckede et al., 2010) to the cost function (Eq. 1) can ensure physically consistent fluxes. But that would alter the algebraic properties of the problem and make the implementation of our methods more complex. With regions that act as buffers against the uncertainties on the LBC, the ND methods proved to acceptably deal with the issue.
Second, all CTMs suffer from weaknesses and errors in their parameterizations and numerical scheme. The induced errors can be systematic and not only random, as suggested in Sect. 3.2. They should then be considered as a bias η in the observational errors ǫ = y 0 −Hx ∼ N (η, R ′ ). Further investigations on the effects of the parameterizations, the resolu-tion and the inputs to the CTM shall be carried out to quantify and fix as much as possible the bias η. More specifically, in Sect. 3.2, we showed that very high diagnosed error variances during the night could be related to systematic statedependent biases in the CTM vertical mixing and in PBL modelling.
Third, aggregating fluxes within bigger regions implicitly implies full correlations of the errors on the background in each region. Despite this strong assumption, our methods are supposed to diagnose the error on the aggregated fluxes that are within the footprint of the network. On the opposite, the aggregated regions that are partly within and partly outside the network footprint will exhibit strongly biased diagnosed errors and increments. Kaminski et al. (2001) studied the issue and found potential errors of the same magnitude as the fluxes themselves. A better choice of the resolution and of aggregated regions considering the prior fluxes and the transport patterns (e.g., Wu et al., 2011) during the window of inversion should significantly improve the results of the methods.
Despite these weaknesses in our methods, the optimal tuple of covariance matrices gives better results than a tuple built on expert considerations: either these expert-built tuples, which are most of the time diagonal, are similar to (R DS , B DS ) that causes inconsistent negative CH 4 fluxes, or the observation errors are enhanced to reduce their impact on the inversion; but in this latter configuration applied to our inversion windows, the corresponding inverted fluxes remain close to the prior ones and the flux uncertainties are not noticeably reduced. Our objectively calculated tuple gives better inversion results, with reduced posterior uncertainties. Moreover, some of the computed error patterns are generic and are transferable to other larger systems. In this study, we chose a particular representation of the complete fullresolution state vector. Most errors represented by the covariance matrix R are independent of this representation . As a consequence most results on R will remain valid in the framework of a full-resolution state vector. The recovery of the errors of the non-aggregated background vector are more ambiguous and only large patterns could be inferred for finer resolutions. Additional hypotheses must be made on the shape of the full-resolution background errors to deduce their values from the aggregated matrix. Our methods can then be seen as a way to simplify and project problems with large state vectors in order to infer the patterns of the errors with relatively small computation costs.

Implications for data selection
The framework we chose allows the explicit computation of the sensitivity of the inversion to each observation. We follow Cardinali et al. (2004) and calculate the influence matrix, which gives the effect of a small change of y 0 on Hx a : S = (R + HBH T ) −1 HBH T . For every observation y 0 j , high observation errors reduce the contribution of the observation morning. One should expect a better confidence in late afternoon results from the CTM, when the vertical mixing is still active (though reduced), than early morning.
We emphasised in Sect. 3.2 the diurnal patterns of the errors on the observations. Further efforts have to be made in modelling the PBL height and the vertical mixing to ensure better quality and reduced bias in the simulations at the end of the night. Our method would then allow a better use of the observations for local and regional inversions.

Conclusions
We inquired into the possibility of precisely and objectively estimating the covariance matrices of the errors on the observations and the background (R and B) that best fit inversion system requirements. A best guess of these matrices with regard to objective criteria is needed in the Bayesian inversion framework, especially for regional studies. To do so, we used algorithms developed in a theoretical framework, but too complex to be tested in full-resolution systems. The translation to a regional configuration was carried out by simplifying the system and reducing the total size of the covariance matrices to allow an algorithmic tuning of R and B that estimates optimum tuples (R, B) in terms of statistical properties.
We tested 3 algorithms of growing complexity (and computational costs) to estimate the optimal tuple (R, B). Unlike other studies, which make strong physical assumptions, such as isotropic spatial correlations in the observation errors or temporal decay of the correlations, we minimised the number of assumptions to keep a better objectivity in our results. In principle, all the patterns of errors can then be recovered specifically to the system during the window of the inversion we focused on.
Amongst other noteworthy patterns, our algorithms retrieved the errors due to the mis-estimation of the planetary boundary layer height in global circulation models, i.e., large errors during the night and lower during the day when the CTM reproduces the atmospheric transport better. This source of errors contribute in a large part (50 % in most plain sites) to the diurnal variability of the observation errors and also causes significant temporal correlations within a 24 h period.
Additionally, our approach does not require a prior filtering of the observations we could consider as ill-simulated by the model. In theory, all the available observations can then be assimilated in the inversion system and not only the ones during early afternoon. However, our study points at probable significant systematic bias in the CTM during the night. The night observations should then be excluded. But the algorithms give objective tools to diagnose the need for efforts to better simulate the atmospheric behaviour during the late night when the observations seem to have the biggest impact on the inversion results. Late afternoon observations were also computed to have a significant influence on the inversion results. A cautious implementation of these observations into an inversion system is expected to enhance the efficiency of the system.
The prospects from this work will be to quantify the uncertainties in our methods and their impact on the optimised fluxes. A dedicated Observing System Simulation Experiment could be carried out in that sense. The computational costs should also be reduced by running our scripts in parallel. In the framework of an inversion system with full temporal and spatial resolution, when variational algorithms are necessary to compute the optimal fluxes, our general method may overburden the computer and memory capacity. Indeed, the limiting factor in our algorithms comes from the diagonal maximisation of the log-likelihood needed to compute the non-diagonal optimal tuple. The maximising algorithm induces computational costs limited by the size of the 'background' vector x b . The computational complexity is at least O [dim(x b )] 3 while full-resolution state space dimension is of several orders of magnitude larger than our reduced state space. Nevertheless, the method could be tested in systems of intermediate complexity to infer additional knowledge on the statistics of the errors.