Note : A new coupled system for global-to-regional downscaling of CO 2 concentration estimation

We introduce a global-to-regional nesting scheme for atmospheric transport models used in simulating concentrations of green house gases from globally distributed surface fluxes. The coupled system of the regional Stochastic Time-Inverted Lagrangian Transport (STILT) model and the global atmospheric transport model (TM3) is designed to resolve atmospheric trace gas concentrations at high temporal and spatial resolutions in a specified domain e.g. for regional inverse applications. The nesting technique used for the coupling is based on a decomposition of the atmospheric concentration signal into a far-field and a near-field contribution enabling the usage of different model types for global (Eulerian) and regional (Lagrangian) scales. For illustrating the performance of the coupled TM3-STILT system we compare simulated mixing ratios of carbon dioxide with available observations at 10 sites in Europe. For all chosen sites the TM3STILT provided higher correlations between the modelled and the measured time series than the TM3 global model. Autocorrelation analysis demonstrated that the TM3-STILT model captured temporal variability of measured tracer concentrations better than TM3 at most sites. Correspondence to: K. Trusilova (ktrusil@bgc-jena.mpg.de)


Introduction
Anthropogenic and biogenic emissions of trace gases at the surface cause large variations of atmospheric concentrations in the Planetary Boundary Layer (PBL), and upstream sources and sinks influence tracer mixing ratios at observational stations.However, the scales on which the transport in this "near-field" often takes place are not resolved by typical grid sizes of global atmospheric transport models.The comparison against available observations is therefore difficult and often requires averaging of the observed quantities in time and space to the model resolution.Such averaging leads to a loss of information about fine scale signals that can be retrieved from trace gas fluxes at high temporal and spatial resolutions by an atmospheric inverse scheme (Lin and Gerbig, 2005;Gerbig et al., 2008).
In order to refine the spatial and temporal resolutions of the predicted concentrations of trace gases in global transport models, domain nesting techniques are often used.Such a zooming technique was implemented in the atmospheric chemistry-transport zoom model TM5 (Krol et al., 2005) and in the nested regional transport-chemistry model HANK (Hess et al., 2000).In order to smooth the transition from the global grid of several degrees to a fine grid of less than 1 • ×1 • typically, often one or more intermediate zoom regions are needed.This multi-domain nesting leads to an increased computational cost of model runs and can be justified only when the whole zoomed area is of interest.However, for a model evaluation or an inverse model study, usually Published by Copernicus Publications on behalf of the European Geosciences Union.
only the model output at a limited number of available observation locations is required.Such a local refinement of the spatial resolution for predicted variables may be achieved by using Lagrangian transport models (Rotach et al., 1996;Lin et al., 2003;Stohl et al., 2005;Pongratz et al., 2008) for calculating the footprint (near-field) and background (far-field) influences for each available observation.Lagrangian transport models simulate transport by following the time evolution of a particle ensemble and interpolating meteorological fields to the sub-grid scale location of each particle.The particle approach can provide significant gains in information as compared to Eulerian-grid models due to the adjustable time step and fine spatial resolution of surface sources and sinks (Lin et al., 2003).
In this paper, we present a new fully coupled system of a global three-dimensional atmospheric tracer model and a stochastic Lagrangian transport model designed for predicting atmospheric trace gas concentrations at fine temporal and spatial resolutions.The coupled system is to be used in the two-step scheme for high-resolution regional atmospheric trace gas inversions described in the work of Rödenbeck et al. (2009).This scheme allows coupling of two independent atmospheric transport models for retrieving surface fluxes of observed trace gases at finer regional spatial resolutions than possible with global models.Rödenbeck et al. (2009) presented the scheme using a hypothetical regional model nested into the global one; here we demonstrate an implementation of the coupled system consisting of a global Eulerian grid model and a regional Lagrangian model to be used in this scheme.In this study we analyse the performance of the coupled system for the European domain and evaluate it against continuous measurements of CO 2 concentrations at 10 European sites.

