Structure of the transport uncertainty in mesoscale inversions of CO2 sources and sinks using ensemble model simulations

We study the characteristics of a statistical ensemble of mesoscale simulations in order to estimate the model error in the simulation of CO2 concentrations. The ensemble consists of ten members and the reference simulation using the operationnal short range forecast PEARP, perturbed using the Singular Vector technique. We then used this ensemble of simulations as the initial and boundary conditions for the meso scale model (Ḿ eso-NH) simulations, which uses CO2 fluxes from the ISBA-A-gs land surface model. The final ensemble represents the model dependence to the boundary conditions, conserving the physical properties of the dynamical schemes, but excluding the intrinsic error of the model. First, the variance of our ensemble is estimated over the domain, with associated spatial and temporal correlations. Second, we extract the signal from noisy horizontal correlations, due to the limited size ensemble, using diffusion equation modelling. The computational cost of such ensemble limits the number of members (simulations) especially when running online the carbon flux and the atmospheric models. In the theory, 50 to 100 members would be required to explore the overall sensitivity of the ensemble. The present diffusion model allows us to extract a significant part of the noisy error, and makes this study feasable with a limited number of simulations. Finally, we compute the diagonal and non-diagonal terms of the observation error covariance matrix and introduced it into our CO 2 flux matrix inversion for 18 days of the 2005 intensive campaign CERES over the South West of France. Variances are based on model-data mismatch to ensure we treat model bias as well as ensemble dispersion, whereas spatial and temporal covariances are estimated with our method. The horizontal structure of the ensemble variance maniCorrespondence to: T. Lauvaux (thomas.lauvaux@lsce.ipsl.fr) fests the discontinuities of the mesoscale structures during the day, but remains locally driven during the night. On the vertical, surface layer variance shows large correlations with the upper levels in the boundary layer ( > 0.6), dropping to 0.4 with the lower levels of the free troposphere. Large temporal correlations were found during the afternoon ( > 0.5 for several hours), reduced during the night. The diffusion equation model extracted relevant error covariance signals horizontally, with reduced correlations over mountain areas and during the night over the continent. The posterior error reduction on the inverted CO 2 fluxes accounting for the model error correlations illustrates the predominance of the temporal over the spatial correlations when using tower-based CO 2 concentration observations.


