Biogeosciences Modelling regional scale surface fluxes , meteorology and CO 2 mixing ratios for the Cabauw tower in the Netherlands

We simulated meteorology and atmospheric CO2 transport over the Netherlands with the mesoscale model RAMS-Leaf3 coupled to the biospheric CO 2 flux model 5PM. The results were compared with meteorological and CO2 observations, with emphasis on the tall tower of Cabauw. An analysis of the coupled exchange of energy, moisture and CO 2 showed that the surface fluxes in the domain strongly influenced the atmospheric properties. The majority of the variability in the afternoon CO 2 mixing ratio in the middle of the domain was determined by biospheric and fossil fuel CO 2 fluxes in the limited area domain (640×640 km). Variation of the surface CO 2 fluxes, reflecting the uncertainty of the parameters in the CO 2 flux model 5PM, resulted in a range of simulated atmospheric CO 2 mixing ratios of on average 11.7 ppm in the well-mixed boundary layer. Additionally, we found that observed surface energy fluxes and observed atmospheric temperature and moisture could not be reconciled with the simulations. Including this as an uncertainty in the simulation of surface energy fluxes changed simulated atmospheric vertical mixing and horizontal advection, leading to differences in simulated CO 2 of on average 1.7 ppm. This is an important source of uncertainty and should be accounted for to avoid biased calculations of the CO2 mixing ratio, but it does not overwhelm the signal in the CO2 mixing ratio due to the uncertainty range of the surface CO2 fluxes. Correspondence to: L. F. Tolk (lieselotte.tolk@falw.vu.nl)


Introduction
Terrestrial carbon uptake is an important process in the global carbon cycle.It removes a substantial part of the anthropogenic emitted CO 2 from the atmosphere (Canadell et al., 2007).A useful method to increase our understanding of the terrestrial CO 2 fluxes is inverse modelling of atmospheric CO 2 mixing ratio observations (e.g.Gurney et al., 2002).In this method the atmospheric signal is used to constrain the surface fluxes using an atmospheric transport model.The results of inversion calculations depend to a large extent on the quality of atmospheric modelling (Stephens et al., 2007).
Therefore, correct simulation of the atmospheric transport, and accounting for the uncertainties, is an important goal in inverse modelling of CO 2 .Atmospheric transport is modelled at increasingly higher resolutions to capture the high spatial and temporal variability in observed CO 2 mixing ratios over the continent.Continental scale studies show that the forward simulation of CO 2 improved by increasing the horizontal resolution from a number of degrees (Gurney et al., 2002) to one degree or less (Geels et al., 2007;Parazoo et al., 2008).Further increasing the horizontal resolution to just a few kilometres in more limited domain studies (Dolman et al., 2006) was shown to improve the CO 2 mixing ratio simulation at observation stations in uneven and coastal terrain, because of the models ability to simulate mesoscale circulations, like sea breezes and topography induced katabatic flows (Nicholls et al., 2004;Riley et al., 2005;Van der Molen and Dolman, 2007;Sarrat et al., 2007;Ahmadov et al., 2009).This also avoids representation errors by resolving a larger part of the variability in the CO 2 mixing ratio (Corbin et al., 2008;Tolk et al., 2008).L. F. Tolk et al.: Modelling regional scale surface fluxes, meteorology and CO 2 Despite these achievements correct modelling of the CO 2 mixing ratios remains challenging.Model intercomparisons of global (Stephens et al., 2007;Law et al., 2008), continental (Geels et al., 2007) and mesoscale models (Van Lipzig et al., 2006;Sarrat et al., 2007) showed discrepancies in the meteorology and CO 2 modelled of different models.In the simulation of CO 2 mixing ratios both advection and entrainment play an important role (Vila et al., 2004;Casso-Torralba et al., 2008) and the quantification of uncertainties in these physical processes is one of the major questions in transport modelling.Comparisons at a coarser scale with observations showed that an erroneous simulation of the advection (Lin and Gerbig, 2005) and of vertical mixing (Gerbig et al., 2008) can lead to uncertainties in the simulated CO 2 mixing ratio of several ppm.Here we study at regional scale the effect of surface flux uncertainties on these transport errors.
In the present study a high resolution simulation is performed with the non-hydrostatic Regional Atmospheric Modeling System (RAMS; Pielke et al., 1992).The performance of the simulation is assessed with meteorological and CO 2 observations.We address a potential source of error in the simulated atmospheric vertical mixing: the simulation of the surface energy fluxes.Its uncertainty is estimated based on a comparison of different models (RAMS, WRF and ECMWF) and by using different parameter values in the surface flux model within RAMS.The model simulations are compared with both surface flux observations and with atmospheric CO 2 mixing ratio observations.We coupled the biospheric CO 2 flux model 5PM to the RAMS atmospheric transport model, in order to study the coupled exchange of energy, moisture and CO 2 .In this framework the impact of the surface energy fluxes on the simulation of atmospheric transport and consequently on the CO 2 mixing ratio is quantified.Novel in our approach is that we distil the impact of the uncertainty in the simulated surface energy fluxes on the atmospheric CO 2 mixing ratio, and that we quantify this CO 2 transport error in a Eulerian approach.
Also, the uncertainty in the CO 2 surface fluxes is addressed.With the coupled RAMS-5PM simulation system these are propagated into a range of CO 2 mixing ratios.This indicates the minimal performance of the atmospheric transport model required for the use in inversion studies, since the uncertainty in the transport modelling should not exceed the uncertainty related to CO 2 surface flux uncertainty.The parameters in the biospheric model 5PM have been optimized in a previous study for a number of eddy correlation flux observations (Groenendijk et al., 2009).Innovative in this study is that we show at mesoscale a realistic uncertainty range of CO 2 mixing ratios due to uncertainties in the CO 2 surface fluxes, based on independently determined a-priori flux estimates.
Finally, we separate the contribution of different CO 2 sources and sinks to the CO 2 mixing ratio at Cabauw, i.e. the influence of the advection of CO 2 into our domain (back-ground contribution), the fossil fuel emissions, sea-air CO 2 exchange, and terrestrial respiration and assimilation fluxes.The relative importance of the different CO 2 contributions indicates which uncertainties in the surface CO 2 fluxes are important and which can be neglected.Additionally, the relative contribution of the near field versus the far field fluxes on the CO 2 mixing ratio is shown, another important subject in regional scale inverse modelling (Zupanski et al., 2007;Lauvaux et al., 2008;Gerbig et al., 2009).
The paper is organized as follows: in Sect. 2 the simulation set-up is described, in Sect. 3 the performance of the model is validated against meteorological observations, Sect. 4 describes the simulated CO 2 fluxes and mixing ratios compared to observations and in Sect. 5 the implications of our results for the interpretation of the observations, future forward and inverse CO 2 simulations are discussed.