Modelling system
The global-regional coupled system of the global Euleriangrid and the regional Lagrangian atmospheric transport models is highly modular and both components can be decoupled at the stage of compilation.The technical implementation of the coupling is consistent and it is based on the principle of decomposition of measured tracer mixing ratio on the "near-field" and "far-field" contributions.This decomposition is possible because of the linearity of the atmospheric continuity equation for a conservative tracer.The near-field part of the signal is resolved by the fine-scale regional model whereas the far-field part -the relatively well-mixed atmosphere background -is resolved by the global coarse-grid model.
The global atmospheric tracer model TM3 (Heimann and Körner, 2003) is a three-dimensional transport model that solves the continuity equation for an arbitrary number of atmospheric tracers.It uses prescribed meteorological fields, which can be obtained from an atmospheric general circulation model or from a weather forecast model.6 hourly meteorological fields from the NCEP/NCAR Global Reanalysis Products (Kalnay et al., 1996) were used for our simulations.
The Stochastic Time-Inverted Lagrangian Transport model STILT (Gerbig et al., 2003;Lin et al., 2003) simulates atmospheric transport by following the time evolution of a particle ensemble backwards in time and calculates footprint elements that represent the sensitivity of the mixing ratio of the tracer at the receptor location to any given surface flux.The model uses meteorological fields from the ECMWF1 operational data at the spatial resolution of 0.25 degree horizontally with 60 vertical levels on the latitude/longitude projection that have been obtained from the ECMWF Data Server.The particle number was set to 100 (per receptor point), providing the statistical error of ∼13% of the regional contribution to the CO 2 mixing ratio at the receptor point (Gerbig et al., 2003).As described by Lin et al. (2003), STILT computes the influence function I (x r ,t r |x,t), which links the local concentration C(x r ,t r ) of a tracer, measured at receptor location x r at time t r to the surface sources S(x,t) of the tracer (in ppm units) emitted at locations x at time t into a shallow layer near the surface: The first term on the right-hand side in Eq. (1), C(x r ,t r ) NF , represents the upstream near-field influence to the concentration at the receptor location.The function I (x r ,t r |x,t) links volumetric sources S(x,t) to the concentrations at the receptor C(x r ,t r ).The layer, which is considered to be in contact with the surface, is as thick as 50% of the local mixing height.The term C(x r ,t r ) NF is calculated along trajectories of particles while they "interact" with the surface layer and the corresponding air parcels exchange tracers with surface sources and sinks.A trajectory of a particle is the path followed by the particle along the mean wind direction starting at the location (receptor) at time t r .The trajectory length is limited by the maximum backward time of calculation, t 0 (x)≤72 (hours) for this study, and by the limits of the regional model domain.
The second term in Eq. (1), C(x r ,t r ) FF , represents the influence of the background concentration of the given tracer to its concentration at the receptor.This background far-field CO 2 concentration is calculated by the global model.We compare the performance of the global TM3 model versus the regional TM3-STILT coupled model within the DoI/PoI.The far-field influence C FF is the same for both model simulations, whereas the near-field (C NF ) contributions resolved by the global and the regional models differ.
The C NF value was sampled from both models in the following way: the sampling level in the TM3 model was defined above the sea level at which the mixing ratio of the corresponding grid cell was taken; the sampling height in the STILT model was defined above the ground level.Values of C NF from both models were compared to C NF extracted from observations.

Calculating the far-field part of the mixing ratio signal
The far-field part of the signal C FF is calculated by the global model in two steps as in the 2-step inversion scheme (Rödenbeck et al., 2009): 1. the global three-dimensional distribution of carbon dioxide mixing ratio in the atmosphere (C glob ) for GD/FP is calculated with the global TM3 model using an inversion-retrieved estimate of the carbon surface flux (F posterior ); 2. the near-field contribution in the global simulation is subtracted from the global CO 2 field C glob .
For retrieving F posterior we used an inversion similar to the global inversion described in the work of Rödenbeck (2005) but with the following prior estimates of global surface flux fields: the anthropogenic emission fields taken from the EDGARv3.2emission database (Olivier et al., 2005;van Aardenne et al., 2005) that was linearly extrapolated to 2006, the net biosphere-atmosphere exchange calculated with the terrestrial ecosystem model BIOME-BGC (Trusilova and Churkina, 2008), and the net ocean-atmosphere carbon flux (F o ) that was specified as the sum of the ocean uptake flux induced by the anthropogenic perturbation as compiled In order to estimate in the European region the far-field contribution C FF only, the near-field regional contribution C NF has to be subtracted from the C glob .The C NF is calculated with the global model within the DoI/PoI using the modified flux input F posterior 0 ={F posterior within DoI/PoI; 0 -otherwise}.Finally, the remaining C FF is calculated as: where A TM3 is linear transport operator of the global model; A TM3 is the transport operator of the global model, A TM3 , modified for the comparison with the regional model in the following way: any tracer emitted within the DoI/PoI and transported beyond the DoI/PoI cannot re-enter it; any tracer field outside DoI/PoI is set to 0 in order to avoid double counting of the same tracer mass contribution to the near field (within the DoI/PoI) and to the far field (outside the DoI/PoI).The near field signal C NF has to be subtracted from C FF in order to avoid double counting of the near field influence in the global atmosphere.
The far-field mixing ratio C FF is calculated using the optimized surface flux estimate F posterior retrieved by a global atmospheric inversion in a previous step.This guarantees that the global flux field, F posterior and, subsequently, the C FF are consistent with the observations globally.With this procedure, C FF is a good coarse scale estimate of the far-field influence given by the global atmospheric transport model.