Introduction
Atmospheric inversions are a widely-used tool for the quantification of surface sources for CO 2 (e.g.Gurney et al., 2002;Rayner et al., 2008), for CH 4 (Bousquet et al., 2006) and for CO (Petron et al., 2002).The theory and applications are described in Enting (2002).In the most common approach, Bayesian Synthesis Inversion, one starts with a prior distribution of surface fluxes.These are used as inputs to an atmospheric transport model (possibly including chemical modification).The transport model simulates concentrations at a set of observing locations.The fluxes are adjusted to optimize consistency with both the observed concentrations and a priori flux information.Under the usual assumption of multivariate normal probabilities, the required degree of consistency is described by covariance functions for the prior fluxes and the model-data mismatch.
The proper specification of these covariances is as important as the prior and data values themselves.In the underlying statistical formulation, the uncertainties normalise quantities such as the model-data mismatch so their specification Published by Copernicus Publications on behalf of the European Geosciences Union.1090 T. Lauvaux et al.: Atmospheric CO 2 modelling: error correlations is vital for a correct solution.In theory, the prior flux error covariances can be estimated more directly than the observation error covariances, since the error flux covariances should represent the statistics of the difference between the true fluxes and our chosen flux prior.Nevertheless, its characterization is limited by the number of flux observations for each ecosystem in the domain.This requires some independent knowledge of the truth which should, in turn, not be used in our later inversion.See Chevallier et al. (2006) for an example of using eddy covariance flux measurements as truth estimations.
This paper concerns the estimation of the model-data mismatch covariance.In the full theory (see Tarantola, 2005, ch.1.)there are separate covariances describing both the model and measurement parts.The measurement covariance describes the statistics of the difference between a measurement and the true value.As with the prior this can often be assessed with independent measurement and one has access to a well-developed discipline of metrology.The model uncertainty describes the statistics of the difference between a simulation and the true value that would be observed if the real atmosphere was forced with the same fluxes as the model.We see immediately that this error is much harder to characterize since there are few cases in which we know what fluxes influenced the real atmosphere.In atmospheric transport inversions the model contribution usually dominates the measurement error (Gerbig et al., 2003;Stephens et al., 2007;Law et al., 2008) and this model error is our focus here.
There are many contributions to the model error.Transport models are mesh-based.Thus there is likely to be a mismatch between a point measurement and the average value over a model grid box.The importance of this depends on the magnitude and structure of concentration variability in the atmosphere and we can attempt to quantify it with highresolution models (e.g.Corbin and Denning, 2006;Corbin et al., 2008) or with spatially dense measurements (e.g.Gerbig et al., 2003).Such direct comparison using radiosonde's was performed at larger scale (Gerbig et al., 2008) but requires a sufficiently large number of observations over the domain.Finite resolution equally confounds accurate description of the input fluxes leading to so-called aggregation error (Kaminski et al., 2001).Next the formulation of the model can be incorrect.In the absence of further data to demonstrate this we are forced to treat the statistics of an ensemble of models as if they are the statistics of the difference between the model and the truth.This is the assumption underlying the series of Transcom studies (e.g.Gurney et al., 2002;Baker et al., 2007;Law et al., 2008).This approach can hide bias in the model ensemble and is commonly criticized for understating model error although Stephens et al. (2007) noted that screening models against available independent data can actually reduce the dispersion of the ensemble.
Finally, even if the physical formulation of the model is close to the truth, uncertainties in the analyses with which it is constrained means that our transport fields are only one realisation of an ensemble consistent with available meteorological information.In the case of mesoscale modelling using meteorological analyses, forecast error growth (e.g.Lorenz, 1982) affecting the spread of the ensemble in addition to analysis uncertainties remains negligible.The impact on atmospheric inversions of the dispersion of this ensemble, and particularly the strong likelihood of correlations in this dispersion has not been previously investigated.A direct estimate of the ensemble statistics is fraught with difficulty since a sufficiently large ensemble is prohibitively expensive to compute.Instead, where available, we use a physical model of error growth and propagation previously applied to the numerical weather prediction problem (short range ensemble prediction system PEARP 1 ).
We apply our approach to an inversion of CO 2 surface fluxes in a limited domain although we stress in the discussion that it is applicable to global problems.We calculate the ensemble statistics of transport error on the domain of the CarboEurope Regional Experiment Strategy (CERES) (Dolman et al., 2006;Lauvaux et al., 2008).At this scale we can compare the ensemble behaviour with other measures of the model-data mismatch such as the variability in measured concentrations (Gerbig et al., 2003).Also, mesoscale inversions are heavily reliant on the ability of transport models to simulate details of concentration variations usually ignored in global inversions.The improvement of atmospheric models now allows the use of high frequency data on limited domains (Sarrat et al., 2007).Although model errors should decrease with increasing spatial and temporal resolution, the complexity of dynamical processes at the mesoscale are likely to produce a complex uncertainty structure.Discontinuities in the dynamical flow at higher resolution induced by the surface properties also suggest complex structures of the model error.
In this study, we propose to assess the structure of the model error by perturbing the synoptic conditions, that affect the mesoscale dynamical structures, using an ensemble of simulations at a larger scale.The ensemble presented in this study remains physically consistent, thanks to the twostep perturbation method that allows us to keep the original dynamical schemes of the meso scale model, here Méso-NH (Lafore et al., 1998).The sensitivity to the model parameters that affect the vertical transport for example, and is of major importance in the transport error, implies to use an ensemble of transport models (changing a parameter is similar to a different model) that limits the interpretation of the ensemble spread and the error structures.Here, we restricted ourselves to the boundary influence which is surely affecting the model we use in the inversion.Previous studies attempted to disturb the model stability by perturbing directly the modelled wind fields (Law et al., 2003), or by exploring the physical parameters of the model (Annan et al., 2005).It is impossible to tell whether a given perturbation of wind fields or parameters is physically realisable while the perturbations we use are generated from an ensemble forecast system designed to explore the plausible range of synoptic variations.
The outline of the paper is as follows: Starting from the mesoscale ensemble of simulations, (i) we estimated the variance with its spatio-temporal correlations; (ii) we modelled the horizontal correlations by using the diffusion equation to filter the noise of our limited size ensemble; and (iii) we tested the combined spatio-temporal observation covariance matrix in our CO 2 flux matrix inversion.We also discuss the implications for network design, and optimal spatial and temporal densities of measurements.

Models and diffusion equation filtering
The major task of this work is to estimate the spatial and temporal correlations in the model-data mismatch using an ensemble of regional CO 2 simulations.In this section we first describe the underlying physical models, the generation of the ensemble and finally the calculation of the various correlation terms, focusing on horizontal correlations.The role of the resulting correlation matrix in the subsequent flux inversions is treated in Lauvaux et al. (2008) and only briefly reviewed here.

Models
The ensemble is based on the global spectral ARPEGE model (Courtier et al., 1991) used with the nominal spectral truncature T358 and a stretching coefficient of 2.4, corresponding roughly to 20 km resolution over France.On the vertical axis, 41 levels describe the atmosphere from the surface to 1 hPa.The model is coupled with a four dimensional variationnal assimilation system (4DVAR) including classical meteorological observations and satellite data.The 11 ARPEGE simulations are run over 102 h, from 6 p.m. the 23 of May, to 12 p.m. the 27 of May 2005.In the Sect.3.3, we discuss the temporal evolution of the error structures over the 102 h of simulation to verify that no significant increase (or decrease) affected the error correlations.The same calculation was done on the correlation lengths in the Sect.3.4 to confirm the absence of trends.
The non-hydrostatic atmospheric mesoscale model Méso-NH (Lafore et al., 1998) was used to simulate the atmospheric dynamics during the same period (23 to 27 of May 2005), over the limited domain of the CERES campaign, a 300*300 km 2 of South western France, (Dolman et al., 2006) at 8 km resolution.65 vertical levels describe the atmospheric column up to 13 km.The boundary conditions from the ARPEGE ensemble simulations are coupled each 3 h to constrain the meso scale model, following the Sommerfield equation for the normal wind velocity components at the boundaries with a constant phase speed (relaxation term) of 20 m s −1 .The reference simulation and the ten different members run over 102 h, on the same period as the ARPEGE simulations.
The land surface scheme ISBA-A-gs (Calvet et al., 1998) was coupled on-line to the atmospheric model Méso-NH to simulate the surface fluxes for water, heat, and CO 2 (Sarrat et al., 2007(Sarrat et al., , 2008)).The land cover map describes the vegetation cover at 250 m resolution from the ECOCLIMAP database to calculate the biogenic CO 2 surface fluxes.The anthropogenic emissions are prescribed by the 10 km resolution inventory from the University of Stuttgart (IER), and the air-sea fluxes by Takahashi et al. (1997).The three-week simulation was run over the CERES domain of 61 points east-west by 51 points north-south at 8 km resolution.The interaction at each timestep between the atmospheric model and the surface scheme implies that the surface energy fluxes are slightly different in each member, amplifying the differences in boundary layer dynamics.Nevertheless the initial ground water content remains identical for all the different simulations limiting the surface flux differences for the different members.
The mesoscale inversion we used to estimate the impact of the observation error structure requires the computation of the adjoint transport.This tracer backward transport was simulated here by the Lagrangian Particle Dispersion Model (LPDM) described by Uliasz (1994).Particles are released from the receptors in a "backward in time" mode with the wind fields generated by the eulerian model Méso-NH.The dynamical fields in LPDM are forced by mean winds (u,v, w), potential temperature, and turbulent kinetic energy (TKE).

