Comparing Lagrangian and Eulerian models for CO 2 transport – a step towards Bayesian inverse modeling using WRF / STILT-VPRM

We present simulations of atmospheric CO 2 concentrations provided by two modeling systems, run at high spatial resolution: the Eulerian-based Weather Research Forecasting (WRF) model and the Lagrangian-based Stochastic Time-Inverted Lagrangian Transport (STILT) model, both of which are coupled to a diagnostic biospheric model, the Vegetation Photosynthesis and Respiration Model (VPRM). The consistency of the simulations is assessed with special attention paid to the details of horizontal as well as vertical transport and mixing of CO 2 concentrations in the atmosphere. The dependence of model mismatch (Eulerian vs. Lagrangian) on models’ spatial resolution is further investigated. A case study using airborne measurements during which two models showed large deviations from each other is analyzed in detail as an extreme case. Using aircraft observations and pulse release simulations, we identified differences in the representation of details in the interaction between turbulent mixing and advection through wind shear as the main cause of discrepancies between WRF and STILT transport at a spatial resolution such as 2 and 6 km. Based on observations and inter-model comparisons of atmospheric CO 2 concentrations, we show that a refinement of the parameterization of turbulent velocity variance and Lagrangian time-scale in STILT is needed to achieve a better match between the Eulerian and the Lagrangian transport at such a high spatial resolution (e.g. 2 and 6 km). Nevertheless, the inter-model differences in simulated CO 2 time series for a tall tower observatory at Ochsenkopf in Germany are about a factor of two smaller than the model-data mismatch and about a factor of three smaller than the mismatch between the current global model simulations and the data.