Calculating the near-field part of the mixing ratio signal
The near-field contribution resolved by the global and the regional models -C NF TM3 and C NF STILT , respectively -are calculated within DoI/PoI with the set of fluxes F upscaled to spatial resolutions of 4 • ×5 • and 0.25 • ×0.25 • for the global model and the regional model: (5) where A STILT is the linear transport operator of the regional model.
For calculating C NF TM3 and C NF STILT , three components of the carbon surface flux were specified with high temporal and spatial resolution: the anthropogenic carbon emission flux (F a ), the land biosphere (F b ) and the ocean (F o ) net carbon exchange fluxes.The F a was derived from the European anthropogenic emission inventory at a resolution of 10 km (Pregger et al., 2007).The F b VPRM was calculated with the diagnostic model VPRM (Mahadevan et al., 2008) for the DoI/PoI at 1/12 • ×1/8 • .The VPRM model was driven with the ECMWF operational meteorological data and used the land cover classification database SYNMAP (Jung et al., 2006).The

Evaluation of the coupled TM3-STILT model
For the evaluation of the models two simulations were performed for the DoI/PoI domain: (1) on the coarse spatial/temporal resolution with the TM3 model and (2) on the fine spatial resolution with the TM3-STILT coupled model.
The time series of the near-field CO 2 mixing ratios, C NF TM3 and C NF STILT , were then compared with the nearfield signal (C NF obs ) extracted from the time series of observations (C obs ) that was estimated as:

Data selection
As the coupled TM3-STILT system is designed for performing regional atmospheric inversions on high spatial resolutions, for the evaluation we chose the sites that will potentially be used in a regional inversion for the European domain (Table 1).The CO 2 concentration data were provided by the measurement institutions listed in the Table 1.The data selection procedure for continuous and for flask data is described in the work of Rödenbeck (2005).
Although the measured and modelled data are continuous and of high frequency (1 h), we chose to analyse the model output only between 10:00 and 17:00 UTC in order not to include night time values because of several considerations: 1. the night-time mixing ratios near the surface, where the stations are located, are very much influenced by small scale meteorological circulation patterns which are mostly not resolved in either global or regional models; 2. the day-time values measured at tall towers or mountain stations are representative of the well-mixed planetary boundary layer and are better represented in atmospheric models.

Subtraction of the seasonal signal
The large-scale seasonal variability of the CO 2 concentrations over Europe is dominated by the surface flux from the entire northern hemisphere, and which is represented primarily in the far-field component of the mixing ratio.Here the synoptic variability is of interest and, therefore, prior to the statistical analysis we subtracted the seasonal component from the observed and modelled time series.Otherwise, including the season cycle in the time series would increase the model-to-data match (due to correlations governed by the large scale seasonal patterns) but such approach would not allow a fair analysis of the role the transport scales play in representing the fine scale variability of trace gas signals.
A harmonic function was used to approximate the seasonal flux component for each time series: where c(t) is time series of measured | modelled concentration values, c 0 ,c 1 ,a i ,b i are constants estimated by the harmonic-fitting algorithm.t is time expressed in days of the year.The seasonal signal approximated by the harmonic function c(t) was subtracted from the time series C NF TM3 , C NF STILT , and C NF obs .The resulting time series X TM3 , X STILT , and X obs were then statistically analysed.

Statistical measures for the model evaluation
For comparing the performance of atmospheric transport models we calculated the correlation coefficient (r) and the Root Mean Square Deviation (RMSD) between modelled and measured time series of CO 2 (only for the near-field influence parts X TM3 , X STILT , and X obs ), as well as the autocorrelation coefficient (R) and the variance (Var) for these time series.All gaps in the observed CO 2 mixing ratios were linearly interpolated in time.The autocorrelation analysis was applied for the same set of lag periods (1-8, 18, 24, 48, 72 and 96 h) to the interpolated time series of measured and modelled data in order to quantify the global and the coupled model's capability to capture the fine scale variability in the measured CO 2 signals.