Generation of the ensemble
The perturbations represented here by singular vectors are added to the large scale simulations with the global spectral model ARPEGE (Courtier et al., 1991) over the whole period and the resulting ARPEGE fields used as boundary conditions for Méso-NH.The optimal perturbations are defined by the maximisation of the following ratio Our approach to generating the ensemble maintains the physical and thermodynamic consistency of the coupled higher resolution model Méso-NH.Even though the parameters of the numerical scheme play a certain role in the mismatch between modeled and observed concentrations at fine scale, the atmosphere shows relevant spatial patterns depending on the meso scale structures rather than the length-scales of the perturbations.The perturbation approach has been demonstrated for the ozone air quality models (Carvalho et al., 2007).We modelled the spatial correlations on the horizontal plane based on the local structure of the variance, to filter at each grid point the noise and extract the significant correlations.This method allows then anisotropic structures on the plane.

Observation error covariance estimation by diffusion equation modelling
We separately estimate the horizontal, vertical and temporal components of the error correlation.In this section, we detail the methodology of Pannekoucke and Massart (2008) applied for the horizontal component only.The structural properties of the ensemble variability are estimated by using a diffusion operator based on a local diffusion tensor.Our starting point is the 11 (1 reference and 10 perturbations) CO 2 fields from MésoNH at 1 h resolution.The estimation of the correlation tensor C is computed in the model space from the ensemble, then mapped to the observation space by the operator P , and weighted by the variances , estimated independently by data comparison, as follows: As the ensemble is of finite size, the estimation is affected by sampling noise (Lorenc, 2003).Thus, the estimation of the local observation error has to be done carefully by extracting the signal and filtering the noise.We do this by assuming some simple model of the spread of the differences among ensemble members and estimating the parameters of this model from the raw correlations calculated from the ensemble.This is done following the work of Pannekoucke and Massart ( 2008): the local correlation function is derived through the estimation of the local anisotropy tensor ν(x) under a local Gaussian approximation, and then the correlation function is modelled around each x with the formulation based on the diffusion equation (Weaver and Courtier , 2001).In that framework, the anisotropy tensor is no more than the local diffusion tensor.The idea of such an error model derives from the particular solution of the homogeneous diffusion equation: if the local diffusion tensor ν(x) is constant over the domain (which is assumed to be the plane), then the solution of the diffusion equation with initial state η(x, t = 0) = δ x (x) where δ, the Dirac distribution, is with x the surrounding points of x, and with and | | is the determinant of .The inverse of can be written as The scales L x and L y correspond to the one-dimensional differential length-scale along the direction x and y (Daley, 1991;Pannekoucke et al., 2008).This particular heat kernel Eq. ( 4) can be seen as a Gaussian correlation function on the plane (except for the normalization term).Thus, it is possible to construct quasi-Gaussian correlation functions as the result of time integration of the diffusion equation.In the more general case, the local diffusion tensor ν(x) is non constant, and it can be approximated locally, under some local homogeneous assumption.This approximation is based on the estimation of the local matrix −1 (x) from the estimation of the local length-scales L x (x), L y (x) and 1/L xy (x).From our ensemble of simulations, we estimate the length-scales L x (x) and L y (x) at each pixel based on the gaussian approximation of the correlations for each surrounding pixel on the two axis x and y, and the terms 1/L xy (x) of the local tensor with a Gaussian function, averaged on the four diagonal directions.Then the local diffusion tensor ν(x) is computed from Eq. ( 5) as half the local tensor (for the particular time integration t = 1).
Correlation matrix modelling constructed with the diffusion equation results from the product of operator C = L 1/2 W −1 L 1/2 , where L 1/2 is the propagator that corresponds to the time integration of Eq. (3) from time t=0 to time t=1/2, and W −1 is a metric tensor.
The diagnosis of the correlation is then achieved through a randomization method (Weaver and Courtier , 2001;Fisher and Courtier, 1995): of Gaussian vectors ζ with zero mean and the idendity as a covariance matrix are diffused from time t=0 to time t=1/2, using the heterogeneous diffusion equation with the local diffusion tensor.The result of this time integration is an ensemble of error vectors The assumed Gaussian correlation functions imply that the estimated correlations are positive.However, with the dynamics and the sampling noise, negative correlations can occur.In such a case, the points associated with negative correlations are not taken into account.Negative correlations affect less than 10% of the estimated correlations.They occur depending on the size of the ensemble, which suggests Biogeosciences, 6, 1089-1102, 2009 www.biogeosciences.net/6/1089/2009/that they are mainly associated with sampling noise.However, some negative correlations could be attributed to the dynamics.The original domain was duplicated two times, equivalent to four symetric domains, to avoid boundary issues.
Note that the covariance tensor C modelled with the diffusion equation has to be normalized in order to obtain a correlation tensor C= C T , where is the diagonal matrix of inverse standard deviations of C.This is a reappearance of the normalization term 2π| | 1/2 that occurs in Eq. ( 4).In the heterogeneous case, the normalization is not constant but it can be either estimated through the randomization method, or it can be approximated by the local normalization as 2 (x) = 2π| (x)| 1/2 (when this estimation is accurate enough).Finally, the above development generates the horizontal correlation for all points in the domain.We need simply to choose points corresponding to the locations of measurements to generate the observation correlation matrix we need.