Introduction
Inverse modeling tools use the atmosphere as an "integrator" to obtain information on the source-sink distribution of CO 2 on different spatial and temporal scales.A common practice is to use a global atmospheric transport model together with a network of atmospheric measurements to estimate the relationship between flux and tracer distributions via inverse techniques.The reliability of the inverse flux estimation depends largely on the quality of the transport represented in the models (Gerbig et al., 2008;Lin and Gerbig, 2005;Stephens et al., 2007;Geels et al., 2007).
Atmospheric transport models can be based on either Eulerian or Lagrangian formulations of the fluid transport process.In the Lagrangian formulation, the motion of fluid elements is described by solving the Lagrangian equations of mass and momentum along the trajectory of the particle/fluid element.In the Eulerian approach, the mass concentration of fluid elements is calculated as a function of space and time instead of calculating trajectories of fluid elements.Despite this difference, these two approaches solve numerically the same partial differential equation; hence theoretically they lead to identical results when the temporal and spatial resolutions are sufficiently increased and if the same parameterization for sub-grid scale transport is applied in each of them (Hanna, 1979;Lee and Stone, 1983).Both approaches are used in the inverse modeling community to estimate source-sink distributions (Gerbig et al., 2003a;Lauvaux, 2008;Rödenbeck et al., 2003) The atmospheric distribution of passive inert trace gases is variable on small scales (both temporal and spatial), caused by strong surface flux variability in the near field and by mesoscale transport phenomena.However, the current global models, with spatial resolutions of no more than 1 • × 1 • , fail to resolve these variations on measured atmospheric CO 2 , which potentially leads to biases in flux estimates (Ahmadov et al., 2009).In order to better represent measurements made in the mixed layer (the lowest 1-2 km of the atmosphere) by stations such as tall towers, the inverse system requires the atmospheric transport models to be set up at high spatial resolution (2-20 km).In addition, the fluxes in the near-field of the observatories can be highly variable (Gerbig et al., 2003b), calling for a-priori fluxes to be specified at high spatial resolution.Recent studies have demonstrated improvement in capturing the variability of observed CO 2 concentrations when increasing the spatial resolution of the transport models (Ahmadov et al., 2007;Pérez-Landa et al., 2007;Sarrat et al., 2007;van der Molen and Dolman, 2007).Gerbig et al. (2003a) describes a receptor-oriented framework using a Lagrangian Particle Dispersion Model (LPDM) together with lateral boundary conditions and a biospheric flux model to derive regional fluxes at high spatial and temporal resolution.The "footprints" (sensitivity of model output (e.g.concentration) to input variables (e.g.surface fluxes)) derived from an LPDM is similar to the adjoint of an Eulerian Transport model (Errico, 1997).Using these footprints has the advantage of resolving the fine structures originating from surface flux variations on scales smaller than the grid size of the meteorological fields used.In the case of the Eulerian approach, the models are affected by numerical diffusion, limiting the resolution to scales larger than the grid size in the underlying meteorology.The framework is thus analogous to a regional adjoint model in a Eulerian framework, providing an alternative to generating and implementing adjoint model code for a Eulerian transport model.
We use a framework similar to that introduced by Gerbig et al. (2003a) which consists of a receptor-oriented transport model driven "offline" by assimilated meteorological fields, an Eulerian "online" transport model and a diagnostic biospheric model to derive regional flux estimates.The receptororiented transport model is the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Gerbig et al., 2003b), the Eulerian transport model is the Weather Research Forecasting (WRF) model (http://www.wrf-model.org/) and the biosphere model is the Vegetation Photosynthesis and Respiration Model (VPRM) (Mahadevan et al., 2008).The term "online" indicates here that the meteorological fields are simulated during the model run, while "offline" refers to the use of previously computed meteorological fields which are read in during the simulation.
The wind fields generated by WRF are used in STILT to calculate ensembles of back trajectories starting at a receptor location (Nehrkorn et al., 2010).Resulting footprints (sensitivities to upstream surface-atmosphere fluxes) are then mapped to high-resolution biospheric fluxes that are prescribed from a diagnostic biosphere model, anthropogenic fluxes, as well as initial/lateral boundary conditions from a global model.This part of the framework -offline Lagrangian modeling -provides time series of CO 2 mixing ratios at the receptor location.The other part of the framework -the online Eulerian modeling -generates 3-D fields of CO 2 concentration, using the same surface fluxes and boundary conditions as the Lagrangian system.Hence the framework allows for a direct comparison of Eulerian (forward) and Lagrangian (adjoint) models to assess the consistency in simulating transport, which is a prerequisite for using STILT for the inverse estimation of fluxes.A schematic representation of the modeling framework is illustrated in Fig. 1.
This paper presents the simulated CO 2 time series generated by the Eulerian and Lagrangian transport models at high resolution for a domain centred over the atmospheric monitoring station Ochsenkopf, located in the Fichtelgebirge mountain chain in northern Bavaria, Germany (50 • 1 48 N, 11 • 48 30 E).The consistency of those two simulations is assessed with special attention paid to the details of horizontal and vertical transport and mixing.The consistency between model parameters and its impact on tracer simulations is examined in detail.The paper is organized as follows: Sect. 2 describes the major components of the modeling framework and the model domain.Results are presented and discussed in Sect.3, exploring reasons for possible discrepancies between modeled mixing ratios from two modeling systems.Section 4 provides the conclusion of this study.

Modeling framework
Here we describe major components of the Eulerian and the Lagrangian parts of the modeling framework: the coupled models WRF-VPRM (Eulerian), which provides spatial and temporal distributions of CO 2 , and WRF/STILT-VPRM (Lagrangian), which simulates the temporal distribution of CO 2 at the observation point (receptor).Including the Eulerian and Lagrangian models in a single framework allows for the quantitative comparison between the two different approaches, while using the same domain, surface fluxes and initial/ lateral boundary conditions.

WRF-VPRM model
As indicated by the name, this coupled model consists of the mesoscale weather prediction model WRF with a passive tracer option from WRF-Chem (Grell et al., 2005) together with the biospheric model VPRM to simulate the distribution of CO 2 .WRF-VPRM has been used in many regional applications and shown remarkable skill in capturing fine-scale spatial variability of CO 2 mixing ratios (Ahmadov et al., 2007;Ahmadov et al., 2009;Pillai et al., 2010).The model implementation is described in detail by Ahmadov et al. (2007).A brief description of the model is given as follows.We use WRF/Chem version 3.0 (hereafter referred to as WRF) with a tagged tracer option to distinguish different components (i.e.biospheric, anthropogenic etc.) of CO 2 .A K-diffusion scheme with heat exchange coefficient -K his used in WRF to account for turbulent vertical mixing of tracers (Grell et al., 2005).Note that the vertical diffusion of meteorological parameters is performed by the boundary layer scheme in WRF.Also note that convective fluxes are not used for tracer transport, but used cumulus parameterization scheme for meteorological parameters.An overview of the WRF physics/dynamics options used here is given in Table 1.
VPRM is a satellite-data based diagnostic biosphere model which uses MODIS (http://modis.gsfc.nasa.gov/)satellite indices as well as observed or simulated meteorological variables to calculate Net Ecosystem Exchange (NEE) at high spatial resolution (Mahadevan et al., 2008).In our model setup, VPRM computes biospheric fluxes utilizing the meteorological variables from WRF and then passes these fluxes to WRF to be transported in the passive tracer mode.SYN-MAP data (Jung et al., 2006) with a spatial resolution of 1 km and 8 different vegetation classes are used to specify the different types of vegetation cover in the domain.The VPRM parameters are optimized against eddy flux measurements for different biomes in Europe collected during the Integrated EU project "CarboEurope-IP" (http://www.bgc-jena.mpg.de/bgc-processes/ceip/).
Fossil fuel emission data at a spatial resolution of 10 km are prescribed from an inventory provided by IER (Institut für Energiewirtschaft und Rationelle Energieanwendung), University of Stuttgart (data available at http://carboeurope.ier.uni-stuttgart.de/)to account for anthropogenic fluxes.Both biospheric and anthropogenic surface flux inputs are projected to the Lambert Conical Cartesian co-ordinate system used by WRF-VPRM.Projection of gridded fossil fuel emissions to the WRF grid is done using mass  (Rödenbeck, 2005), generated by the global atmospheric tracer transport model, TM3 (Heimann and Koerner, 2003), based on optimized fluxes transported at a spatial resolution of 4 • × 5 • , and a temporal resolution of 3 h (ana96 v3.0, http://www.bgc-jena.mpg.de/∼ christian.roedenbeck/download-CO2-3D/).Analyzed meteorological fields from the ECMWF model (http://www.ecmwf.int/),at a temporal and horizontal resolution of 6 h and approximately 25 km respectively, serve as initial and lateral meteorological boundary conditions for the WRF-VPRM.
We use the nesting option with a horizontal resolution of 6 km (parent) and 2 km (nested) as well as 41 vertical levels (lowest layer at about 18 m).Each simulation day starts at 18:00 UTC of the previous day, and continues with hourly output for 30 h.The first 6 h are used as meteorological spinup time.The initial conditions of the tracer concentrations are prescribed from the previous day of the simulation except for the first day of simulation where TM3 fields are used as mentioned above.The lateral boundary conditions are specified from TM3 fields.

WRF/STILT-VPRM model
STILT is a Lagrangian Particle Dispersion Model, which simulates ensembles of particles representing air parcels of equal mass, transported backward in time from an observation point (receptor) by mean winds and sub-grid turbulent winds (Gerbig et al., 2003b).The model has been used extensively in regional simulations and inversion studies for different greenhouse gases (Lin et al., 2003(Lin et al., , 2004;;Gerbig et al., 2003b;Miller et al., 2008;Gourdji et al., 2010;Göckede et al., 2010).A brief description of STILT is given as follows.

D. Pillai et al.: Comparing Lagrangian and Eulerian models for CO 2 transport
We used the STILT repository version 608, checked out on 3 July 2009.The turbulent flow is modeled as a Markov chain, where particles are transported at each time step using following equation: where u is the turbulent component of the mean velocity vector u, u is a random vector drawn from a normal distribution with a width equal to the variance of the vertical velocity (σ w ), t is the time step, and R is an autocorrelation coefficient which determines the standard random walk for the turbulent velocity components for each time step.R is expressed as: where T L is the Lagrangian time-scale in the horizontal (u) or vertical direction (w) that determines the degree to which particles keep the memory of previous motion.T L is set to zero for a random walk and large T L represents the advection by mean wind.Profiles for T L and σ w are derived from meteorological fields (Gerbig et al., 2003b).Here STILT footprints are driven by meteorological fields from the high-resolution mesoscale model, WRF (hereafter referred as "WRF/STILT" to indicate that STILT is driven by WRF meteorology) (Nehrkorn et al., 2010).The WRF-VPRM source code is modified to output the meteorological variables required to drive STILT.WRF/STILT computes changes in the tracer concentration C (x r , t r ) at the receptor location x r measured at time t r as the sum of changes in the tracer concentration at the receptor due to surface fluxes F in the domain V between initialization time t 0 and t r (denoted as "C surface (x r , t r )) and the contribution from the initial tracer field C(x, t 0 ) (denoted as "C bg (x r , t r )) (Gerbig et al., 2003b).C (x r , t r ) is expressed as: Here I (x r , t r |x, t) is the influence function which links spatially and temporally resolved surface source or sinks S (x, t) to the tracer concentration at the receptor and is expressed as: for a given number of particles (N tot ) released from the receptor and particle density ρ(x r , t r |x, t) at location x and time t.
The tracer concentration at the receptor due to fluxes F , denoted as C surface (x r , t r ), is expressed as: where f (x r , t r |xj, yk, ti) is given by Here h represents the column height into which the tracer is diluted (half of the planetary boundary layer (PBL) height in the current application), ρ is the column averaged air density and mair is the molar mass of air.The interested reader is refereed to Gerbig et al. (2003a) and Trusilova et al(2010) for more details.The term f (x r , t r |xj, yk, ti) links surface fluxes to concentration changes at the receptor and is denoted as the "footprint".The footprint derived here is analogous to a numerical version of the adjoint for WRF transport (Lin et al., 2003;Gerbig et al., 2003a;Errico, 1997).
We release 100 particles from a receptor point and WRF/STILT transports particles backward in time for a maximum of 3 days or until particles leave the domain.A maximum of 3 days is sufficient for all particles to leave the outer domain.WRF/STILT is used with a nested option where the wind fields are provided at the spatial resolution of the WRF inner domain (2 km × 2 km) until the particles leave the inner domain and afterwards at the spatial resolution of the parent domain (6 km × 6 km).The footprints are calculated according to Eq. ( 5) and are gridded to a maximum resolution of 2 km × 2 km.The horizontal size of the grid cells for resolving the footprint is dynamically adjusted according to the increase in footprint area in order to save computation time, as well as to avoid under-sampling of surface fluxes when the statistical probability of finding a particle in particular grid box becomes smaller (Gerbig et al., 2003b).
The surface fluxes including the VPRM biospheric fluxes, simulated at a spatial resolution of 2 km × 2 km, and the IER (anthropogenic) fluxes, interpolated to 2 km × 2 km are coupled to the transport according to Eq. ( 4) in order to estimate the associated surface flux contributions to the concentration field at the receptor (C surface (x r , t r )).The total CO 2 concentration at the receptor (C (x r , t r )) is calculated by adding the global background tracer distribution -C bg (x r , t r ) -to C surface (x r , t r ) as given in Eq. (3a), where the lateral tracer boundary conditions are prescribed from the TM3 global model.Note that we used the same surface fluxes and initial/lateral boundary conditions as given in Sect.2.1.The asterisk symbol denotes the position of a flat region which is described in Section 3.   and these variations are captured by models in most of the cases.Both models produce similar results in predicting CO 2 concentrations as indicated by the squared correlation coefficient, R 2 = 0.65.A summary of the statistics calculated from the model simulations for different tower levels is given in Table 2.As evident from the summary statistics, the models also produced similar results for other model levels in the boundary layer.Inter-model differences, e.g. for the 163 m level with a standard deviation of about 1.8 ppm, are about a factor of two smaller than the standard deviation of modelobservation differences.A similar result is also obtained at Ochsenkopf for other seasons (not shown).It is possible that the standard deviation, used as a measure of model performance, is dependent on time scale of analysis if there is difference in smoothness between model predictions relative to each other and relative to observations.that case, the standard deviation cannot necessarily represent the similarity between two model predictions and between model predictions and observations.In order to examine this, we carried out the model comparison at different time scales and the results are shown in Table 3.The consistency among the results at different time-scales confirms that there is no significant dependence of time scale of analysis on model performance relative to each other, and relative to the observations, which validates the use of standard deviation as the measure of models' discrepancy here.The evaluation of these models against observations indicates the model bias (calculated as mean of the difference over the whole period) in opposite signs, which explains the larger inter-model bias compared to the model (model-measurement) bias (see Table 2).
Note that these values cannot be used as a measure of discrepancies between models or between model and observations since it can be averaged out for longer time series.Not using convective fluxes for tracer transport in both models can likely be the reason for the large model-observation differences as compared to the inter-model differences, if there is an impact on observed CO 2 from convective transport by clouds.However, the comparison of model performance for non-convective periods (time series excluding the data where convective rainfall is greater than 0.5 mm) has showed no reduction in the standard deviation of model-observation differences (not shown) and this confirms that there is no impact of deep convection on these mismatches.
The discrepancies between the simulations, albeit smaller than the model-observation differences, prompt further investigation, especially since both models are driven with same meteorological and surface flux fields.The inter-model comparison is also carried out for another location (49 • 48 N, 11 • 45 E) with relatively flat terrain (451 m a.s.l.) within the nested model domain in order to examine whether the complex topography over Ochsenkopf plays a role on causing part of these mismatches (not shown).Note that no CO 2 concentration measurement is available at this location.A similar performance obtained for the flat region, with no better inter-model agreement, suggests that these discrepancies are not related to the complex orography over Ochsenkopf.Possible factors that can cause these discrepancies are (1) differences in turbulent transport within the boundary layer, (2) potential violation of mass conservation in the driving meteorology due to discrepancies in coordinate transformations during data processing procedures, (3) sensitivity to the number of particles used in the Lagrangian model, (4) differences in input flux fields and (5) different numerical parametrization of advection or convection in the models.Time reversibility of STILT (Gerbig et al., 2003b) and mass conservation in STILT when using WRF wind fields (Nehrkorn et al., 2010) have been confirmed for this setup, ruling out the lack of mass conservation as a possible reason.The results show no sensitivity to the number of particles used in STILT, giving rise to negligible bias (0.02 to 0.04 ppm) between STILT simulations with 1000 particles instead of 100.This confirms that the discrepancy is not caused by the choice of number of particles in the STILT.In the following sections, we examine which of the remaining factors contribute to the differences between the model simulations.

Mixing height Parameterization
Differences in the way turbulent transport within the mixed layer is represented in the models can lead to differences in the vertical distribution of surface flux influences and thus to differences in tracer mixing ratios.Also differences in the mixing height (z i ) have been shown to cause differences in CO 2 mixing ratios of several ppm during the growing season (Gerbig et al., 2008).Hence it is appropriate to examine the consistency of vertical mixing and associated turbulence parameterized in the models.In the current set-up, WRF derives z i based on the bulk Richardson number method.WRF/STILT also uses the bulk Richardson number method locally to calculate z i , utilizing profiles of atmospheric variables (temperature and wind) and their gradients provided by WRF.
The bulk Richardson number Ri b in both models is calculated as follows: where g is the gravitational acceleration, z is the height above ground and θ v is the virtual potential temperature.u and v refer to lateral wind components.The subscripts s and k refer to the lowest and k th model levels.The mixing height is defined as the first level at which Ri b becomes greater than the critical Richardson number Ri c (set to be 0.25).
The comparison of mixing height derived within WRF with those from WRF/STILT at the Ochsenkopf tower site for August 2006 (Fig. 4a) reveals that z i derived within WRF is found to be lower than that of WRF/STILT in certain periods of the nocturnal boundary conditions.Mixing heights from the two models are not in perfect agreement, with a squared correlation coefficient of R 2 = 0.65.This discrepancy was found to occur prevalently during cloudy conditions, with mixing height fields simulated by WRF exhibiting a strong spatial variability for periods with broken cloud cover, so that slight differences in horizontal interpolation of the meteorological fields result in large differences in diagnosed mixing heights.Indeed, removal of cloudy periods improved the inter-model agreement significantly (R 2 = 0.91).However an average z i discrepancy of about 35 % (when using all data) cannot be neglected from causing corresponding deviations in tracer mixing ratios.The sensitivity of simulated CO 2 concentrations to the inter-model difference in the parameterization of mixing heights is tested by using WRF derived z i in WRF/STILT, and the results are compared with the standard WRF/STILT set-up.Surprisingly, the comparison between standard and modified z i set-up in WRF/STILT reveals only slight differences (see orange dotted and black lines in Fig. 3).The probable reason for the smallness of the difference between these modeled tracer concentrations is the existence of simulated patchy mixing height fields as generated by WRF at a spatial resolution of 2 km 2 km.In the case of patchy spatial patterns of z i , an air parcel which was once in the boundary layer at one time step can be in the free troposphere (FT) at the next time step.Hence the mixing height, which usually acts as a barrier for vertical mixing, cannot act as such a barrier for such spatially variable mixing height fields, when advection over small distances can turn mixed layer air into FT air and vice versa.Differences in profiles of the variance of turbulent vertical velocities between standard and test runs (WRF-derived z i ) in WRF/STILT (not shown) are negligible, which also indicates that the local z i differences cannot significantly affect the turbulent mixing of tracers.In summary, this confirms that the differences in simulated CO 2 concentrations between WRF and WRF/STILT are not caused by differences in mixing heights.The summary statistics of this test run are also included in Table 2.

Biospheric fluxes, VPRM
Discrepancies in the biospheric fluxes between the modeling systems can cause differences in simulated CO 2 .The biospheric fluxes (GEE and Respiration) at the receptor location, derived from both modeling systems, are diagnosed for the proper treatment of the given meteorological fields (temperature and radiation) and the VPRM input parameters.Fig. 4b shows the simulated biospheric fluxes at the tower site for the period of August 2006 and suggests that the fluxes are consistent between the modeling systems (GEE: R 2 = 0.9, Respiration: R 2 = 0.98).A 10 % (2 %) deviation in simulated GEE (Respiration) between models is caused by the temporal interpolation of radiation and temperature fields within STILT.The possible differences in CO 2 concentrations caused by the flux discrepancy of about 10 % are estimated to be only 0.1 ppm (bias).This estimation was based on simulations of CO 2 concentrations generated by STILT with 10 % enhancement in GEE.Hence it is indicated that the model differences cannot be attributed to biospheric flux discrepancies.

Advection scheme: WRF and WRF/STILT
Another factor which can induce inter-model discrepancy is related to the differences in the details of vertical mixing and advection (shear) of both models, as their combination is responsible for horizontal spread in simulated plumes.The analysis of the vertical structure of tracer transport can give more insight.For this purpose we utilized the observations of CO 2 vertical profiles obtained during the Ochsenkopf aircraft campaign, which can provide a qualitative assessment on the model simulations.Both models, in general, are able to capture the vertical distribution of CO 2 variability relatively well and showed similar performance (Pillai et al., 2011).Here we have chosen the 20 October 2008 as an extreme case where WRF-VPRM and WRF/STILT-VPRM showed larger deviations.An elevated concentration of CO 2 was found in the valley south of Ochsenkopf (hereafter referred simply as the valley) during this period at around 10:00 UTC (i.e.before the full development of the convective mixed layer) (see Fig. 5).WRF-VPRM predicted a large contribution from fossil fuel fluxes (determined by using the tagged tracer CO 2.foss ) and simulated higher total CO 2 concentration in the valley, consistent with observations, while WRF/STILT-VPRM failed to capture this large accumulation of CO 2 in the valley.However, WRF/STILT-VPRM reproduced the CO 2 accumulation when the contribution from advected fossil fuel emissions (CO 2.foss ) is replaced with that given by WRF-VPRM (not shown).This result reveals that the "missing" accumulation of CO 2 in WRF/STILT-VPRM during this particular period is due to the failure in capturing the advection of the fossil fuel contribution in WRF/STILT-VPRM.
It is perhaps a surprising result when considering other periods in which WRF-VPRM and WRF/STILT-VPRM showed similar results on simulating CO 2 concentrations.Hence it is appropriate to investigate further the causes of this large discrepancy in simulating advection of tracers.The following section explores this by conducting different model sensitivity tests in WRF/STILT-VPRM, assuming WRF-VPRM gives fairly good predictions in this specific period based on its more reasonable performance in the above case study.

Pulse release experiment in WRF/STILT and WRF
A more comprehensive comparison of advection of tracers in both models can be studied by following the simulated transport of a plume emitted at a given location.This can give a vivid picture on any possible deviation of advection between models.Hence we attempted to release an emission pulse from a location where a strong potential influence of surface fluxes (as determined by STILT footprints for the above mentioned extreme period) and a relatively strong fossil fuel emission source is found, in order to quantify the effect of this emission on downstream CO 2 concentrations as simulated by both models.In this way, one can assess the  The emission source is defined in such a way that it emits a "pulse" with a total concentration of S conc at a particular time t.To simulate the pulse in STILT, N tot particles were released from the emission point (48.5 • N, 11.0 • E, release point 8 m a.g.l.) at 04:00 UTC and transported forward in time.Note that the likely source location (the above spatial co-ordinate) was found by looking at the STILT footprint at 04:00 UTC on 20 October 2008, predicted using backward when particles were transported backward from a receptor point (valley in Fig. 5a) at 10:00 UTC.As STILT was not able to capture the source region, the footprint was interpreted loosely with the aid of emission map and largest nearby fossil fuel emission sources were assessed.The resulting tracer concentration at a specified location down-stream R conc STILT is given by: where N R is the particle density at the receptor after taking into account air density differences between source and receptor locations.The receptor boxes were defined along the WRF/STILT particle trajectory locations at a given time with a horizontal dimension of 6 km × 6 km.The vertical dimension of the receptor boxes was roughly equal to the thickness of each WRF vertical layers and was placed at the respective WRF vertical level.In this way, one can reproduce STILT plume distributions with a grid cell size of 6 km × 6 km and with vertical levels corresponding to WRF.
In WRF the pulse emission was implemented as a tagged tracer flux field with a spatial resolution of 6 km × 6 km   (corresponding to the spatial resolution of plume simulations generated by STILT) and with a single non-zero value entry at the prescribed source pixel for time t = 04:00 UTC.S conc in Eq. ( 6) was given by the corresponding tagged tracer concentration at the first model level (∼8 m above ground) of the source pixel in WRF.
A comparison of the WRF/STILT and the WRF simulated plumes at 10:00 UTC, when the enhanced CO 2 was measured near Ochsenkopf, is shown in Fig. 6a and d for an atmospheric column from surface to ∼190 m (pressureweighted column average of lowest six model levels).The WRF simulated plume reached the aircraft location at 10:00 UTC with its northern edge, while the STILT simulated plume just misses it.The models show considerable differences in shape and advection of the plume under these nocturnal conditions.A relatively larger horizontal (east-west) spread of plume is simulated in WRF/STILT when compared to that in WRF.Also notably, STILT transported the plume much faster from the source point, without leaving any presence of plume close to the source location.
The above result provides a clear indication that the interaction of wind shear and turbulent diffusion is simulated differently in WRF and WRF/STILT.Note that the turbu-lent transport is realized using K-diffusion in WRF, while a stochastic process (Markov chain) is used in WRF/STILT (see Sects.2.1 and 2.2).In STILT, the spread of the plume is largely controlled by the rate at which the plume is turbulently mixed into the residual layer above the mixed layer, where wind speed and direction are different.In fact, limiting σ w (i.e.vertical turbulent velocity variance) to 1 cm s −1 in STILT reduces the east-west-extent of the plume dramatically (Fig. 6b).A reason for the unusually large values of σ w in the FT of up to 100 cm s −1 might be the coupling between the high-resolution meteorological fields (6 km horizontal, and 10 vertical levels below 2 km) used to drive the particles with the parameterization for σ w developed for coarser resolution fields.Similarly, the Lagrangian decorrelation time scale T L has some control on the residence time at low levels, where winds are slower, after release of the plume.It thus controls how strongly the plume is flushed away by advection.Indeed, reducing T L by a factor of ten causes the plume intensity close to the source location to increase significantly (Fig. 6c), and tends to result in a plume distribution that closer matches the one given by WRF.The intensity of FT turbulence determines the dissipation rate and the dilution of plume in the boundary layer.T L determines the turbulent  mixing between different vertical levels, i.e. larger T L causes the plume to be more efficiently transported by turbulent motion away from surface to altitudes with stronger mean-wind and thus stronger horizontal advection.
The above results of sensitivity tests reveal that one can expect differences in WRF/STILT-VPRM and WRF-VPRM simulations of CO 2 , corresponding to different tracer advection even though the same meteorological fields, surface fluxes, and vertical mixing are used.This explains the inter-model mismatches seen in the tower data analysis described in Sect.3. The discrepancy, albeit small, is largely caused by the model-to-model difference in advecting both biospheric and anthropogenic signals (not shown).The intermodel differences can be particularly high when signatures from strong sources such as fossil fuel emissions are transported from larger distances and dominate over those from local sources (e.g.biospheric fluxes).To achieve a better match between the Lagrangian and Eulerian transport models, a refinement of the parameterizations determining the profiles of σ w and T L in STILT is required.

Sensitivity to model resolution
Ideally both model dynamics should converge when spatial and temporal resolutions are further increased.Sensitivity of the inter-model mismatch to models' spatial resolution is examined by performing the model simulations at different horizontal resolutions such as 6, 12, 24 and 48 km. Figure 7 shows the dependence of inter-model or modelmeasurement mismatch on models' horizontal resolution for Ochsenkopf tower at a measurement level of 90 m during August 2006.The inter-model or model-measurement mismatch shows sensitivity to the horizontal resolution and the model-to-model or model-to-measurement agreement deteriorates with decreasing horizontal resolution of the models.Similar result is also found for other tower heights (not shown).Hence it is expected that the inter-model differences become smaller at higher resolution ( 2 km) and in this case the remaining differences between the models can be interpreted as due to the different parameterizations.However, it can be also speculated that at a very high resolution, the models are less sensitive to the difference in parameterizations since the important small scale (local) meteorological features are better resolved and not averaged out by the model grid.In this case the inter-model mismatch will only depend on the difference in numerical advection algorithms between Lagrangian and Eulerian models.

Summary and conclusions
In this paper, we have presented high-resolution simulations of CO 2 from online Eulerian (WRF-VPRM) and offline Lagrangian (WRF/STILT-VPRM) modeling systems for a domain over Ochsenkopf, Germany.The consistency between Eulerian and Lagrangian transport models in parameterizing turbulent mixing and in transporting CO 2 as a tracer is assessed using identical meteorological fields and surface fluxes.
Overall, the models show similar performance in predicting CO 2 concentrations at Ochsenkopf with high inter-model correlations.The factors which cause the discrepancies between the models have been investigated by comparing different model parameters.The inter-model difference in local z i is found to have a negligible impact on simulated CO 2 concentrations between models due to the presence of a leaky The consistency of advection schemes in WRF and WRF/STILT is examined by simulating CO 2 concentrations along an aircraft trajectory in the Ochsenkopf aircraft campaign.Both models provided similar results for most of the cases as described in Pillai et al., 2011; however we found a short period when WRF and WRF/STILT showed a large deviation in their simulation of the contribution of fossil fuel fluxes at one of the Ochsenkopf valleys.Further analysis on identifying the sources of these discrepancies suggests that the WRF/STILT predictions are highly sensitive to two model parameters -the vertical velocity variance (σ w ) and the Lagrangian time-scale (T L ) that were developed for coarse resolution fields; hence a further refinement of σ w and T L is required in STILT when driving with highresolution meteorological fields.However no firm conclusions can be drawn about the relative merits of different advection schemes used in the models.This calls for a more detailed inter-comparison study using different model resolutions since the inter-model or model-measurement mismatch shows strong sensitivity to spatial resolutions.
Nevertheless, the similarity of the results provided by WRF and WRF/STILT at high-resolution as well as the fact that the inter-model differences are a factor of two smaller than the model-observation differences and about a factor of three smaller than the mismatch between the current global model simulations and the observations, suggests the usefulness of STILT as an adjoint model of WRF.To achieve the definitive proof to justify the use of STILT as an adjoint of WRF, one would further need to carry out quantitative analysis of error characteristics between the models and to perform successful inversion using this model framework.

2. 3
Model domain and period of simulationsWRF-VPRM simulations of CO 2 and meteorological fields are carried out for the period from 2 to 30 August 2006, and for a single day in 2008 (20 October 2008) for a domain centered over Ochsenkopf in Germany (see Fig.2).The outer and inner domains have a total area of about 600 km × 600 km and 250 km × 250 km respectively.A period during summer (August 2006) is chosen as a case for the comparison of transport models because an increase of biological activity as well as strong variability of diurnal patterns of surface fluxes and mesoscale transport can be expected.The period on 20 October 2008, was chosen due to the availability of vertical profiles of CO 2 concentrations measured during an aircraft campaign(Pillai et al., 2011).These profiles provide a quality assessment on the performance of the transport models and also assist in finding the potential source of any model mismatch.Accordingly, WRF/STILT-VPRM simulations of CO 2 are carried out for different receptor locations corresponding to either different vertical levels of the Ochsenkopf tower or to the flight-track during the airborne measurement campaign at Ochsenkopf.

Figure 2 .
Figure 2. Model Domains showing topography (altitude in meters): The rectangle inside the domain indicates the boundaries of nested domain with 2 km × 2 km resolution.The domain outside the nested domain is with 6 km × 6 km resolution.The elevation data is from USGS elevation model-GTOPO30s-with spatial resolution of approximately 1 km.5The asterisk symbol denotes the position of a flat region which is described in Section 3.

Fig. 2 .
Fig. 2. Model Domains showing topography (altitude in meters): The rectangle inside the domain indicates the boundaries of nested domain with 2 km × 2 km resolution.The domain outside the nested domain is with 6 km × 6 km resolution.The elevation data is from USGS elevation model-GTOPO30s-with spatial resolution of approximately 1 km.The asterisk symbol denotes the position of a flat region which is described in Sect.3.

Figure 3 .
Figure 3.Comparison of observed and simulated CO2 concentrations (3-hourly averages) at 163 m over the Ochsenkopf tower site for August 2006.The magenta dotted line denotes the WRF/STILT-VPRM prediction when prescribing mixing height from WRF.

Fig. 3 .
Fig. 3. Comparison of observed and simulated CO 2 concentrations (3-hourly averages) at 163 m over the Ochsenkopf tower site for August 2006.The magenta dotted line denotes the WRF/STILT-VPRM prediction when prescribing mixing height from WRF.

Figure 4 .
Figure 4. Time-series of (a) mixing height (z i ) in meters with inset showing the diurnally averaged z i simulated by WRF-VPRM and WRF/STILT-VPRM for August 2006.(b & c) Inter-model comparison of Gross Ecosystem Exchange (GEE) and Respiration fluxes 5 (both are in the units of 2 1 mole/m s μ

Fig. 4 .Figure 4 .
Fig. 4. Time-series of (a) mixing height (z i ) in meters with inset showing the diurnally averaged z i simulated by WRF-VPRM and WRF/STILT-VPRM for August 2006.(b and c) Inter-model comparison of Gross Ecosystem Exchange (GEE) and Respiration fluxes (both are in the units of µmole m −2 s −1 ) simulated for the same period.The orange line denotes the one-to-one line.The bias, standard deviation of the difference (stdv) and squared correlation coefficient (R 2 ) are indicated in the figure panels.

Fig. 4 .Fig. 4 .
Fig. 4. Time-series of (a) mixing height (z i ) in meters with inset showing the diurnally averaged z i simulated by WRF-VPRM and WRF/STILT-VPRM for August 2006.(b and c) Inter-model comparison of Gross Ecosystem Exchange (GEE) and Respiration fluxes (both are in the units of µmole m −2 s −1 ) simulated for the same period.The orange line denotes the one-to-one line.The bias, standard deviation of the difference (stdv) and squared correlation coefficient (R 2 ) are indicated in the figure panels.

Figure 5 .
Figure 5.The altitude-distance cross-section showing CO 2 distribution around Ochsenkopf during the DIMO aircraft campaign on 20 October 2008: (a) Observations (b) WRF-VPRM, (c) WRF/STILT-VPRM and (d) aircraft track colored with flight altitude range (see the legend box labeled as "AGL"; units: km).In a-c: The time of the flight is given in purple at the top X-axis.In d: the cumulative distance (see the legend box labeled as "distance"; units: km) flown by the aircraft is labeled with the asterisk (*) symbol.

Fig. 5 .
Fig. 5.The altitude-distance cross-section showing CO 2 distribution around Ochsenkopf during the DIMO aircraft campaign on 20 October 2008: (a) Observations (b) WRF-VPRM, (c) WRF/STILT-VPRM and (d) aircraft track colored with flight altitude range (see the legend box labeled as "a.g.l."; units: km).(a-c): The time of the flight is given in purple at the top X-axis.(d): the cumulative distance (see the legend box labeled as "distance"; units: km) flown by the aircraft is labeled with the asterisk (*) symbol.

Figure 6 .
Figure 6.Total contribution of a pseudo-emission source on downstream concentrations as predicted by WRF/STILT (forward) under different model parameter set-ups and by the WRF model.CO 2 concentrations for a pressure-weighted atmospheric column from the surface to 190 m simulated by a) WRF/STILT (control run) b) WRF/STILT with free 5 tropospheric turbulence reduced to 1 cm s -1 c) WRF/STILT with free tropospheric turbulence reduced to 1 cm s -1 and the Lagrangian time scale reduced by a factor of 10: T L = 0.1*T L-original and d) WRF are shown.The "+" symbol (in magenta) denotes the source point and the "×" symbol (in magenta) denotes the approximate location of the valley (the aircraft location at 10 UTC) where a large CO 2 concentration was observed 10 (see Fig. 5).The colour-bar indicates CO 2 concentrations in units of ppm.

Figure 6 .
Figure 6.Total contribution of a pseudo-emission source on downstream concentrations as predicted by WRF/STILT (forward) under different model parameter set-ups and by the WRF model.CO 2 concentrations for a pressure-weighted atmospheric column from the surface to 190 m simulated by a) WRF/STILT (control run) b) WRF/STILT with free 5 tropospheric turbulence reduced to 1 cm s -1 c) WRF/STILT with free tropospheric turbulence reduced to 1 cm s -1 and the Lagrangian time scale reduced by a factor of 10: T L = 0.1*T L-original and d) WRF are shown.The "+" symbol (in magenta) denotes the source point and the "×" symbol (in magenta) denotes the approximate location of the valley (the aircraft location at 10 UTC) where a large CO 2 concentration was observed 10 (see Fig. 5).The colour-bar indicates CO 2 concentrations in units of ppm.