Simulation period and domain
We performed simulations with the Regional Atmospheric Modeling System (RAMS) for 22 days in June 2006.In this time of year the biogenic assimilation fluxes during daytime of CO 2 were large.This period was selected because it covers a number of meteorological regimes with different wind directions and frontal passages, influencing the atmospheric properties and carbon exchange.South-easterly winds coincided with clear sky conditions, while northerly and southwesterly winds caused more cloudy conditions.
A two way nested grid was used (Fig. 1) centred on the Netherlands at 52.25 • N and 5.2 • E, with a 320×320 km domain at 4 km resolution nested in a 640×640 km domain at 16 km resolution (Table 1).The centre of the domain is relatively flat, with a maximum elevation of ∼100 m.The southeast part of the domain has more orography, up to ∼500 m.The dominant land use types in the area are cultivated lands (crops and grasslands) and urban areas.Large cities and industrial areas of the Netherlands, Belgium and the German Ruhr Area are within the domain.To the north and the west the Netherlands is bounded by the North Sea.

Simulation setup
The atmospheric simulations were performed with the nonhydrostatic mesoscale model RAMS (Pielke et al., 1992), which has already been used to simulate the behaviour of CO 2 in the atmosphere in a number of studies (e.g.Denning et al., 2003;Nicholls et al., 2004;Sarrat et al., 2007;Ter Maat and Hutjes, 2008).The version used in this study is B-RAMS-3.2,including adaptations to secure mass conservation (Medvigy et al., 2005;Meesters et al., 2008).For vertical diffusion, we built in the turbulence scheme as described in Sect. 2 of Hong and Pan (1996) (corrections: division by h should be added in their Eq.( 4), and removed in their Eq.( 9)).This scheme uses a non-local turbulence parameterization within the convective boundary layer.Inclusion of this within the medium range forecast (MRF) model has been shown by Troen and Mahrt (1986), Holtslag et al. (1995) and Hong and Pan (1996) to simulate the daytime boundary layer structures more realistically than local mixing schemes.In our simulations with RAMS we applied this non-local scheme also to the CO 2 -transport.Cumulus convection was not parameterized in the simulations.The surface energy fluxes were simulated using Leaf-3 (Walko et al., 2000).The vegetation Leaf Area Index (LAI) of the MODIS database was used.
Meteorology, soil temperature and soil moisture were initialized with ECMWF analysis data (Uppala et al., 2005).In order to be consistent with the RAMS soil wilting point (wp) and field capacity (fc), the ECMWF soil moisture (η) was scaled towards RAMS soil variables based on a soil wetness index (SWI): Optimized CO 2 mixing ratio fields at 1×1 • resolution from CarbonTracker Europe (Peters et al., 2009) were used for initial and boundary conditions of the CO 2 mixing ratio.The simulations were nudged every 3 h to CarbonTracker CO 2 mixing ratios and every 6 h to the ECMWF analysis meteorology with a nudging relaxation time scale of 900 s.The nudging extended inward from the lateral boundary by 5 grid cells and the centre of the domain was free of nudging.CO 2 fluxes from fossil fuel burning were included in the simulations based on the IER database at 10 km resolution (http: //carboeurope.ier.uni-stuttgart.de/).The CO 2 fluxes from the coastal sea inside the domain were calculated based on climatologic estimates of the partial pressure of CO 2 in the sea (Wanninkhof, 1992;Takahashi et al., 2002).Biospheric CO 2 surface fluxes were modelled with 5PM (Groenendijk et al., 2009).This model was coupled to RAMS through the radiation, the temperature and humidity of the canopy air and influences the CO 2 mixing ratio at the lowest atmospheric level.The CO 2 assimilation does in this model not depend on the energy fluxes of RAMS through the stomatal conductance (Collatz et al., 1991) or on Leaf Area Index (LAI) (Sellers et al., 1996).In 5PM the photosynthesis is calculated following Farquhar et al. (1980), where photosynthesis is either limited by the carboxylation rate, which is enzyme limited, or by the light limited RuBP regeneration rate.The most important assimilation parameters in this model are the maximum carboxylation capacity (Vc max ) and the light use efficiency (α).Respiration was calculated with the relationship by Lloyd and Taylor (1994): where R 10 is the respiration rate at a reference temperature of 10 • C, E 0 is the activation energy divided by the universal gas constant, T 0 is a constant of 227.13 K and T is soil temperature.For further specifications of 5PM see Groenendijk et al. (2009).Groenendijk et al. (2009) optimized the parameters of this model (Vc max , α, R 10 and E 0 ) for the full canopy based on a large number of Fluxnet observations (Baldocchi et al., 2001).We applied parameter values optimized for the temperate zone, for the period of May-July for all years (Table 2).The parameter values used in our simulations were kept constant in time.Simulations were performed with CO 2 fluxes calculated based on the best guess parameter values.For respiration and assimilation of the most abundant vegetation species (crops and grass) we also simulated fluxes using the upper and the lower parameter values within the standard deviation of the parameter estimate.In this way a range of CO 2 mixing ratios was simulated based on the different CO 2 flux parameter settings.In the rest of this work, we report the range of uncertainties, i.e. the difference between the highest and lowest values in our set of simulations.Further specifications on the design of the simulations are given in Table 1.

Observations
A large number of observations of the atmospheric properties and the surface fluxes were available for model validation.Data from continuous CO 2 mixing ratio measurements, performed by a Licor 7000 with a precision of 0.05 ppm, and  Groenendijk et al. (2009).Vc max is the full canopy maximum carboxylation capacity, α the light use efficiency for the full canopy, E 0 / is the respiration activation energy divided by the universal gas constant, R 10 is the respiration rate at 10 • C and n sites indicates the number of sites used in the optimizations of the parameters.Where uncertainty ranges are shown the best guess, upper and lower estimates of the parameters are used in the simulations.In between brackets the parameter values that returned the best CO 2 mixing ratios. Assimilation meteorological data from the tall tower at Cabauw at a height of 20 m, 60 m, 120 m and 200 m were used.Also, atmospheric observations for temperature, humidity, wind speed and direction were available at 110 synoptic 2 m stations over the Netherlands and from the radiosondes that were released twice a day at De Bilt, which is about 25 km north-east of the Cabauw site.Observations of the surface fluxes were available for sensible heat, latent heat and CO 2 fluxes from eddy correlation measurements (Aubinet et al., 2001;Dolman et al., 2002;Wilson et al., 2002;Jacobs et al., 2007;Braam, 2008;Aubinet et al., 2009).Additionally, scintillometer measurements provided extra sensible heat flux measurements over a horizontal path of 0.35-5 km (De Bruin et al., 2004).The locations are specified in Fig. 1 and in Table 3.
3 Results: Meteorological performance of the model

Consistency of the simulation in time
The simulated period of 22 days covered a number of different weather regimes with different wind directions, solar radiation and temperature (Fig. 2).East to south-easterly winds brought clear weather, with increasing temperature and relative large day-night temperature amplitudes, like at 9-14 June 2006 (doy 160-165).In general, northern and Here we show the results of the standard simulation used in this study.For this simulation, the standard RAMS settings have been modified to obtain a more realistic Bowen ratio, as will be described in Sect.3.2.The model reproduced the synoptic variations over the full 22-day period without any re-initialization of the simulation (Fig. 2 and Table 4).This was achieved by prescribing boundary conditions from reanalysis products such as ECMWF meteorology and Car-bonTracker CO 2 mixing ratios.A change in the large-scale atmospheric situation was thus passed on to the inner domain for which RAMS simulated, mimicking the effect of a re-initialization.A large advantage of not needing to reinitialize the RAMS model over multi-week periods is that mass continuity of tracers and a balance of the physical equations for energy and water was ensured.
A comparison of the statistics for the first and last half of the period showed the consistency of the model performance in time (Table 4).Hourly temperature (T) and humidity (q) were simulated comparably well in both periods.Radiation showed a better performance in the first half of the period, which can be attributed to the occurrence of clouds in the second half of the period, rather than to a drift of the simulated meteorology with time.Incoming solar radiation and its reduction by clouds was mostly simulated with the correct amplitude and frequency.However, the exact location of the clouds and subsequently the timing of the radiation reduction sometimes deviated from the observations, as was also seen in similar mesoscale model studies (Denning et al., 2003;Van Lipzig et al., 2006;Parazoo et al., 2008).

Uncertainties in the surface energy fluxes
Surface energy fluxes are important drivers of processes in the atmosphere, influencing amongst others the atmospheric T, q and vertical mixing.Uncertainties in the surface energy fluxes thus may be an important source of uncertainty in the simulation of the atmospheric properties and are addressed in this section.
We compared the simulated energy fluxes with eddy covariance and scintillometer measurements (Fig. 3) and made a comparison with two other models: WRF and ECMWF.Additionally, we studied the sensitivity of the simulated sensible (H) and latent (LE) heat fluxes to changes in the surface flux calculation by Leaf3, and its effect on the atmosphere.This revealed that the simulations at some days either captured the observed surface energy fluxes, or the observed T and q vertical profile in the planetary boundary layer (PBL), but could not reconcile both.The atmospheric observations and the comparison with other models suggested a higher Bowen ratio (i.e. the ratio of sensible to latent heat flux, β) than simulated with standard Leaf3 settings.Possible options for such an increased β are shown in Table 5 and described below, with a focus on T and H because of their importance for atmospheric vertical mixing.
The energy fluxes showed a large variation between low (crops, grasslands) and high (forest) vegetation types (Fig. 3, note the different scale on the y-axis for low and high vegetation).These differences were driven by differing vegetation characteristics such as the low aerodynamic resistance and stomatal conductance in forest and were reproduced in the RAMS simulations.
With the standard Leaf3 vegetation characteristics (green line in Fig. 3) most of the eddy correlation observations were captured reasonably well.The scintillometer observations at Maas-en-Waal and Haarweg were slightly underestimated.With these settings the PBL T was underestimated and q overestimated at clear days with eastern winds, in comparison with radiosonde observations (Fig. 4 and Table 5), the Cabauw tall tower and synoptic 2 m observations (not shown).
We compared our findings with the ECMWF (Uppala et al., 2005) forecast simulations and with a simulation performed with WRF (Skamarock et al., 2008) using the MRF PBL scheme and either ECMWF or NCEP boundary conditions for the same domain and period.Both simulations Table 4. Statistics of the simulation, in comparison with observations of the potential temperature (T) and CO 2 at 20 m and 200 m, humidity (q) at 2 m and the incoming shortwave radiation (rshort) at Cabauw, for the full period, the first and the second half of the simulation period.Humidity is compared to the simulated canopy air humidity, CO 2 and T with simulated atmospheric values.

Full period
First half Second half matched the observed T well (not shown), but also failed to match the observed H for grass and crops (light and dark blue lines in Fig. 3).
Sensitivity tests showed that the strength of LE and H depended strongly on (1) the water availability, (2) the minimal stomatal resistance, which determines the plants' resistance to transport of water and CO 2 and (3) the fraction of the surface that is vegetated (Table 5).
(1) Decreasing the soil moisture content led to a strong increase in β, because of the decrease in evaporation from the soil and transpiration from the plants.(2) Also a doubling of the minimal stomatal resistance, a rather uncertain parameter in Leaf3 (Walko et al., 2000), for grass and crops from 100 sm −1 to 200 sm −1 increased β. (3) With standard Leaf3 settings the vegetation fraction of grass and crops did not exceed ∼70%, even with a high LAI.We increased the vegetation fraction to 90% when the LAI is larger than 1, in line with for example the settings in the ECMWF surface model.This increase led to a considerable increase in H and the atmospheric T (Table 5).
Simulations with these adapted settings returned an overestimation of H at most of the low vegetation sites (e.g. with adaptations (2) and (3), red line in Fig. 3), but returned a rather well simulated PBL T and q (red lines in Figs. 2 and 4).We will refer to the simulation with increased stomatal resistance and vegetation fraction and standard soil moisture as the "high β simulation", while the simulation with standard Leaf3 vegetation characteristics will be referred to as the "low β simulation".These two cases however do not capture the full range of uncertainty in beta, as even the low beta simulation sometimes still overestimates H and underestimates LE (Fig. 3 and Table 5).Decreasing for example the stomatal resistance as suggested for Cabauw by Jackson et al. (2003) would further reduce β in such situations (Table 5).The uncertainty in surface energy fluxes may thus even be larger than inferred here from the atmospheric observations, indicating that our uncertainty estimates may be conservative.
We tested the sensitivity of our findings to the moment of initialization.When initialized 6 h, 18 h or over 5 days in advance, the simulated noon temperatures deviated 2.0, 2.5 or 3.1 K from the observations, respectively (Table 5) using the settings yielding the low β.Hence, the largest part of the temperature underestimation built up within a few hours and was rather independent of the moment of initialization.For the high β simulation the results were robust to a change of the moment of initialization (not shown).
The fact that the surface flux observations and the atmospheric observations both suggest different optimal β's indicated an uncertainty in what the correct β should be in the simulations for the full domain.Further discussion on this will be presented in Sect. 5. We will use the simulation with the high β (and hence lowest T and q bias in the atmosphere) to investigate the structure of the simulated CO 2 fields.This simulation corresponds to the standard run that was discussed in Sect.3.1.In the next section the effects of the uncertainty in the surface energy fluxes on the atmospheric vertical mixing are addressed.

Atmospheric temperature and humidity profiles
The results of the simulations with (1) low β and (2) high β were compared with the radiosonde observations in De Bilt (e.g.Fig. 4).Generally a well-mixed PBL developed that has a lower potential temperature and is moister than the free troposphere.At clear nights cooling near the surface led to a shallow (200 m) and stable PBL.
Simulations with the MRF turbulence scheme showed a better performance than the standard RAMS Mellor Yamada turbulence scheme (Fig. 4).Still, the height of the PBL was not always captured correctly and the jump of T and q was less pronounced than observed.Increasing the vertical resolution to 60 instead of 25 vertical layers in the lower 3 km of the atmosphere did not change this.This may be due to the limited horizontal resolution of 4 km, or to uncertainties in the parameterization of vertical transport, like lack of subsidence, too much entrainment or an incorrect free tropospheric lapse rate of T or q.The free tropospheric values are generally assumed reasonable due to the use of ECMWF analysis boundary conditions.
Simulation of the nocturnal PBL is even more challenging.The atmospheric stability during clear nights was systematically underestimated by the model, which is a common feature for most atmospheric transport models (Geels et al., 2007;Gerbig et al., 2008).In the simulation the nocturnal surface signal reached up to ∼400 m while in the observations it was limited to ∼200 m.The poorly simulated height of the PBL will cause discrepancies in the mixing ratios which are not directly related to the magnitude of the CO 2 fluxes.
The depth of the PBL is amongst others driven by the surface sensible heat fluxes.Uncertainty in these, as described in Sect.3.2, will thus lead to uncertainties in the PBL depth.At days when the PBL height was clearly defined the rootmean-square error (RMSE) of the noon PBL height was for In black the observations.The red band is the range of CO 2 fluxes simulated with a spread of biosphere model parameters (Table 2).Yellow shows the simulation that best fits the CO 2 mixing ratio observations.both simulations ∼350 m.The simulation with a relative low β (1, green line in Fig. 4) showed a mean bias of ∼−100 m, while the simulation with a higher β (2, red line in Fig. 4) had a positive mean bias of ∼75 m.The uncertainty in the surface fluxes can thus explain part of the uncertainty in the simulated PBL height as for example indicated for ECMWF simulations in Gerbig et al. (2008).

Results: CO 2 fluxes and mixing ratios
The framework of the simulated meteorology as described in the previous sections allows us to study the coupling between the surface CO 2 fluxes and the atmospheric CO 2 mixing ratios.First, we will address the uncertainties in the CO 2 fluxes.Secondly, we will show how these propagate into a range of simulated CO 2 mixing ratios which is compared to the observations at the Cabauw tall tower.Finally, the contribution of the different CO 2 fluxes and the background CO 2 to the total CO 2 mixing ratio is unravelled.

CO 2 flux variability
The parameters optimized in the biosphere model 5PM showed a rather large variability in time and space, which is reflected in the uncertainty in the average CO 2 flux parameter values (Table 2).The range of CO 2 fluxes (Fig. 5) simulated with these parameter settings, which were optimized based on European observations for the temperate zone, agreed well with observations and literature values (Jacobs et al., 2007) for our much smaller domain.
CO 2 fluxes simulated with the parameters that gave an unbiased result compared to the observed CO 2 mixing ratio at 200 m at Cabauw are indicated in yellow in Fig. 5  Fig. 6.CO 2 concentration at 4 heights at Cabauw tower.In black the observations.The red band is the range of CO 2 mixing ratio simulated with a spread of CO 2 flux parameters (Table 2), it reflects the effect of uncertainties in the surface CO 2 fluxes on the CO 2 mixing ratios.brackets in Table 2.Note that in a follow-up study we intend to formally estimate the parameters of the 5PM model based on the forward results presented here, rather than simply selecting an unbiased set of values.
Respiration fluxes show a soil temperature driven diurnal and synoptic variation, which is represented by the model (Fig. 5a).Additionally, the observed respiration fluxes showed differences between sites in magnitude and diurnal amplitude, not observed in T and not included in the model.These kind of variations are meant to be implicitly included in the parameter uncertainties.The observed respiration at Cabauw, Horstermeer and Lonzee are in the lower part of the uncertainty range, while the observations at Lutjewad are often near the top of the range.
Assimilation fluxes were inhomogeneous in space and time as well (Fig. 5b).Generally, the uptake by crops was higher than by grass, as was correctly simulated.Also the observed assimilation reduction at days with limited radiation was captured by the model (within the limitations of the RAMS radiation calculations).Observed assimilation at Lutjewad is relatively high and near the maximum of the uncertainty range, while Lonzee and Cabauw are near the minimum and Horstermeer is in the middle of the range.
Both simulated respiration and assimilation at Lonzee were in the beginning of the period outside the uncertainty range from observations, most likely because the CO 2 flux calculations were not LAI dependent.To overcome this the biosphere model should be extended with spatially explicit, time-varying LAI (e.g.Sellers et al., 1996), which should also give a better representation of the spatial variability of the assimilation fluxes.
Our uncertainty estimates agreed with variability in Dutch grass sites estimated by Jacobs et al. (2007).The standard deviations of their respiration parameters were almost the same as in our study.Their photosynthesis model, and therefore these parameters are different, but the range is also comparable to the range used in this study, with 20-60% variations on the parameters in the GPP model.They concluded that within small regions with relatively uniform climatic conditions the variability may be similar to the one observed at European larger scales, as is applied here.

CO 2 mixing ratios
Uncertainties in the CO 2 respiration and assimilation fluxes had a significant influence on the CO 2 mixing ratio when the air had passed land areas.The different CO 2 flux parameter settings returned a range of simulated CO 2 mixing ratios (Fig. 6).This range varied in time between 1 ppm and 25 ppm, with a 22-day average of 11.7 ppm in the well mixed afternoon PBL (at 200 m height).Small ranges were related to northern or south-westerly wind directions (Fig. 2  when the air predominantly originated from over sea and the land signal is suppressed, while a broader range occurred when the continental signal was large, i.e. with south-easterly wind, low wind speeds or frontal passages.This broad range indicates that the atmospheric mixing ratio potentially contains much information about the CO 2 fluxes within the simulation domain.This agrees with studies from Lauvaux et al. (2008) and Zupanski et al. (2007) and shows the potential for future inversion on this temporal and spatial scale.
As mentioned previously, the uncertainties in the simulation of the meteorology discussed in Sect. 3 give rise to uncertainties in the simulated CO 2 transport.Here we give an overview of the impact of those features on the simulation of the CO 2 mixing ratios.
The uncertainty in the surface energy fluxes (Sect.3.2) and consequently the vertical mixing (Sect.3.3) results in an uncertainty in the CO 2 mixing ratio.To quantify this the results of the simulations with relative low (1) and high (2) simulated β, were compared for CO 2 (Fig. 7).These two simulations (with the same CO 2 flux parameter settings) returned a difference in the afternoon CO 2 mixing ratio of on average 1.9 ppm.
This was to a small extent due to changes in the CO 2 fluxes between the two simulations.The RMSE of the total flux change was 0.69 µmol m −2 s −1 .Changes in the respiration caused by the temperature difference (RMSE=0.63 µmol m −2 s −1 ) were systematic.Assimilation fluxes slightly changed due to changes in cloud cover formation (RMSE=0.29 µmol m −2 s −1 ), and their difference over the full period was near to zero.Hence, the effect of respiration change on the CO 2 mixing ratio was relatively large compared to the assimilation change.Compensating for this by adjusting the R 10 left a difference of 1.7 ppm.This uncertainty in the CO 2 mixing ratio is caused by the difference in vertical mixing and horizontal advection, due to the surface flux uncertainty.
Other difficulties were related to the simulation of the nocturnal CO 2 mixing ratio, which is up to now, because of the known large uncertainties in the transport model vertical mixing schemes in stable conditions, not used in inversion studies (Gurney et al., 2002;Stephens et al., 2007;Geels et al., 2007).Our simulations confirm that the simulation of nocturnal CO 2 is biased, because of the simulation of a too deep nocturnal PBL (see Sect. 3.3).The absolute nocturnal CO 2 mixing ratio accumulation was not simulated correctly, leading to a low R 2 of the CO 2 mixing ratio time series at 200 m (Table 4) and simulated mixing ratios at 60-200 m that were during some nights totally outside the range of simulated CO 2 mixing ratios.As such it is clear that the representation of the nocturnal boundary layer in mesoscale models requires improvement.
Nevertheless a significant part of the diurnal variation, especially at lower sample levels, is captured by simulations.During the nights CO 2 accumulates near the surface and mixing ratios of over 450 ppm were seen at 20 m (Fig. 6).In the morning the PBL became unstable and the signal of the lower layers was mixed onto higher levels.This was for example reflected by early morning peaks at 200 m height, in the observations as well as the simulations.During a number of nights a sudden drop of simulated and observed mixing ratios was seen in the mixing ratio at 20 and 60 m around midnight, this is most likely due to the formation of a low level jet (e.g.Bosveld et al., 2008)  uncertainty is the simulated location and the timing of the cloud cover (Sect.3.1).At partly clouded days, this led to an error in simulation of the CO 2 fluxes and the depth of the PBL, and consequently biases in the CO 2 mixing ratio.Therefore, the exact simulated mixing ratios at these days should be regarded as more uncertain.Frontal passages may cause a comparable misrepresentation of the simulated concentrations (e.g. at 25 June; doy 176).

Different tracer signals at Cabauw
To study the relative importance of the different sources and sinks influencing the CO 2 mixing ratios we separated the simulated CO 2 mixing ratios at Cabauw into contributions from (a) background CO 2 entering through the lateral boundaries and (b) fluxes from within the RAMS domain.The latter were further separated into contributions from assimilation and respiration of different vegetation types, sea-air fluxes and fossil fuel emissions, which were included in the simulation as separate atmospheric tracers (Figs. 7 and 8).Assimilation and respiration fluxes were important in determining the CO 2 mixing ratio.They had an average influence during the day of −10.5 and 7.8 ppm respectively, with peaks up to ∼30 ppm at 200 m (Figs. 7 and 8a).In the nocturnal PBL the respiration tracer showed peaks up to ∼60 ppm at 20 m (not shown).The total biosphere signal, i.e. the sum of the respiration and assimilation, was because of the cancelling opposite signs more modest and had an average influence of 3.2 ppm.
Generally, the crops tracer was the most abundant assimilation tracer, even though Cabauw is a grassland site (Fig. 8b).This was due to the relatively high assimilation of crops compared to grass (Fig. 5b) combined with the large area of crops in the domain (Fig. 1).
In our domain the fossil fuel was also an important contributor to the total signal.Plumes with high fossil fuel tracer concentrations originating from the industrial sources moved over the domain.The major sources were the Ruhr Area, southeast of the Netherlands, the ports in the southwest of the Netherlands and in Belgium, and smaller diffuse sources found over the total domain.Afternoon values in the well mixed PBL varied between 2-8 ppm and the average mixing ratio of the fossil fuel fluxes was almost as large as the biospheric signal whereas its afternoon variance is half that of the biospheric signal (Fig. 8c).
The contribution of the sea-air CO 2 fluxes to the total signal was very limited at these time scales (Fig. 8a and  c), although on a global scale the exchange of CO 2 by the oceans plays an important role.The average flux of ∼0.02 µmol m −2 s −1 at the North Sea was a negligible compared to the continental fluxes in the domain.
CO 2 mixing ratios entering the domain through the lateral boundaries showed a strong diurnal cycle over land, and a much smaller cycle over the sea.When the wind direction was from the east or south the low daily continental values in the background mixing ratio, caused by assimilation over the continent outside of our simulation domain, reached Cabauw.This was for example seen at 12 and 13 June (doy 163 and 164) when the background concentration was reduced by ∼10 ppm.Because the influence on the mixing ratios in the middle of the domain can be considerable, the quality of the boundary conditions in a limited domain simulation is important.

Discussion and conclusions
We simulated three weeks in June 2006 with the B-RAMS-3.2mesoscale model at 4 km resolution.The simulations were able to reproduce the observed time series of the meteorological variables and CO 2 mixing ratios satisfactorily for most of the three week period.The model performance showed no drift and comparison with data remained acceptable throughout the full simulation.We found only limited sensitivity of the model performance to the moment of initialization.This, combined with the absence of significant drift, shows the possibility of a non-stop simulation without divergence of the results from the observations.Hence, in our simulations a re-initialization of the meteorology was not necessary.The advection of the ECMWF meteorological and CarbonTracker CO 2 mixing ratio boundary conditions over the domain prevents a runaway of the results in a dynamical and continuous way.This opens the way towards seasonal inversions at a high resolution.
The simulations with the RAMS MRF scheme showed a better performance than the standard Mellor Yamada scheme.This is in line with the findings of Holtslag et al. (1995) and confirms the importance of the parameterizations of turbulence and entrainment (Vila et al., 2004;Casso-Torralba et al., 2008) for the simulation of the atmospheric profiles.
Also a realistic simulation of surface energy fluxes is important for the simulation of atmospheric transport.Comparison with other models (WRF, ECMWF) and observations revealed a discrepancy between the simulations and the surface observations on the one hand and the atmospheric observations on the other hand.For a number of days the observed T profile could not be reconciled with H observations, something also seen in previous studies (e.g.Holtslag et al., 1995;Ek and Holtslag, 2004).At those times, the observed PBL height could also not be reproduced with a simple mixed layer model (Vila and Casso-Torralba, 2007) based on the observed H.
The atmospheric observations suggests that the total Bowen ratio (β) over the full domain may be higher than indicated by the surface observations.This may be due to heterogeneity of the energy fluxes within one vegetation type (e.g.Baldocchi et al., 2001) which is not included in our land surface scheme.Also, the limited amount of observations over grass and crops, especially in the eastern part of the domain, and energy balance closure problems (Wilson et al., 2002) of ∼30% for the Cabauw surface flux observations (Braam, 2008) may add uncertainty to the total surface energy flux estimate over the domain.Moreover, (freestanding) houses, trees and roads may lead to a different domain averaged β than observed with the surface observations over vegetated terrain only.Hence, optimizations for a single site (e.g.Jackson et al., 2003) may lead to another estimate of the surface energy fluxes than optimizations with atmospheric properties, which reflect larger scale processes and fluxes over a larger area (e.g.Uppala et al., 2005).
The uncertainty in the sensible heat flux adds uncertainty to the simulation of the vertical mixing, i.e. the simulated shallow and deep convection and the PBL height, and consequently horizontal advection.A comparison of simulations with standard Leaf3 settings with a relative low β (suggested by the surface observations) and with adjusted settings with a relative high β (suggested by atmospheric observations and the ECMWF and WRF simulations) was made.
The surface energy flux and the resulting atmospheric transport uncertainty caused a difference in the simulated CO 2 mixing ratio of on average 1.7 ppm.The estimate of the surface flux uncertainty used in our work is conservative, and represents a minimal value to take into account in future inversion studies.It is in the same order of magnitude as other estimated CO 2 modeling errors, such as due to misrepresentation of smaller scales (0.5-3 ppm; Van der Molen and Dolman, 2007;Tolk et al., 2008), and much larger than the measurement accuracy.
The change in energy fluxes led to a difference in the noon PBL height of ∼22%.A comparable sensitivity of the PBL height to changes in the surface energy fluxes (20%) was found by Vila and Casso-Torralba (2007) in a simple mixed layer study for Cabauw.Gerbig et al. (2008) showed at a lower model resolution that uncertainty in the PBL height of 40% gave a uncertainty of 3.5 ppm in the CO 2 mixing ratio.The ∼20% change in PBL height and consequent 1.7 ppm range due to surface energy flux uncertainty found here can thus explain a part of those errors.
As a consequence of the change in turbulent mixing, the horizontal advection also changed.The wind change had a standard deviation of 0.85 ms −1 and 0.98 ms −1 in u and v direction, respectively.This is about 40% of the random error found by Lin and Gerbigat 80 km resolution.It suggests that the simulated change in CO 2 mixing ratios is importantly influenced by changed advection resulting from the changes in turbulent mixing.
Hence, changing horizontal and vertical transport due to uncertainties in surface energy fluxes can explain an important part of the errors found in previous studies at a coarser resolution.To avoid biased CO 2 mixing ratio estimates, a comparison with observations of the simulated wind direction, wind speed and PBL height and, if needed, adjustment of the surface fluxes like described in Sect.3.2 is recommended as first step in an inversion.
The surface CO 2 fluxes in the domain strongly influenced the simulated CO 2 mixing ratios at Cabauw, causing by far the most of afternoon CO 2 mixing ratio variability (Fig. 8c).A realistic variation of parameter settings for the calculation of CO 2 fluxes (Groenendijk et al., 2009) resulted in a range of simulated CO 2 mixing ratios that spans on average 11.7 ppm.This atmospheric signal will in future inversion studies be used to constrain the surface CO 2 fluxes.The rather broad range indicates the potential for inversions, even though transport errors are in the order of several ppm.It confirms, with a complementary approach, the findings of Lin and Gerbig (2005) and Gerbig et al. (2008) and is in line with previous studies that stress the importance of the near field fluxes for the CO 2 mixing ratios over the continent (Zupanski et al., 2007;Lauvaux et al., 2008;Gerbig et al., 2009).
Although grass is the dominant vegetation type near Cabauw (Fig. 1) its atmospheric CO 2 mixing ratio is strongly affected by the assimilation of crops.This is due to the large magnitude of crops assimilation, the large area covered with crops, advection of the atmospheric signal and at wind still days entrainment from the residual boundary layer (Casso-Torralba et al., 2008;Vila et al., 2004).Hence, the information in the atmospheric mixing ratio measurements is not www.biogeosciences.net/6/2265/2009/Biogeosciences, 6, 2265Biogeosciences, 6, -2280Biogeosciences, 6, , 2009 limited to the very local fluxes within the nearest tens of kilometres.We conclude therefore that the scale of tens to hundreds of kilometers is convenient for future inversions of the atmospheric CO 2 mixing ratio signal at the tall tower of Cabauw.
Another important contributor to the CO 2 mixing ratio at Cabauw are the fossil fuel CO 2 fluxes.At these small scales uncertainties in the timing and exact location of the fluxes are important.The general assumption in global inversions that the fossil fuel fluxes are well known (Gurney et al., 2002) may be true for aggregated values in space and time (annual country totals) but is certainly not true for the scales in time and space that are modeled here.Because of the relative importance of the fossil fuel atmospheric signal, uncertainties in the timing and magnitude fossil fuel fluxes should be quantified and taken into account in future regional inversion studies.
The time series of grass and crops assimilation tracers as simulated for Cabauw were correlated (r=0.67).This suggests that the ability of the observed atmospheric CO 2 to distinguish between these vegetation types will be limited.The respiration and assimilation flux signals cancel each other during the day (correlation=−0.80),providing a relatively modest net atmospheric signal that may not constrain the two fluxes separately.During the night, when assimilation stops, the contribution of the respiration tracer to the total CO 2 mixing ratio becomes large compared to the contribution of the assimilation tracers.Potentially, nocturnal mixing ratios will be able to provide us therefore with a constraint on the division between respiration and assimilation (Ahmadov et al., 2009).
However, simulation of the nocturnal PBL is an important and long known source of uncertainty (e.g.Geels et al., 2007).The stability of the atmosphere at clear nights was systematically underestimated in our simulations which will lead to biased CO 2 mixing ratio estimates.Before simulated nocturnal CO 2 mixing ratios can be used in inversion studies they must at least be corrected for the PBL height.More importantly, the simulation of the nocturnal PBL should be improved as indicated for example by Steeneveld et al. (2008).The ability of the simulations to capture variations in the atmospheric stability at small temporal scales is promising.Because of the potential high value of the nocturnal mixing ratios in separating assimilation and respiration fluxes we plan to focus more work on this issue in the future.
Summarizing, the influence of the surface CO 2 and energy fluxes on the simulated atmospheric CO 2 mixing ratio, the temperature and humidity is large, especially at days with a continental footprint.This shows that atmospheric observations potentially contain much information about these fluxes at the scale of our simulation, i.e. at a spatial scale of tens to hundreds of kilometres.Most of the variability in the CO 2 mixing ratio is caused by fluxes within the domain, mainly by biospheric fluxes.Also the fossil fuel CO 2 fluxes play a role and their uncertainty should be taken into account in in-versions for such an urbanized and industrialized area.Difficulties identified in the simulation of CO 2 mixing ratios that reduce the information content of the simulated mixing ratio are (a) the systematic underestimation of the stability of the nocturnal PBL at clear nights which may lead to a biased CO 2 estimate, (b) incorrect timing of cloud formation, (c) uncertainty in the diurnal PBL height due to uncertainties in the parameterization of vertical transport and (d) the uncertainties in the driving of the atmospheric mixing by the surface energy fluxes.We quantified the latter and show it is with ∼1.7 ppm an important source of uncertainty in the CO 2 mixing ratio in the afternoon well-mixed PBL.Besides these shortcomings the atmospheric mesoscale simulation was shown to simulate the meteorological situation over the Netherlands non-stop for three weeks with reasonable accuracy.This, combined with the large simulated range of atmospheric CO 2 mixing ratios due to the spread in the CO 2 flux parameter settings provides a promising starting point for future inversion studies at the mesoscale.

Fig. 1 .
Fig. 1.Simulation domain with 2 nested grids.The star indicates Cabauw, the square the location of the radiosonde release, the triangles indicate the scintillometers and the dots the eddy correlation observations.

Fig. 2 .
Fig. 2. Observed and simulated time series at the Cabauw of (a) short wave radiation, (b) potential temperature, (c) wind speed and (d) wind direction at 200 m.

Fig. 3 .
Fig. 3. Sensible heat flux for crops, grasslands and forest at 11 and 12 June 2006 (doy 162 and 163), for locations and full names see Table3.The black dots indicate for Maas en Waal (MeW) and Haarweg (Haa) scintillometer observations and for the other sites eddy correlation observations.Green indicates the RAMS simulation with a low Bowen ratio, and red with a high Bowen ratio (see text for explanation), dark blue is the WRF simulation and light blue is the ECMWF forecast simulation.

Fig. 4 .
Fig. 4. Simulated and observed potential temperature (a) and humidity (b) profiles at De Bilt, and CO 2 mixing ratio (c) profiles at Cabauw, 11 June 2006, at 12:00 (left) and 24:00 (right).Black indicates the radiosonde (a), (b) or tall tower (c) observations, red the simulation with a high Bowen ratio, green the simulation with a low Bowen ratio.Blue indicates the simulated temperature and humidity with a low Bowen ratio, but with the Mellor Yamada instead of MRF turbulence scheme.

Fig. 5 .
Fig. 5. CO 2 respiration (a) and assimilation (b) fluxes for 2 grass (Ca1 and Hor) and 2 crops (Lut and Lon) sites.In black the observations.The red band is the range of CO 2 fluxes simulated with a spread of biosphere model parameters (Table2).Yellow shows the simulation that best fits the CO 2 mixing ratio observations.

Fig. 7 .
Fig. 7. Contribution of the background CO 2 to the mixing ratio at Cabauw 200 m (a) and the cumulative contribution of the CO 2 fluxes in the domain (b).The variation in the CO 2 mixing ratio is mainly determined by the fossil fuel, respiration and assimilation fluxes, where atmospheric signal from assimilation by crops dominates over grass assimilation.

FFFig. 8 .
Fig. 8. Average contribution of the different tracers to the diurnal CO 2 mixing ratio (a) and (b) and its variance (c) at 200 m at Cabauw.(a)shows the influence of the assimilation (assim), respiration (resp), sea and fossil fuel (FF) fluxes.In (b) the assimilation flux influence is separated by vegetation type: urban vegetation, broadleaf forest (Blf), needle leaf forest (Nlf), crops and grass.In (c) the variability is shown which is due to variations in the biospheric (bio), fossil fuel (FF), and sea fluxes, and in the background mixing ratio (BG).

Table 1 .
Specification of the simulation settings.

Table 2 .
CO 2 flux parameter settings based on

Table 3 .
Observation specifications.EC is eddy correlation measurements, Sc is scintillometer, RS is radiosonde and TT is tall tower.

Table 5 .
Overview of the sensitivity tests.Results are given for 11 June 12:00.Simulation settings: Pre-sim is the time simulated before this moment, SM is the soil moisture content, rs min is the minimal stomatal resistance of the low vegetation, veg.frac. is the vegetation fraction.Results: Bowen ratio (β), sensible (H ) and latent (LE) heat flux, all for Cabauw; and observed potential temperature (T ) for De Bilt at 500 m height (radiosonde observation), and its difference (T sim −T obs ) with the simulated value.Sim.ID identifies the two simulations further used in this study (for more information see the text).The bold numbers indicate the change compared to the preceding simulation. ), , for example at 11, 12 and 19; doy 162, 163 and 170.Another important source of www.biogeosciences.net/6/2265/2009/Biogeosciences, 6, 2265-2280, 2009