Vertical and temporal correlations
Unfortunately, there is no theory for vertical and temporal error correlations as well developed as that used in the previous section for horizontal correlation.We therefore calculate vertical and temporal correlations directly using the 11 ensemble members at our disposal.These correlations are hence more subject to sampling noise than the horizontal.

Pseudo-flux inversion and spatial error correlations
The models used in this study and the inversion framework that illustrates the impact of the filtered correlation lengths are described in Lauvaux et al. (2008).Empirical analysis will be based on the error reduction Err that depends only on the observation and the prior flux uncertainties, but not on the prior CO 2 flux nor concentration data: with σ the error associated with B the background error covariance matrix, and to A the posterior flux error covariance matrix estimated by the following equation: with H the jacobian of the transport matrix, and R the observation error covariance matrix.The error reduction presented here corresponds to the optimization of the 6-day flux averages over 18 days for nighttime and daytime.Due to the limited number of observations (two towers observing hourly CO 2 concentrations, ie 2 * 144 observations) and the large number of unknowns (2 * 61 * 51), the 6-day average inversion ameliorates the problem of under-constraint but implies complete temporal correlation of the fluxes over the period, consistent with the study of Chevallier et al. (2006).
The covariance matrix of model-data mismatch R combines spatial and temporal correlations that we extract separately from the ensemble simulations.The spatial correlation between the two tower locations is computed after the diffusion modelling step described in the Sect.2.3, using 400 independent realizations.

Combining temporal and spatial correlations
We calculate temporal correlations from the first three days of the simulation period, averaged over the continental part of the domain.At each hour h, the temporal correlation is estimated with each following hour t + n with 1 < n < 24.
For each hour, results are averaged for the three possible days to produce an average diurnal cycle.
The correlation matrices, C t for the temporal component and C s for the spatial component, are combined into C the correlation matrix in the observation space via: This implies similar weights for the two components, which seems reasonable without extra information.The final R observation error covariance matrix is then constructed as defined in the Eq. ( 2) by R = C T with the diagonal matrix of standard deviations estimated from the model-data mismatch.

Results
In this section we first study the structure and magnitude of the covariances derived from our ensemble.Then we investigate the effect of these, separately and in combination, on an inversion.

Atmospheric CO 2 variability over the domain
Considering the 11 simulations at the two measurement locations of the 2005 CERES campaign (Biscarosse and Marmande towers), we estimated the variability of the different members compared to the observed CO 2 concentrations.In Fig. 1b, at the Marmande site, the ensemble spread shows a large diurnal cycle corresponding to different stability conditions during the nights between the different members of the ensemble, and similar modelled CO 2 concentrations due to well-developed convective boundary layers during the day.Maximum values of CO 2 during the night correspond to larger spread of the ensemble, but the largest modeled concentrations are smaller than the observed peaks.
At the Biscarosse location (Fig. 1a), the first two days show little spread in diurnal variability, but the last two days correspond to a large range of simulated concentrations, without any diurnal cycle.Figure 2a shows the spatial pattern of the larger CO 2 variability region including the location of the Biscarosse tower.This well-defined structure mirrors    During the day, the spatial structure of CO 2 variance is induced by meso scale processes like the changing of the direction of the "Autan wind" over the Mediterranean sea.During the night (Fig. 2b), the CO 2 spatial patterns reflect the topography of the river valleys or mountain ranges in addition to the local flux variability.The nighttime standard deviations are larger compared to the daytime, up to 20 ppm over the continent (9 ppm during the day).

Vertical structure of the error covariances
For each level we computed the ensemble variability and the correlation with the lowest level.Both show clear diurnal patterns (Fig. 3).During daytime, the ensemble standard deviation of CO 2 decreases linearly up to 1500 m from 2 ppm to 1 ppm, then remains constant until 3000 m, and decreases again in the mid troposphere.The associated correlation coefficients show a similar vertical profile, decreasing linearly up to 1000 m from 1 to 0.5-0.6,then constant until 3000 m about 0.  in the first 500 m from 5 ppm to 1.5 ppm.Between 1000 m to 3000 m, the profile is constant around 1 ppm, and then decreases in the mid troposphere.Contrary to the daytime profiles, the vertical correlations decrease from 1 to 0.6 in less than 200 m, and remain constant (0.4-0.5) up to 3000 m.Compared to the nighttime standard deviation, the shape of the correlation profile suggests some residual variability in the lower troposphere from the daytime, but the surface layer is mostly influenced by the CO 2 surface flux variability and not by the lower troposphere, which is more correlated during daytime.The implication of such correlations for future inversions will be discussed in Sect. 4.