Fig. 6 .
Fig. 6.Total contribution of a pseudo-emission source on downstream concentrations as predicted by WRF/STILT (forward) under different model parameter set-ups and by the WRF model.CO 2 concentrations for a pressure-weighted atmospheric column from the surface to 190 m simulated by (a) WRF/STILT (control run) (b) WRF/STILT with free tropospheric turbulence reduced to 1 cm s 1 (c) WRF/STILT with free tropospheric turbulence reduced to 1 cm s 1 and the Lagrangian time scale reduced by a factor of 10: T L = 0.1 * T L−original and (d) WRF are shown.The "+" symbol (in magenta) denotes the source point and the "×" symbol (in magenta) denotes the approximate location of the valley (the aircraft location at 10:00 UTC) where a large CO 2 concentration was observed (see Fig.5).The colour-bar indicates CO 2 concentrations in units of ppm.

Figure 7 .
Figure 7. Sensitivity of inter-model or model-measurement mismatch to models' spatial resolution at Ochsenkopf tower at a measurement level of 90 m for August 2006: in terms of (a) squared correlation coefficient (R 2 ) and (b) monthly bias (WRF/STILT-VPRM minus WRF-VPRM or model minus measurement).The error bar in (b) denotes the 5 monthly standard deviation of the differences.