Results and discussion
Figure 2 shows as an example the de-seasonalized time series of the observations together with the simulations by TM3 and TM3-STILT at a flat-terrain station, a mountain station, and a coastal site.The TM3 model produces a smoother time series than TM3-STILT and fails to capture many of fine-scale features of the observed signals, in particular at the Monte Cimone site (CMN).

Analysis of time series using Taylor diagrams
In order to evaluate the TM3 and TM3-STILT capabilities to capture fine-scale features of modelled signals we make use of Taylor diagrams.The Taylor diagram (Taylor, 2001) shows graphically how closely a pattern matches observations.The similarity of two patterns is qualified in terms of their correlation (angle), their standard deviation (or relative standard deviation -radius), and their root mean square deviation (also root mean square error -distance from the location of pattern to the location of observation on the diagram).
Figure 3 shows Taylor diagrams that demonstrate how well the TM3 and TM3-STILT models are able to simulate the CO 2 concentrations at given locations; the match to the observation is shown by the position of each dot.The statistical characteristics used for the Taylor diagrams of Fig. 3 are given in Table 2.
For all sites the TM3-STILT model showed a better correlation of the modelled CO 2 time series to observations than the TM3 model (Table 2).Especially large differences were found for the station at Mace Head (MHD), for which the r TM3,obs =0 and r STILT,obs =0.63.The reason for this large difference can be explained by the capability of the TM3-STILT model to resolve the fine-scale transport in this coastal region.
The variance of the TM3-STILT-simulated time series matches the variance of the observations closer than the TM3-simulated time series.Good performance of TM3-STILT was achieved for CMN, KAS, PUY, OXK, and SCH as illustrated by the orange dots located close to the Rel.SD=1 line in the Fig. 3.For the BIK, CBW, HUN, and MHD the Var STILT matches the variance in observations closer than Var TM3 , however, some of the variability was still missed by the fine-scale model.For the PRS station the regional model produced time series which strongly overestimated the observed variability.This could be explained by an inaccurate placement of the receptor point in the Lagrangian model over the surface with complex terrain: Plateau Rosa is a high-altitude mountain station with a very complex local meteorology that is not adequately represented in either TM3 or TM3-STILT.
The RMSD TM3−STILT was significantly smaller than RMSD TM3 for four stations: BIK, CBW, KAS, and MHD.These sites have rather homogeneous surrounding land cover and, except for KAS, are located on flat terrain.Under these conditions the surface flux uncertainly has a minor influence on the atmospheric signal and TM3-STILT with the fineresolved transport over-performes the TM3 model.On the other hand, at four sites -CMN, PRS, PUY, and SCH -RMSD TM3−STILT was greater than RMSD TM3 .This indicates a greater average mismatch between the TM3-STILT model and the observations.However, this mismatch originates primarily from the larger variability of the TM3-STILT signal predicted for these sites, however the corresponding time series of TM3-STILT are still better correlated to the observations than the estimates by the global model TM3.

Autocorrelation analysis
The autocorrelation analysis of the modelled and measured time series was performed in order to demonstrate the capabilities of the models to represent the statistical properties of the observed CO 2 signals.Because of the daytime selection procedure in the model output, the autocorrelation coefficients in the time series of model outputs and observations was calculated only for lag times of 1-8, 18, 24, 48, 72 and 96 h.Since observations were not always available at 1hourly intervals, linearly interpolated values at the respective time lags were used.The associated 95% confidence intervals for the estimated autocorrelation coefficients were calculated using the Fisher transformation (Shen and Lu, 2006).
The confidence intervals are shown as error bars of the autocorrelation coefficients in Fig. 4. The narrow confidence intervals illustrate that the autocorrelation coefficients calculated for the models and the observations become significantly different already from a lag of 1 h for all sites, except for the PUY and BIK sites.
Figure 4 shows the autocorrelation coefficients and their confidence intervals for the modelled and measured time series of CO 2 mixing ratio.For all time series the autocorrelation coefficients were high for short time lags of 1-3 h.As we considered only daytime observations, this strong dependence of neighbour values can be explained by small changes of CO 2 mixing ratio in 1-3 h during the day.For the larger time lags of 4 to 8 h the autocorrelation in observation was slightly decreased for modelled and measured time series.
For time lags greater than 8 h the TM3 time series still exhibits too high autocorrelation indicating too smoothed out day-to-day variations of tracer mixing ratios (Fig. 4, for BIK, CBW, CMN, MHD, OXK, PRS, PUY, SCH) in contrast to the observations.The autocorrelation coefficients of the TM3-STILT time series matches much more closely the measurements (Fig. 4) and they are decreasing for time lags greater than 8 hours at almost all sites (except HUN and KAS, Fig. 4), in much better agreement with the observations.This result demonstrates that the TM3-STILT model has a better capability to capture the variability of local fine-scale tracer time series.
The poor behaviour of the autocorrelation curve at the KAS station may be explained by errors in positioning the receptor point at this mountain site.The STILT model requires setting the receptors vertical coordinate as the height above ground, which can be easily measured or estimated.However, the receptor point is placed above the smoothed terrain in the model.Therefore, the model output at this sampling height may actually correspond to a different height above sea level (where the observations refer).The positioning of the sampling height remains the major problem for modelling of mountain regions.