Temporal error correlations
We estimated the temporal evolution of the error structures during the 102h of simulation to avoid any trend due to the perturbation method.During the four days of simulation, temporal correlations for each day were similar in spite of ponctual differences due to specific daily conditions.Then we estimated the spatially averaged temporal correlation with the next 24 h over the whole domain.Figure 4 illustrates the temporal correlations between different reference hours with the next 24 h, starting from midday to 9 p.m. in Fig. 4a, and from midnight to 9 a.m. in Fig. 4b.The two figures show basically similar 24-h variations, with negative correlation values between the day and the night.The case starting at midday (Fig. 4a) shows non-significant correlation with the rest of the afternoon (about 0.2) whereas the later starting cases show correlation values greater than 0.5.For correlations with the next day, the strongest correlation reached 0.4, for the 9 p.m. starting case.The nighttime cases (Fig. 4b) decrease to negative values faster, in less than five hours for midnight and 1 a.m. starting cases.The nocturnal stability is reached later, from 4 a.m. to 9 a.m.cases showing larger values during the first hours (more than 0.4).As with the daytime cases, the next night shows small values of correlations (maximum 0.4).Even though nighttime error structures are inferred by more static parameters (orography, surface properties), the sensitivity to perturbations is much larger and implies noticeable differences from night-to-night induced by small variations of the conditions.

Estimation of the correlation lengths over the domain
The diurnal variability of the ensemble variance gives rise to two separate sub-ensembles corresponding to the diurnal length-scales (from midday to 9 p.m.), and the nocturnal length-scales (from 10 p.m. to 7 a.m.).As for the temporal correlations, we compared the estimated correlation lengths for each day and verified that no trend affected the error structures, showing the independence of the structures to the perturbations.We then modelled the spatial lengthscales of the two sub-ensembles using the diffusion equation described in Sect.2.3.The raw diurnal length-scales (Fig. 5a) show high values over the Mediterranean sea and the Atlantic ocean, up to 150 km, but no clear signature over the continent.After the diffusion modelling (Fig. 5b), the low length scale values map the mountain areas (Pyrenees, and Massif Central) more clearly.Higher values over the sea and the ocean remain but low values along the sea shores were filtered.In the Golfe du Lyon (between the Eastern Pyrenees to the Rhone river Estuary), the length-scales originally about 50 km decrease to 30 km.On the average, the highest values over the continent are smaller after the diffusion modelling, decreasing from 80 km to 60 km.For the nocturnal length-scales, averaged initial values are smaller than nocturnal values (Fig. 5c), about 20 to 30 km over the continent, whereas oceanic length-scales are similar to diurnal values (Fig. 5d).We notice similar effects concerning mountain areas with lower values after diffusion modelling (Pyrenees).The noise over the continent is reduced after the diffusion modelling, between 0 to 10 km.On the contrary, to the North of Bordeaux, length-scales are higher than initially.We compute the correlation coefficients between Marmande and Biscarosse towers during the day and the night.The diurnal correlation is about 0.14 and the nocturnal correlation 0.05, corresponding to the length-scales of 40 to 50 km during the day, and 20 to 30 km during the night in Fig. 5.