Fig. 7 .
Fig. 7. Sensitivity of inter-model or model-measurement mismatch to models' spatial resolution at Ochsenkopf tower at a measurement level of 90 m for August 2006: in terms of (a) squared correlation coefficient (R 2 ) and (b) monthly bias (WRF/STILT-VPRM minus WRF-VPRM or model minus measurement).The error bar in (b) denotes the monthly standard deviation of the differences.

Table 1 .
Overview of model set-up used in WRF.

Table 2 .
Summary statistics of inter-model and measurementmodel comparisons for different model levels (in meters) at the Ochsenkopf tall tower observatory for August 2006: Abbreviations: WRF-STILT: WRF-VPRM simulations minus WRF/STILT-VPRM simulations; WRF-STILT.zi wrf : WRF-VPRM simulations minus WRF/STILT-VPRM simulations using mixing height prescribed from WRF; Obs-STILT: Observations minus WRF/STILT-VPRM; Obs-WRF: Observations minus WRF-VPRM.mean: mean of the differences between models or between measurement and model (measurement minus model); sd: standard deviations of the differences between models or between measurement and model (measurement minus model); R 2 : squared correlation coefficient between models or between measurement and model.The model simulations are performed at a spatial resolution of 2 km × 2 km.

Table 3 .
Summary statistics of inter-model and measurement-model comparisons for different time scales of analysis.A tower level of 163 m at Ochsenkopf tall tower observatory for August 2006 is used for this analysis.See Table 2 caption for abbreviations.

www.atmos-chem-phys.net/12/8979/2012/ Atmos. Chem. Phys., 12, 8979-8991, 2012 8990 D. Pillai et al.: Comparing Lagrangian and Eulerian models for CO 2 transport boundary
layer top, as parameterized by the models.The complex terrain over Ochsenkopf does not seem to have a role on causing part of these inter-model mismatches since no better inter-model agreement is obtained for a flat region.