Conclusions
This paper demonstrates a practical implementation and evaluation of the spatial-temporal nesting of a Domain of Interest (DoI) and a Period of Interest (PoI) resolved by a regional model into a global atmospheric transport model.The nesting scheme, previously described in the work of Rödenbeck et al. (2009), was designed for performing high-resolution regional atmospheric trace gas inversions.The fine-resolution regional atmospheric transport model of the Lagrangian type was nested into the Eulerian-grid coarse resolution global model for one year over the European domain.The performance of the nested system versus the global model analysed for predicting time series of CO 2 mixing ratios at ten measurement stations.
The nested high-resolution model system is able to resolve both atmospheric transport and the source/sink patterns of CO 2 on fine spatial (less than 1 • ) and temporal (1 h) scales and to capture to a considerable degree the fine scale variability in measured signals that is not resolved by the global model.
The illustrative evaluation of the performance of the coupled TM3-STILT model is of course limited.On the one hand, the far field contribution by its construction from a global inversion matches the atmospheric observations outside the DoI.On the other hand, the near-field component of the simulation is dominated by the sources within the DoI.These sources have not been optimized in any way and are expected to contain errors in their spatial and temporal distribution.These errors are added to the atmospheric model transport errors and contribute to the mismatch at the observation stations.The statistical analysis performed here attempts to isolate the former error, but this cannot be done unequivocally.
TM3-STILT constitutes a computationally efficient modelling system for the representation of the high-resolution atmospheric transport of conservative atmospheric tracers over a complex domain such as Europe.In the future, the implemented system is to be used in regional atmospheric inversions for the optimization of high-resolution source flux estimates of CO 2 and other long-lived greenhouse gases.
Two domains were defined for the model simulations: the Global Domain (GD) with 72×48 grid cells on the spatial resolution of 4 • ×5 • , and the regional domain of interest (DoI) covering central and western Europe from 12 • W to 35 • E and from 35 •N and 62 • N with 108×188 (latitude by longitude) grid cells corresponding to a spatial resolution of 0.25 • ×0.25 • of the tracer flux field.The layout of the two grids is shown in Fig. 1.The full period of simulations (FP) extends from 1 January 2005 to 1 January 2007 for the global model.The period of interest (PoI) extends from 1 March to 31 October 2006.The time between 1 January 2005 and 1 March 2006 is the spinup time used in the global model to acquire the initial background distribution of the trace gas mixing ratios at the beginning of PoI.

Fig. 1 .
Fig. 1.The global and the regional domains.
F a and F b were aggregated to the resolution 0.25 • ×0.25 • for the DoI/PoI.The F o , which was used in the global inversion discussed above, was re-sampled to 0.25 • ×0.25 • within the DoI/PoI.The corresponding input files contain carbon emission fluxes [mol m −2 3 h −1 ] on the latitude-longitude projection.

Fig. 2 .
Fig. 2. Time series of the measured and modelled near-field CO 2 signal for year 2006 at three locations: BIK -a flat-terrain continental tall tower site in eastern Poland, CMN -a mountain station in Italy, and MHD -coastal station in western Ireland.The seasonal component of the signal has been subtracted from the time series.

Fig. 3 .
Fig. 3. Taylor diagrams for modelled and measured time series of CO 2 mixing ratio at 10 European sites.The Relative Standard Deviation (Rel.SD) was calculated as the ratio of the SD of modelled time series to the SD of observations (SD mod /SD obs ).The black point indicates the perfect match between the model and the observations.

Table 2 .
Statistical characteristics of modelled time series (global mean and seasonal cycle subtracted): r -correlation coefficient, RMSD -Root Mean Square Deviation between observed and modelled time series, and Var -variance of the measured and modelled tile series.