Impact of the modelled covariances on the inversion
After the diffusion modelling of the length-scales and the estimation of the temporal correlations, we construct the error covariance matrix R in two steps.simulations, the operational ensemble forecast perturbs the critical regions controlling the synoptic conditions.The under estimation of the variance remains but the covariances are well-structured as shown in figure 5 after diffusion modelling.
The error correlation length-scale used in the paper does not correspond to the classical definition used in trace gas inversions, i.e. the distance at which the correlation reaches 1/e.Here we use the differential length-scale (Daley , 1991, p110), related to the curvature of the correlation function at the origin.The estimation of this length-scale is simple in the isotropic case as it is related to the power spectrum of the auto-correlation (with the help of the Wiener-Kintchine theorem).Directional approximations of this length-scale, with emphasis on there sampling distributions, exist (Pannekoucke et al., 2008).The estimation is slightly biased with an over-estimation of the truth (Pannekoucke et al., 2008).The sampling distribution has an heavy tail with a positive skewness.It results that, in average, the sampling noise leads to an estimation larger than the truth (see e.g.Pannekoucke et al. (2008), Fig. 7).That is, the smooth of the raw lengthscales modeled with the diffusion operator appeared as beneficial revealing clearly the initial true correlation structures (Pannekoucke and Massart, 2008).
Two sub-ensembles were derived from the temporal correlations, corresponding roughly to nocturnal and diurnal periods (each of them limited to a few hours).The length-scales (L x and L y ) of the correlations are then associated with elliptic functions on the horizontal plane.The ellipse is defined by its longer axis denoted L, the smaller axis denoted l, with the anisotropy 1 − l L , and the angle of the anisotropy θ.As the ensemble variations are flow dependent, preferential directions of anisotropy are expected over the domain.For example, diurnal circulation dominated by advection or local thermal gradients constrains the orientation of the largest correlations.
We have calculated the correlations using the whole period.This represents a trade-off.We need a correlation model which is stable and represents the internal dynamics of the model (suggesting longer periods of integration) but on the other hand the correlation structure really should depend on flow regime.Our approach averages over several regimes and hence averages several directions of anisotropy (fig.8 (a)), and makes the magnitude of anisotropy spatially homogeneous (fig.8 (b) ).When comparing the length-scales L y (North-South) at the tower locations, the values are smaller (less than 30km during the day), which corresponds to an East-West orientation of the anisotropy.In this region, the Atlantic ocean forces the error structure at Biscarosse, and the Garonne river valley does likewise at Marmande.
The above results have a range of impacts on current and future inversions.Recall that the errors we assign to an observation in an inversion are implicitly dependent on our ability to model that observation.We have shown that a part of such model errors shows a complex structure of correlations.In general, the existence of positive correlations in the modeldata mismatch reduces the amount of information available from those observations (Section 3.5).This is not universal however.A correlation in the model-data mismatch implies a set of preferred directions for this mismatch.Mismatches in that direction will be deweighted (reducing available informaion) however those in orthogonal directions will be highlighted.This explains the mixed patterns in Fig. 8 (b).
The spatial correlation lengths we obtained do not suggest www.biogeosciences.net/bg/0000/0001/Biogeosciences, 0000, 0001-15, 2009 ) is finally applied to the three inversion segments.On the first period, the three different matrices are applied, and compared to the diagonal matrix case.Figure 6 shows the error reduction for daytime fluxes in the period (from 24 May to 30 May) (Fig. 6a) compared to the three different correlation cases (Fig. 6b, c, and d).The "spatial only" set-up weakens the error reduction from 0% to 10%, from the original value of the error reduction (Fig. 6b).For the "temporal only" set-up (Fig. 6c), we observe the decrease of the error reduction around Marmande tower, but also increases and decreases over the rest of the domain.We discuss in Sect. 4 the origin of the positive and negative impacts.Finally, the impact of the full spatio-temporal covariance matrix is shown on Fig. 6d, highly similar to the "temporal only" case.This similarity is explained by the larger temporal correlations (up to 0.9) compared to the spatial correlations (0.05 to 0.15).
The combined correlation matrix is then used for the two other periods.Figure 7 shows the error reduction for the two periods (from 30 May to 5 June, and 5 to 10 June).We see that the third period shows a different surface response compared to the two first periods, dominated by northern winds.We suppose here that spatial correlations established over the first four day period are the same for the three periods (18 days in total), even though atmospheric dynamics changes.This assumption has limited impact considering the small spatial correlations compared to the temporal component.The two periods show the increase of the error reduction over the Atlantic Ocean, corresponding to the sea breeze circulation of the afternoons and the northern coastal winds.On the contrary, local influence around the towers shows smaller error reductions for the two later periods.

Discussion
It is likely that our estimates of by the model variability of the ensemble are too low.Our perturbations were limited to the boundary conditions for Méso-NH.Our realisations used the same internal dynamical schemes.Larger variability is expected if we explored the full space of input parameters and formulations.The existing biases peculiar to our meso scale model affect the whole ensemble of simulations.Systematic errors in the models affecting the nighttime build-up for example are due to different problems.First, the vertical resolution should be the first limitation to investigate.Second, nocturnal boundary layer conditions usually neutral are not well parameterized.An ensemble of simulations considering this issue could study its sensitivity.Finally, the surface conditions are also affecting the energy budget that drives the a large impact on most current observing networks since it is rare to have observations spaced less than 100km apart.An immediate exception is measurements from aircraft.It appears that, even if the transport model we use for the inversion allows sampling at the resolution of a few km, spatial error correlation may reduce the effective sampling density.We cannot know the true import of this until we calculate the impact of horizontal correlations for various observational densities.
A less immediate but stronger consequence is the inversion of data retrieved from upcoming satellite missions (e.g.Crisp et al., 2004;Hamazaki et al., 2004).The spatial density of these measurements is likely to impinge on the correlation lengh-scales of model error.Again, we need to understand how the ensemble variability projects into the columnintegral that will be measured by these instruments.
Considering tall tower based measurements with several levels, the vertical covariances we described here (see sec.3.2) show significant correlation between the first model level with the higher levels, up to 0.5 during the day and about 0.4 during the night.The use of several observation levels implies then model covariances that reduce the independence of the data, and so the constraint on the fluxes, but which we need to defined properly to avoid unrealistic corrections.Incorrect representation of the boundary layer height, for example, could lead to opposite increments between two levels.The same arguments hold for airborne profile data (Lauvaux et al., 2008).
For the limited observing network used here, temporal error correlations dominate spatial correlations.With the trend towards increasing numbers of continuous measurement sites this will remain the case for the surface network.It is already clear that models at coarse resolution perform poorly when simulating high-frequency observations over continents (e.g.Geels et al., 2007;Law et al., 2008).High resolution models (e.g.Lauvaux et al., 2008;Law et al., 2008) perform better. Biogeosciences, 0000, 0001-15, 2009 www.biogeosciences.net/bg/0000/0001/dynamics, especially in our case where the surface model is coupled online to the atmospheric model.The ground water content and vegetation description errors induce an important part of the nocturnal transport error.Two main difficulties remain to explore these error contributors.First of all, perturbing the internal parameters means exploring different transport models.Generating such an ensemble, while maintaining the physical realism of the simulation, is a difficult task indeed.Second, modelled CO 2 concentrations are affected by flux model errors.Variance is under-estimated in this context.However our chief aim in the ensemble modelling is the correlation structure and this seems more robust than the modelled variances.As an illustration, after 18 hours of simulation, the length-scales are about 50 km and remain constant during the next three afternoons.At the same time, the initial perturbations included in the boundary conditions are growing throughout the first 48 h.The identified length-scales are oriented differently depending on the meteorological situation but the norm of the error is constant in time.Additionally to these aspects, we defined the prior er-ror covariance matrix as uncorrelated and homogeneous over the domain, that is an over simplification of the truth, but would perturb the interpretation of the results if more realistic.Forthcoming inverse systems using CO 2 concentrations (not pseudo-data experiment) have to include the error flux correlations that could, for example, limit the impact of transport error correlations over the oceans (oceanic flux uncertainty is reduced compared to land fluxes).
The limited size of the ensemble could also act as a reduction of the true variance.The study of Talagrand et al. (1997) showed the linear decrease of the Brier Skill Score due to the finite size of the ensemble.Considering the correlation space, the number of singular vectors determines the estimation of the sensitivity.In our study, the 16 singular vectors selected over Western Europe are designed for short range forecasts over France.Even though we used a limited number of simulations, the operational ensemble forecast perturbs the critical regions controlling the synoptic conditions.The under estimation of the variance remains but the covariances are well-structured as shown in Fig. 5   Even here, however, it appears that in such models errors are highly correlated in time.
Finally the totality of these results has serious implications for quantitative network design such as that described by Kaminski and Rayner (2008).The model-data mismatch covariance matrix is a key input into such studies (see, e.g.Gloor et al., 2000).Our results suggest a minimum tion of surface measurements of about 50km before correlations reduce the marginal benefit of new observations.This is neither a function of model resolution nor of the footprint of measurements but of limitations in the simulation of transport.
This work hints at a promising partnership between studies in Numerical Weather Prediction and in atmospheric inversion.In general the characterisation of model error in NWP is much more advanced than in the field of inversions.It would be interesting to see how these errors appear at larger scales, both in space and time.We also need, as already noted, analogous developments for vertical and temporal correlations.

Conclusions
We have used an an ensemble prediction system coupled to a mesoscale transport model to estimate spatial and temporal correlations in the model-data mismatch for CO 2 inversions.Horizontal correlation lengths are of order 50km.There are strong vertical correlations in the boundary layer, particularly during the day.Temporal correlations are stronger than spatial and can last for most of a day.Taking account of these correlations reduces the effective information content of the mesoscale observations we use.The correlations also imply limits on the useful density of future observations.www.biogeosciences.net/bg/0000/0001/Biogeosciences, 0000, 0001-15, 2009 The error correlation length-scale used in the paper does not correspond to the classical definition used in trace gas inversions, i.e. the distance at which the correlation reaches 1/e.Here we use the differential length-scale (Daley, 1991, p110), related to the curvature of the correlation function at the origin.The estimation of this length-scale is simple in the isotropic case as it is related to the power spectrum of the auto-correlation (with the help of the Wiener-Kintchine theorem).Directional approximations of this length-scale, with emphasis on there sampling distributions, exist (Pannekoucke et al., 2008).The estimation is slightly biased with an over-estimation of the truth (Pannekoucke et al., 2008).The sampling distribution has an heavy tail with a positive skewness.It results that, in average, the sampling noise leads to an estimation larger than the truth (see e.g.Pannekoucke et al. (2008), Fig. 7).That is, the smooth of the raw lengthscales modeled with the diffusion operator appeared as beneficial revealing clearly the initial true correlation structures (Pannekoucke and Massart, 2008).
Two sub-ensembles were derived from the temporal correlations, corresponding roughly to nocturnal and diurnal periods (each of them limited to a few hours).The length-scales (L x and L y ) of the correlations are then associated with el-liptic functions on the horizontal plane.The ellipse is defined by its longer axis denoted L, the smaller axis denoted l, with the anisotropy 1 − l L , and the angle of the anisotropy θ.As the ensemble variations are flow dependent, preferential directions of anisotropy are expected over the domain.For example, diurnal circulation dominated by advection or local thermal gradients constrains the orientation of the largest correlations.
We have calculated the correlations using the whole period.This represents a trade-off.We need a correlation model which is stable and represents the internal dynamics of the model (suggesting longer periods of integration) but on the other hand the correlation structure really should depend on flow regime.Our approach averages over several regimes and hence averages several directions of anisotropy (Fig. 8a), and makes the magnitude of anisotropy spatially homogeneous (Fig. 8b).When comparing the length-scales L y (North-South) at the tower locations, the values are smaller (less than 30 km during the day), which corresponds to an East-West orientation of the anisotropy.In this region, the Atlantic ocean forces the error structure at Biscarosse, and the Garonne river valley does likewise at Marmande.
The above results have a range of impacts on current and   from those observations (Sect.3.5).This is not universal however.A correlation in the model-data mismatch implies a set of preferred directions for this mismatch.Mismatches in that direction will be deweighted (reducing available information) however those in orthogonal directions will be highlighted.This explains the mixed patterns in Fig. 8b.
The spatial correlation lengths we obtained do not suggest a large impact on most current observing networks since it is rare to have observations spaced less than 100 km apart.An immediate exception is measurements from aircraft.It appears that, even if the transport model we use for the inversion allows sampling at the resolution of a few km, spatial error correlation may reduce the effective sampling density.We cannot know the true import of this until we calculate the impact of horizontal correlations for various observational densities.
A less immediate but stronger consequence is the inversion of data retrieved from upcoming satellite missions (e.g.Crisp et al., 2004;Hamazaki et al., 2004).The spatial density of these measurements is likely to impinge on the correlation lengh-scales of model error.Again, we need to understand how the ensemble variability projects into the columnintegral that will be measured by these instruments.
Considering tall tower based measurements with several levels, the vertical covariances we described here (see Sect. 3.2) show significant correlation between the first model level with the higher levels, up to 0.5 during the day and about 0.4 during the night.The use of several observation levels implies then model covariances that reduce the independence of the data, and so the constraint on the fluxes, but which we need to defined properly to avoid unrealistic corrections.Incorrect representation of the boundary layer height, for example, could lead to opposite increments between two levels.The same arguments hold for airborne profile data (Lauvaux et al., 2008).
For the limited observing network used here, temporal error correlations dominate spatial correlations.With the trend towards increasing numbers of continuous measurement sites this will remain the case for the surface network.It is already clear that models at coarse resolution perform poorly when simulating high-frequency observations over continents (e.g.Geels et al., 2007;Law et al., 2008).High resolution models (e.g.Lauvaux et al., 2008;Law et al., 2008) perform better.Even here, however, it appears that in such models errors are highly correlated in time.
Finally the totality of these results has serious implications for quantitative network design such as that described by Kaminski and Rayner (2008).The model-data mismatch covariance matrix is a key input into such studies (see, e.g.Gloor et al., 2000).Our results suggest a minimum separation of surface measurements of about 50 km before correlations reduce the marginal benefit of new observations.This is neither a function of model resolution nor of the footprint of measurements but of limitations in the simulation of transport.
This work hints at a promising partnership between studies in Numerical Weather Prediction and in atmospheric inversion.In general the characterisation of model error in NWP is much more advanced than in the field of inversions.It would be interesting to see how these errors appear at larger scales, both in space and time.We also need, as already noted, analogous developments for vertical and temporal correlations.

Conclusions
We have used an an ensemble prediction system coupled to a mesoscale transport model to estimate spatial and temporal correlations in the model-data mismatch for CO 2 inversions.Horizontal correlation lengths are of order 50 km.There are strong vertical correlations in the boundary layer, particularly during the day.Temporal correlations are stronger than spatial and can last for most of a day.Taking account of these correlations reduces the effective information content of the mesoscale observations we use.The correlations also imply limits on the useful density of future observations.

Fig. 1 .
Fig. 1.Simulated CO2 (dotted lines) from the 11 simulations over the four days at Biscarosse (a) and Marmande (b) towers, compared to the observed concentration (solid lines)

Fig. 2 .
Fig. 2. CO 2 variability (standard deviation of the ensemble, in ppm) over the domain during the afternoon of the 26 of May 2005 (a) and during the night between the 26 to the 27 of May (b) (Marmande: diamond, Biscarosse: triangle) 5, and finally decrease to insignificant values in the upper troposphere.At night, the vertical profile of the standard deviation shows a different shape, decreasing rapidly Biogeosciences, 6, 1089-1102, 2009 www.biogeosciences.net/6/1089/2009/

Fig. 3 .
Fig. 3. Vertical CO 2 variances (left side) and associated vertical correlations (right side) over the domain during nighttime (solid line) and daytime (dashed line)

Fig. 4 .
Fig. 4. Temporal correlation (averaged over the whole domain and the first 3 days of the simulation) between midday to 9pm with the next 24 hours (a) and between midnight to 9am with the next 24 hours (b).Starting hours are indicated on the different lines, describing the temporal correlation averages of the first three days for each of the 24 different hours.

Fig. 4 .
Fig. 4. Temporal correlation (averaged over the whole domain and the first 3 days of the simulation) between midday to 9 p.m. with the next 24 h (a) and between midnight to 9 a.m. with the next 24 h (b).Starting hours are indicated on the different lines, describing the temporal correlation averages of the first three days for each of the 24 different hours.

BiogeosciencesFig. 5 .
Fig. 5. Spatial length-scales (in km) over the domain during the day and night, computed from the raw ensemble (a and c respectively)) and modelled using the diffusion equation (b and d respectively) .

Fig. 5 .
Fig. 5. Spatial length-scales (in km) the domain during the day and night, computed from the raw ensemble (a and c, respectively) and modelled using the diffusion equation (b and d, respectively).

Fig. 6 .
Fig. 6.Error reduction (in %) for the first 6-day inversion period (24 May to 30 May) using the uncorrelated covariance matrix (a) and the difference with the 'spatial only' (b), 'temporal only' (c), and combined (d) cases

Fig. 6 .
Fig. 6.Error reduction (in %) for the first 6-day inversion period (24 May to 30 May) using the uncorrelated covariance matrix (a) and the difference with the "spatial only" (b), "temporal only" (c), and combined (d) cases

Fig. 7 .
Fig. 7. Error reduction (in %) for the periods (30 May to 5 June, and 5 June to 10 June) using the uncorrelated covariance matrix (a) and (c), and their differences compared to the fully correlated matrix (b) and (d)

Fig. 8 .
Fig. 8. Angle (θ ) showing the main direction of the error (a) and the norm representing the ratio between the long and the short axis of the ellipse (1 − l L ) (b) of the anisotropy modelled with the diffusion operator.The two towers are the red circles on the maps.
after diffusion modelling.