Mesoscale simulations of atmospheric CO2 variations using a high‐resolution model system with process‐based CO2 fluxes

A new coupled high‐resolution biosphere–atmosphere model (TerrSysMP‐CO2) is applied to simulate mesoscale and diurnal variations of atmospheric CO2 mixing ratios. The model is characterized by process‐based parametrization calculating atmospheric dynamics and biogenic processes considering the prognostically varying CO2 content at the surface. An advanced parametrization of soil respiration is used distinguishing between heterotrophic and autotrophic respiration and explicitly considering the effect of varying soil moisture. In addition to biogenic CO2 fluxes, high‐resolution anthropogenic emissions are included in the simulations.


Introduction
The relevance of increasing atmospheric CO 2 contents for the estimation of future global warming is well known. Thus, an accurate quantification of the CO 2 budget and the distribution of CO 2 sources/sinks is crucial for precise climate simulations both at global and regional scales. To achieve a realistic spatio-temporal distribution of CO 2 fluxes, two contrasting approaches are used. On the one hand, in the 'top-down' method terrestrial CO 2 fluxes are derived from a combination of observed CO 2 variabilities and (global) atmospheric transport models, i.e. inverse modelling (e.g. Gerbig et al., 2009;Peters et al., 2010). On the other hand, in the 'bottom-up' approach field-scale observations of CO 2 exchange, e.g. eddy-covariance (EC) measurements are upscaled to derive the net CO 2 exchange of ecosystems at larger scales. However, due to several uncertainties, a large mismatch exists between local flux measurements and fluxes derived from inverse modelling (e.g. Dolman et al., 2006;Ter Maat et al., 2010). In particular, the representativeness of CO 2 fluxes at field scale is not necessarily guaranteed, and mesoscale atmospheric circulations influencing measured CO 2 mixing ratios are not resolved with global models used for inverse modelling (van der Molen and Pérez-Landa et al., 2007). To bridge the gap between local and global scale, this study deals with the spatio-temporal variability of CO 2 fluxes and atmospheric CO 2 mixing ratios at the mesoscale.
In the planetary boundary layer (PBL), high horizontal gradients were observed in the mesoscale CO 2 distribution in flight campaigns (e.g. Dolman et al., 2006). Therefore, mesoscale models with high spatial resolutions are necessary to properly simulate these gradients caused by variable atmospheric flow as well as by heterogeneous vegetation Smallman et al., 2013;Oney et al., 2015). The question arises as to the main controlling factors (e.g. land use heterogeneity, synoptic flow, mesoscale flow, anthropogenic emissions) explaining CO 2 gradients at the regional and mesoscale where limited observations are available (Nicholls et al., 2004;Ter Maat et al., 2010). This knowledge is required to improve the upscaling of local flux measurements to mesoscale distributions of CO 2 fluxes.
In the last decade, CO 2 mixing ratios and fluxes have been investigated with several mesoscale biosphere-atmosphere models. Nicholls et al. (2004) have analyzed CO 2 variations in the Great Lakes region (USA) and have identified diurnal, local and regional-scale variability caused by varying turbulent mixing and orographically induced atmospheric flows (e.g. katabatic winds). Moreover, Ahmadov et al. (2007) and Sarrat et al. (2009) have compared the results of two model systems with in situ (CO 2 tower, EC stations) and remote-sensing (aircraft) observations in a rather flat region in southwestern France. Both models indicate horizontal patterns in the atmospheric CO 2 distribution resulting from different dominant land use and mesoscale circulations (land-sea breezes) in that particular region. In an intercomparison study, five biosphere-atmosphere models with different levels of complexity have been compared to estimate the uncertainty of mesoscale simulations caused by different parametrization . Moreover, Tolk et al. (2009), Ter Maat et al. (2010) and Smallman et al. (2013) have studied CO 2 fluxes and atmospheric CO 2 patterns over the Netherlands and Scotland by coupling regional atmospheric models to more advanced biosphere models.
Some of these models describe the biosphere with highly simplified diagnostic relations to satellite indices  or with idealized fluxes . Other studies use process-based models but with relatively coarse grid spacings ranging from 8 km (Sarrat et al., 2009) to 4 km (Tolk et al., 2009;Ter Maat et al., 2010) which may be not appropriate for the application in regions with complex terrain. Another limitation of most studies is a coarse resolution of anthropogenic emissions, e.g. 1 • Ter Maat et al., 2010) or 10 km Tolk et al., 2009;Sarrat et al., 2009).
In contrast, in this study we developed and applied a new coupled atmosphere-biosphere model using process-based approaches for biosphere processes. Innovative in this study is that for soil respiration we have included RothC , a commonly applied carbon turnover model (e.g. Jones et al., 2005), in the model instead of simple temperature relationships used in most former studies. Moreover, the fully prognostic spatio-temporal varying atmospheric CO 2 content is used to calculate photosynthesis and transpiration. As additional CO 2 sources we integrate very-high-resolution (1 km) time variable anthropogenic emissions into the simulations. Combined with fine grid spacings (atmospheric transport: 1.1 km, biosphere: 500 m), this model captures most relevant processes to perform precise simulations in regions with complex terrain and densely populated areas.
This study has the following objectives. The new model system is used to simulate mesoscale and diurnal variations of atmospheric CO 2 mixing ratios. For the high-resolution simulations a heterogeneous region was selected to analyze the impact of different processes (e.g. atmospheric flow, complex terrain, land use) on atmospheric CO 2 patterns and we hypothesize that mesoscale circulations induced by the topography control the spatial CO 2 distribution. The results provide valuable information for regional-scale inverse modelling using coarser grid resolutions. Moreover, by comparing simulated vertical CO 2 profiles with observations at several heights on a meteorological tower, the influence of varying intensities of turbulent mixing on the vertical CO 2 distribution is studied for different weather conditions. Finally, we show the effect of using two different plant crop physiological parameter sets on canopy fluxes and the resulting meteorological conditions in the PBL.
In section 2, the modelling system is described as well as the selected domain, lateral boundary conditions and observations used for verification. In section 3, the model results are compared with meteorological and CO 2 mixing ratio measurements and EC fluxes with a special focus on the performance in different turbulent situations. In section 4, the model is applied for the analysis of heterogeneous canopy fluxes and mesoscale spatiotemporal variations of atmospheric CO 2 mixing ratios. The results are summarized and discussed in section 5.

Numerical model and observations
2.1. Numerical modelling system: TerrSysMP-CO 2 The modelling system is based on the Terrestrial Systems Modeling Platform (TerrSysMP; Shrestha et al., 2014) consisting of the regional numerical weather prediction (NWP) model COSMO (Consortium for Small-scale Modelling, version 4.21), the Community Land Model (CLM, version 3.5) and the hydrological model ParFlow. Since in this study only the atmospheric model (COSMO) is coupled to the terrestrial model (CLM), ParFlow is not further considered. Significant improvements at the regional climate scale have been shown by coupling the climate version of COSMO to the advanced land surface model CLM (e.g. Davin et al., 2011;Davin and Seneviratne, 2012). In contrast, TerrSysMP couples the NWP version of COSMO to CLM. Only the most important features of TerrSysMP are described, but a more detailed description of the model components and their coupling can be found in Shrestha et al. (2014). The components of TerrSysMP are coupled via the coupler OASIS3 (Valcke, 2013) using the multiple-executable approach. This means that each model has its own executable, allowing different model time steps and grid spacings of each model component. For the numerical investigation of CO 2 fluxes and spatio-temporal variations of mesoscale CO 2 patterns in the atmosphere, TerrSysMP was extended by a fully prognostic two-way coupling of CO 2 between COSMO and CLM, in the following referred to as TerrSysMP-CO 2 . COSMO prognostically calculates the atmospheric transport of CO 2 (section 2.1.1). The resulting time-variable atmospheric CO 2 content is used for the calculation of CO 2 assimilation and respiration by CLM (sections 2.1.2 and 2.1.3). In turn, the atmospheric CO 2 concentration is influenced by local CO 2 tendencies consisting of both biogenic net CO 2 fluxes provided by CLM, and anthropogenic emissions (section 2.1.4). A schematic of important CO 2 related interactions and submodels of TerrSysMP-CO 2 is shown in Figure 1; the abbreviations are explained below.

Atmospheric transport
The atmospheric state as well as the transport of CO 2 in the atmosphere (advection, turbulence, convection) is calculated by an updated version (version 4.21) of the non-hydrostatic limitedarea NWP model COSMO (Baldauf et al., 2011). In contrast to the original configuration, an advanced version of the positive definite one-dimensional flux-form advection scheme of Bott (1989) is used (Schneider and Bott, 2014). Vertical turbulent mixing is based on the level-2.5 scheme of Mellor and Yamada (1982) and adapted to the COSMO model (Raschendorfer, 2001). For high-resolution simulations (1.1 km), deep cumulus convection is assumed to be explicitly resolved. For the nesting simulations, the hybrid mass-flux convection scheme HYMACS (Kuell and Bott, 2008;Kuell and Bott, 2011) is used for shallow, mid-level and deep convection because at a grid spacing of 2.8 km deep convection is only partly resolved (Uebel and Bott, 2015). The transport of CO 2 mixing ratios was included as a passive fluid tracer and feedbacks on radiative transfer is neglected.

CO 2 assimilation and leaf respiration
CO 2 assimilation by the canopy (A can ) is determined with the biogeophysical parametrization of CLM (Dai et al., 2003;Oleson et al., 2010) coupling a biochemical photosynthesis model with a plant physiologically based stomatal conductance model. Leaf stomatal resistance (r st ), controlling the exchange of water and CO 2 between the leaves and the ambient air (i.e. transpiration and photosynthesis), is calculated similarly to the Ball-Berry approach as described by Collatz et al. (1991). The response of r st (the inverse of the stomatal conductance g st ) to photosynthesis A and environmental conditions is parametrized by In this equation, m is an empirical plant functional type (PFT)dependent scaling factor, c s and e s are the CO 2 partial pressure and water vapour pressure at the leaf surface, respectively, e * i is the saturation vapour pressure inside the leaf at the vegetation temperature T v , p atm is the atmospheric pressure at the surface and g st,min is the minimum stomatal conductance when A is zero.
When CLM is coupled to COSMO, the environmental conditions p atm and e s are prognostic variables provided by the atmospheric model. In TerrSysMP, e s is derived from the atmospheric specific humidity (q v ) at the COSMO surface level ( Figure 1). However, c s is diagnostically calculated from a constant CO 2 mixing ratio in the ambient air. In contrast, the two-way coupling of CO 2 in TerrSysMP-CO 2 additionally includes c s as a prognostic variable which is derived from the specific CO 2 content (q CO 2 ) at the surface level. Thus, in TerrSysMP-CO 2 , r st (Eq. (1)) is not only influenced by the spatio-temporal varying atmospheric humidity, but also by varying atmospheric CO 2 contents. The stomatal control of the H 2 O and CO 2 exchange between the canopy and the atmosphere is then calculated more consistently.
Leaf photosynthesis (A) is calculated with a modified version of the biochemical model of Farquhar et al. (1980), but with an advanced representation of the maximum rate of carboxylation (V c,max ) (Thornton and Zimmermann, 2007). V c,max controls the activity of A, e.g. the capacity utilization limitation (i.e. 0.5V c,max ) as upper limit of A. It is highly variable among different PFTs (Table 1). For more details on the calculation of V c,max and A, the reader is referred to Oleson et al. (2010). Leaf (dark) respiration (R L ) was included in TerrSysMP-CO 2 as R L = 0.015V c,max according to Collatz et al. (1991).
In CLM the upscaling from leaf level to canopy level is realized with the vertical canopy integration scheme of Thornton and Table 1. V c,max (μ mol(CO 2 ) m −2 s −1 ) at the top of the canopy (including N limitation) at T v = 25 • C for a selection of plant functional types ( Figure 3 Zimmermann (2007) which explicitly considers structural and functional characteristics of the canopy. This one vegetation layer scheme uses a 'two-big-leaf' approach distinguishing between sunlit and shaded leaves. Sunlit leaves receive (and absorb) unattenuated direct beam solar radiation (↓S dir ) and diffuse solar radiation (↓S dif ) whereas shaded leaves receive (and absorb) scattered diffuse solar radiation only. V c,max , A and R L are calculated separately using average absorbed photosynthetically active radiation (PAR) for the two leaf types. A can (and R can L ) are then determined by multiplying by the sunlit and shaded leaf area indices L (e.g. A can = A sun L sun + A sha L sha ).

Soil respiration
CLM (version 3.5) includes a carbon-nitrogen (CN) model that can be optionally used for a prognostic calculation of the carbon and nitrogen cycles (Thornton and Rosenbloom, 2005;Oleson et al., 2010). Instead of using this model for soil respiration, we integrated the RothC-26.3 (RothC) model  into CLM. This carbon turnover model simulates the decomposition of organic plant material in the mineral soil. Similar to the carbon turnover approach in the CN model, RothC utilizes a carbon (C) pool concept to calculate the decomposition by soil micro-organisms. Soil organic matter is divided into five C pools: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), humified organic matter (HUM) and inert organic matter (IOM). The decomposition of each C pool (except the fully decomposed IOM) is calculated by an exponential decay: (2) C i represent the C pools DPM, RPM, BIO and HUM and C i,new the C pools after decomposition during the time interval t. f (T soil ) considers the sensitivity on soil temperature (T soil ) and f (h) describes a moisture reduction depending on the soil water pressure head (h). The plant retainment factor c r = 0.6 considers the deceleration of decomposition when plants grow and λ i are the corresponding optimum decomposition rates for pool i (adopted from Coleman and Jenkinson, 2005). Whereas the decomposition rate is independent of soil texture, the partitioning of decomposed material into CO 2 , HUM and BIO depends on the clay content of the soil . Equation (2) is solved for the CLM soil layers l = {1, ..., 7} (0-83 cm depth). The sum of CO 2 produced in this process is released to the atmosphere as heterotrophic respiration (R hetero ). For the sensitivity of microbial activity on T soil (f (T soil )) and on soil moisture (f (h)), the assumptions ofŠimȗnek and  are used. f (T soil ) is described by an exponential dependency formulated as T ref is the reference temperature of RothC (9.25 • C ) where f (T soil ) = 1.0 and we assume Q 10 = 2.1, leading to realistic respiration rates at soil temperatures that usually occur in the temperate zone. Moisture reduction is expressed as: In Eq. (4) h 2 = −1 m is the pressure head for optimal soil respiration, h 3 = −10 5 m is the pressure head when production ceases, and h 1 is equal to the air entry pressure. h 1 depends on the soil texture and is diagnostically predefined in CLM for the given soil type. The values of h 2 and h 3 are adopted from Suarez anď Simȗnek (1993). Eqs. (3) and (4) are solved for each CLM soil level l.
The initial partitioning of total organic carbon (TOC) into the RothC pools is obtained by the pedotransfer functions of Weihermüller et al. (2013), assuming equilibrium dependent on the clay content and the TOC content only. We use TOC profiles measured at more than 500 locations within the model domain, provided by the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen (LANUV, 2012). By averaging all measurements having the same land use, representative TOC depth profiles were determined for each PFT which were linearly interpolated to the soil levels of CLM (Uebel, 2016, gives more details). The main advantages of this method are the elimination of a model spin-up, linked with tremendous computational costs, and the application of measured TOC profiles instead of a simulated equilibrium carbon distribution (final result of the spin-up) as is needed for the CN model.
In forests, above-ground litter and the forest floor (O horizon) are further C pools. The decomposition (R O ) is calculated with the same equations as for the mineral soil (Eqs (2)-(4)). The ground temperature (T g ) is used instead of T soil , clay is set to zero and h of the uppermost soil layer is used. With the relationships of Zhang et al. (2008), the decomposition rates of litter and in the O horizon were estimated as λ litter = 0.4 year −1 and λ O = 0.13 year −1 , respectively. We assume a TOC content of 9.0 t ha −1 for above-ground litter and 19.8 t ha −1 for the O horizon as found in an analysis of the Hessisches Landesamt für Umwelt und Geologie for forests in Hesse, Germany (Wegener, 2008).
Additionally, we developed a rather simple relationship between (below-ground) autotrophic respiration (R auto ) and mean rates of canopy photosynthesis and leaf respiration (A can , R can L ). It is based on the assumption of Ryan (1991) who pointed out that total plant respiration (R can L + R auto ) is about 50% of photosynthesis for many plant species, i.e.
A can and R can L are monthly means, determined by a 1 year simulation of CLM driven with hourly COSMO-DE (2.8 km grid resolution) analyses for the year of interest. A linear interpolation between two months leads to R auto for the simulated day. Optionally, it would also be possible to simulate daily or weekly mean autotrophic respiration with Eq. (5). Moisture variability is taken into account by the moisture reduction factor r auto which is determined by vertically averaging Eq. (4) weighted with the effective root fraction e root (l) of each CLM soil level l: Finally, to obtain the actual autotrophic respiration (R auto ) adjusted with the actual moisture conditions in the root zone, R auto is normalized with the monthly mean reduction r auto and then multiplied by the actual moisture reduction r auto : In summary, compared to similar model systems, an advanced parametrization of soil respiration is included in TerrSysMP-CO 2 distinguishing between heterotrophic and autotrophic respiration. In previous models, total soil respiration has been either diagnostically determined by constant fluxes for different land-use classes  or by simple functions of R soil on temperature, e.g. a linear relation Pillai et al., 2011), a Q 10 function (Sarrat et al., 2009), or the approach of Lloyd and Taylor (1994) (Tolk et al., 2009;Ter Maat et al., 2010). In contrast, we additionally consider the effect of varying soil moisture and apply the novel concept of using measured TOC profiles for the calculation of R hetero .

Anthropogenic emissions
In addition to biogenic CO 2 fluxes, anthropogenic CO 2 emissions strongly influence the atmospheric CO 2 content. We use a preliminary version of the TNO-CAMS CO 2 inventory, being consistent with an updated version of the TNO MACC II emission inventory (based on Kuenen et al., 2014), provided by the Netherlands Organisation for Applied Scientific Research, TNO (H. Denier van der Gon, personal communication, 2013). This high-resolution gridded dataset is based on yearly official national reports of emitted air pollutants. The emissions are split into several Source Nomenclature for Air Pollution (SNAP) sectors; among these power generation, non-industrial combustion, industrial combustion, road transport, other mobile sources and waste treatment are relevant CO 2 sources.
For Central Europe the data were further disaggregated to a resolution of 1.0 km by the Rheinisches Institut für Umweltforschung an der Universität zu Köln (P. Franke, E. Bem and J. Klimpt, personal communication, 2013) by applying highresolution proxy data containing information on the geographical locations of emissions of each SNAP sector. Figure 2 shows the emissions of the inner domain of TerrSysMP-CO 2 (section 2.2.1).
Each SNAP sector has a characteristic temporal variability depending on the season of the year, the day of the week and the time of the day, e.g. non-industrial combustion is distinctly higher in winter than in summer, whereas road transport mainly depends on the time of the day and on the day of the week. To account for these temporal variations, the emission time factors used in the LOTOS-EUROS chemistry-transport model (Schaap et al., 2005) are applied for each of the CO 2 emitting SNAP sectors. This results in hourly anthropogenic emissions which are used as additional source of atmospheric CO 2 to be read in COSMO. To consider the elevated release and additional plume rise of industrial and power plants, these emissions are vertically spread between model levels 43 and 47 (≈ 95-500 m). In general, very detailed spatiotemporal information on anthropogenic emissions and their genesis is available for the numerical simulations.

Model domain and grid spacings
For high-resolution mesoscale simulations of atmospheric CO 2 mixing ratios, the atmospheric component of TerrSysMP-CO 2 (COSMO) uses a rotated spherical grid with a horizontal grid spacing of 0.01 • (≈ 1.1 km). In the vertical a stretched grid of 50 model layers is used, starting with a layer thickness of 20 m at the surface. This ensures a relatively high vertical grid resolution in the PBL with 8 and 16 layers below 500 and 2000 m, respectively. To account for the high spatial heterogeneity of the land surface, the vegetation and soil model component (CLM) has a smaller horizontal grid size of 0.005 • ×0.00775 • (≈ 500 m, regular geographical coordinates, lat/lon). The soil hydrological processes are calculated by CLM using a stretched vertical grid of ten levels down to a depth of 287 cm. The domain covers an area of 167×167 km (150×150 COSMO grid points) including the western part of Germany (DE, Deutschland), parts of Belgium (BE), the Netherlands (NL) and Luxemburg (LUX). This region is very heterogeneous in orography, land use and anthropogenic emissions (Figures 3 and 2). In the central and southern part, the domain is characterized by the hilly terrain of the Eifel with several mountain ridges separated by narrow valleys (Figure 3(a)). Further mountainous regions are east of the rivers Rhine (Rhein) and Moselle (Mosel) (i.e. Bergisches Land, Hunsrück). In contrast, in the northern and northwestern part, the terrain is predominantly flat. Figure 3(b) depicts the land cover in terms of PFTs of CLM based on MODIS land cover data (Shrestha et al., 2014). In mountainous regions and in the Moselle valley the vegetation is dominated by broadleaf and needleleaf forest and by some grassland. The flat terrain is mainly characterized by agriculture and urban areas, the latter having a rather high fraction in this particular domain (13.5%). For all PFTs the predefined monthly values of leaf and stem area index (LAI and SAI) of Shrestha et al. (2014) are used. The classification of the soil is based on the FAO-UNESCO (1975) Soil Map of the World as used in the external parameter set of the operational COSMO model. The predominant soil texture is loam followed by sandy loam (e.g. Rhine valley, Belgian Eifel) and sand (northwestern domain). Small areas are classified as clay-loam.
Along the river Rhine is one of the highest population densities of Central Europe; Bonn, Cologne and Dusseldorf are in this region. Further urban regions are Maastricht, Aachen and Liège ( Figure 3). This is reflected in high anthropogenic emissions (cf. Figure 2) mainly caused by urban traffic, residential and industrial combustion. In contrast, in the rural Eifel region anthropogenic emissions are low. By far the highest CO 2 emissions are produced by three lignite-fired power plants ('+' in Figure 2) in North Rhine-Westphalia (NRW), being responsible for about one third of all fossil fuel burning in NRW, as well as by power plants along the river Maas.

Initialization and lateral boundary conditions
As for all limited-area NWP models, TerrSysMP-CO 2 needs initial conditions (ICs) and lateral boundary conditions (LBCs) for standard atmospheric state variables in COSMO as well as for CO 2 mixing ratios. For this, the inner domain is nested into larger domains having a grid spacing of 0.025 • (≈ 2.8 km) both for COSMO and CLM. The grid spacing is rather fine compared to other studies which use grid sizes of 6 to 48 km for their nesting domains (e.g. Ahmadov et al., 2007;Tolk et al., 2009;Ter Maat et al., 2010, and others). Depending on the predominant direction of the atmospheric flow, three different nesting domains provide LBCs for the inner domain ( Figure 4). The domains are large enough to ensure that at a moderate atmospheric flow (≈ 10 m s −1 in the lower atmosphere), in a 24 h forecast the PBL of the inner domain is not influenced by lateral boundary effects of the outer domain. The nesting simulations are driven by hourly COSMO-EU model analysis data (grid size: 0.0625 • ) provided by the DWD and by hourly anthropogenic emissions. Similar to Ter Maat et al. (2010), atmospheric CO 2 is initialized as horizontally (and vertically) homogeneous mixing ratios adjusted to the actual atmospheric mean of the considered day being derived from the mean seasonal CO 2 variation, measured in 102.5 m height on a tower in Jülich (section 2.3).
For the case-study analyzed in section 4, the nesting simulation is started with a lead time of 24 h to the high-resolution run in order to provide a reasonable heterogeneous three-dimensional CO 2 distribution as IC for the inner domain. The nesting simulation is started at 1800 UTC because at this time of the day the horizontally and vertically homogeneous CO 2 distribution is a good approximation of the real CO 2 distribution due to a deep and well-mixed PBL (Tolk et al., 2009). The simulation of the period investigated in section 3 is initialized with the same method. However, without nudging or other data assimilation methods (which inhibit a free model prediction) this 8-days period has to be split into several short-time simulations. For each of the consecutive simulations, the simulated atmospheric CO 2 fields of the previous runs are used as IC of atmospheric CO 2 mixing ratios. The nesting runs are restarted every 24 h to ensure that the initial CO 2 distribution resulting from the predicted meteorological situation does not deviate too strongly from the meteorological fields of the COSMO-EU analyses driving the model.
For the initialization of CLM, a multi-year model spin-up is required. For that, CLM was run stand-alone for 7 years, driven by hourly COSMO-DE analyses, until a dynamical equilibrium was reached. This spin-up provides realistic soil moisture and temperature profiles needed for the calculation of soil water dynamics and heat transfer in the soil as well as all other state variables needed for CLM.

Observations
The model performance is analyzed by comparing the simulation results with EC measurements of latent (LH) and sensible (SH) heat fluxes and CO 2 fluxes. The EC stations are mapped in Figure 3(a) (×): Merzenhausen (ME, 93 m asl, agriculture), Selhausen (SE, 105 m, agriculture), Niederzier (NI, 102 m, grassland), Rollesbroich (RO, 515 m, grassland), Wüstebach (WU, 610 m, spruce forest). All EC stations were equipped with the LI-7500 Open Path Gas Analyzer (LI-COR Biosciences, Lincoln, USA) measuring high-frequency CO 2 and H 2 O mixing ratios as well as with the CSAT3 3D sonic anemometer (Campbell Scientic, Logan, USA). For the verification of TerrSysMP-CO 2 , 30 min quality flagged averages are used (Mauder et al., 2013) to filter out unrealistic CO 2 , LH and SH fluxes. Moreover, at the Selhausen site soil respiration was measured with the LI-8100 Automated Soil Gas Flux System (LI-COR).
Additionally, meteorological and CO 2 mixing ratio measurements were made on the 124 m tall tower of the Forschungszentrum Jülich GmbH, referred to below as the Jülich tower (red ✚). The tower is located at the eastern side of a 90 m × 40 m wide clearing surrounded by a small and about 25 m tall broadleaf forest. Most of the surrounding flat terrain is characterized by arable land use. The Jülich tower has seven platforms (12.5, 20.0, 32.5, 52.5, 82.5, 102.5, 120.0 m) equipped with meteorological instruments. Similar to observations at Cabauw (The Netherlands; Tolk et al., 2009) or at Beromünster (Switzerland; Satar et al., 2016) we measured CO 2 mixing ratios at several heights. Three LI-840A CO 2 /H 2 O Gas Analyzers (LI-COR) were installed at 12.5, 32.5 and 52.5 m height. At 102.5 m above ground level (agl), the Fast Greenhouse Gas Analyzer CO 2 CH 4 H 2 O (Los Gatos Research, Mountain View, USA) operationally measures CO 2 mixing ratios. This height is appropriate to precisely estimate CO 2 contents in the well-mixed PBL in the daytime (Haszpra et al., 2015).

Verification of the model with observations
In this section, the performance of TerrSysMP-CO 2 is evaluated by comparing simulated meteorological conditions and atmospheric CO 2 variability with observations at the Jülich tower. Also, net ecosystem exchange (NEE) is compared with fluxes from EC stations. In order to capture diverse atmospheric states, we selected a period of eight consecutive days (3-10 June 2014) with different weather conditions. In the surroundings of the Jülich tower, more than 50% of the arable land is cultivated with cereal crops. With only one crop type (PFT = 15) in the current version of TerrSysMP-CO 2 , plant physiological parameters of a cereal crop are most appropriate. Thus, we replaced the default plant physiological parameters of CLM for this PFT by the parameters of Sulis et al. (2015) adapted for extensively fertilized winter wheat in Central Europe.

Meteorological situation
We first describe the synoptic situation of this period. On 3 June, sunny phases alternated with showers in the afternoon. In the morning of 4 June ahead of a cold front, a band of light stratiform rain passed the domain followed by convective rain in the afternoon. On both days the temperatures were rather low and the wind continuously increased to a moderate breeze ( Figure 5). The cold front passage in the late evening of 4 June induced intense showers together with a moderate to strong southwesterly wind which continued on 5 June. To the rear of the front, the weather was temporarily characterized by partial cloudiness and rather cool and dry air masses. In the following days (6-8 June), with a southerly flow very warm and moist air masses were advected causing significant increases of temperature and humidity. The sky was completely clear on 6 June, and fair weather conditions dominated on 7 and 8 June. At noon on 9 June, Jülich was influenced by clouds of a supercell storm and in the evening a very severe squall line passed NRW with damaging wind gusts. Both convective systems are not captured by the model, probably because the ICs underestimated the convective activity. Instead, at night strong showers were simulated in the western part of the domain. On 10 June, again fair weather conditions occurred with some light showers in the evening.
In order to understand the behaviour of simulated and observed CO 2 mixing ratios at the Jülich tower, we compared simulated meteorological variables with the corresponding measurements on the tower ( Figure 5, Table 2).
The 2 m temperatures were simulated realistically (RMSE = 1.95 • C , R 2 = 0.88). The underestimated daytime temperatures on 8 June were caused by simulated advection of cool air masses from cloudy regions in the western part of the domain. On several nights the 2 m temperatures were overestimated, explaining the small warm bias. This is of particular  importance in the following analysis (section 3.3). The simulated temperature values at 120 m agl were very similar to the observations (RMSE = 1.87 • C ) but with a small cold bias. The radiation budget was simulated fairly well (RMSE = 108 W m −2 ; bias = 25.4 W m −2 ; R 2 = 0.70). Even the radiation attenuation by cloudiness on 4 June was captured by the model. The only significant deviations from the measured radiation were underestimated cumulus clouds on 5 and 10 June.
The wind speed at 120 m agl agreed rather well with the observations (RMSE = 1.43 m s −1 ; R 2 = 0.58). Even brief effects (e.g. windless conditions in the evening of 7 June) are simulated. In order to evaluate whether the synoptic flow was simulated realistically, a 6 h moving average was applied in the statistical analysis to filter out the strong fluctuations of the wind speed. The small bias of 0.5 m s −1 can be explained by a reduction of the observed wind speed caused by the forest surrounding the tower.
The specific humidity was rather well predicted (RMSE = 1.31 g kg −1 ; bias = −0.77 g kg −1 ; R 2 = 0.79), but with a small negative bias, especially on 8 and 9 June. This is possibly one of the reasons for the missing convective activity on 9 June.
Finally, it should be mentioned that this analysis is strongly affected by the missing squall line which caused significant deviations from the observations on 10 June. Almost all statistical scores are improved by excluding the last 24 h from the calculations. In particular, solar radiation is improved (RMSE = 99 W m −2 ; bias = 14.4 W m −2 ; R 2 = 0.76) but also the RMSE of wind speed (1.11 m s −1 ) and temperature (2 m: 1.80 • C ; 120 m: 1.62 • C ) are distinctly lower for the first 7 days of the simulation. The rather low correlation of the wind speed (Table 2) results from the missing heavy gust front of the squall line whereas without the last 24 h the correlation is significantly better (R 2 = 0.73).
In general, the meteorological conditions at the Jülich tower are well represented by TerrSysMP-CO 2 . An additional comparison with observations at DWD operational synoptic stations in other regions within the domain confirms that TerrSysMP-CO 2 simulates the synoptic situation realistically with deviations similar to the Jülich tower (not shown).

Verification of simulated NEE
The simulated NEE of different PFTs is compared with EC fluxes (Figure 6). For needleleaf trees, NEE is compared with fluxes above a spruce forest near Wüstebach. The simulated and observed CO 2 fluxes show a good agreement (RMSE = 5.63 μ mol(CO 2 ) m −2 s −1 ; bias 1.05 μ mol(CO 2 ) m −2 s −1 ). This holds in particular for 3 June and for 7-10 June. Hence, the plant physiological parameters seem to be appropriate for a typical needleleaf forest in Central Europe. Only on 5 June the measured CO 2 fluxes tended to be more negative than in the simulation, but high wind speeds on this day led to highly fluctuating (and possibly erroneous) EC fluxes. The small positive bias mainly results from too low measured nocturnal CO 2 release rates due to accumulation below the measurement height (40 m).
Similarly, the simulated daytime NEE is in accordance to observed CO 2 fluxes at the meadow near Rollesbroich (RMSE = 5.16 μ mol(CO 2 ) m −2 s −1 ; bias 0.95 μ mol(CO 2 ) m −2 s −1 ) which has a measured grass height of about 25 cm. On most days, the measured CO 2 uptake is in the range of the simulated NEE or slightly stronger. Nighttime respiration cannot be analyzed due to a lack of reliable EC fluxes.
In 2014, the Merzenhausen field was cultivated with winter barley. In contrast to forest and grassland, the measured NEE is not captured by the model (RMSE = 12.37 μ mol(CO 2 ) m −2 s −1 ; bias −5.52 μ mol(CO 2 ) m −2 s −1 ), simulating a much stronger CO 2 uptake in the daytime than observed. A closer look indicates that the measured CO 2 uptake decreased to the end of the period. The deviations have two reasons. First, winter barley has a phenological cycle earlier by several weeks than the mean cycle of (cereal) crops assumed by the model. Second, in this particular year, unusually high temperatures in spring induced the timing of phenological stages earlier than usual by about 10 days. Both effects cause a shift in the phenological stage of the Merzenhausen field compared to TerrSysMP-CO 2 using a 10 year mean for the seasonal course of LAIs for all crops in this region. Hence, in contrast to the model, barley fields in the flat terrain around Jülich already approached senescence, explaining both the observed decreasing daytime NEE and the overestimation of simulated NEE. A comparison of simulated NEE with a winter wheat (or rye) field, which is less senescent than barley, would probably lead to a much better agreement ( Figure 11). Unfortunately, measurements for winter wheat are not available in 2014. Due to the earlier senescence of all cereal crops in 2014, the simulated daytime NEE seems to be slightly overestimated. Considering the high fraction of arable land in the flat terrain, this leads to lower CO 2 contents within the PBL than in reality (see next section).
The Selhausen site was cropped with sugar beet. Even utilizing cereal crop parameters, NEE agrees better than at the barley field (RMSE = 8.66 μ mol(CO 2 ) m −2 s −1 ; bias −0.70 μ mol(CO 2 ) m −2 s −1 ). The observed CO 2 uptake increased from 3 to 10 June due to rapidly growing plants. Between 27 May and 11 June, the vegetation height increased from 20 to 35 cm. On 8-10 June the simulated fluxes coincide well with the EC fluxes, whereas especially on 6 June the simulation overestimates the CO 2 uptake in the afternoon. The wilting of sugar beet is more sensitive to low soil moisture and low atmospheric humidity than is estimated with the cereal crop parameters (i.e. higher r st than simulated), which explains the deviations on this day.
The model efficiency, expressed by R 2 , is reasonable for the Wüstebach forest and for grassland (Rollesbroich) but seems to be weak for agriculture (not shown). In particular, the barley field (Merzenhausen) was not representative for cereal crops within the domain due to being on average several weeks later in senescence. Moreover, the statistical analysis is based on strongly fluctuating half-hourly EC fluxes (Figure 6), which explains the lower R 2 for canopy fluxes than for atmospheric state variables because these short-time fluctuations cannot be simulated. For a reliable estimation of the model performance, a statistical analysis of mean diurnal NEE cycles over a much longer time period and for several representative plants would be necessary.

CO 2 mixing ratios and PBL conditions
Diurnal CO 2 variations are height-dependent and strongly differ from day to day. Figure 7 depicts simulated and measured CO 2 mixing ratios at the Jülich tower for 3-10 June 2014. Measured CO 2 concentrations showed high diurnal variations on 3, 4, 7 and 8 June, but lower variations on the other days (Figure 7(a)). During the night, the highest CO 2 contents occurred at 12.5 m agl whereas the nocturnal increase became  significantly lower at higher altitudes. In the daytime, CO 2 mixing ratios at all measurement heights were similar (due to a well-mixed PBL). At night on 5, 9 and 10 June, similar CO 2 contents were measured at all altitudes but with slightly stronger fluctuations at 12.5 m on 9 June. Two days ('day 1' and 'day 2') with contrasting behaviour at night are analyzed in more detail below.
The simulated CO 2 contents at the lower four model levels and at 982 m are shown in Figure 7(b). TerrSysMP-CO 2 captures the decreasing CO 2 amplitudes with height, however the vertical CO 2 gradients partially differ from the measurements. As in the observations, for 5 June the daily amplitude is lowest and almost not height-dependent. The simulated boundary layer height (BLH), calculated with the bulk Richardson number approach (Seibert et al., 2000) using a critical value of 0.25, indicates that a shallow thermally stable PBL does not evolve during this night ( Figure 8). Instead, strong winds to the rear of a cold front ( Figure 5) cause efficient and deep turbulent mixing. On 5, 9 and 10 June the daily amplitudes are simulated faily well whereas on 3, 4 and 6-8 June the nocturnal CO 2 increase near the surface is underestimated by about 1/3 (Figure 7(c)). On these nights, shallow nocturnal PBLs formed (BLH < 100 m) characterized by weak winds and a very stable thermal stratification. The latter is partially underestimated by the model, which explains the deviations (see below). The only exception is the night of 10 June (after the missing squall line), showing again a stable thermal stratification (low BLH) resulting in overestimated vertical CO 2 gradients. With increasing BLH every morning vertical gradients of simulated and observed CO 2 mixing ratios diminish. The CO 2 contents about 1 km agl are decoupled from the diurnal variations below, except for some short periods in the afternoon when the BLH exceeds this height level. Figure 7(c) compares the simulated and observed CO 2 mixing ratios near the surface. The 12.5 m platform of the Jülich tower can be considered to be within the forest canopy. The abovedescribed underestimated CO 2 mixing ratios on some nights are apparent. However, except for 7 and 9 June, the simulated daytime CO 2 contents agree well with the observations. Thus, both the model bias of −3.80 ppmv and the RMSE (8.67 ppmv), calculated for 0800-1930 UTC (i.e. convective PBL in the daytime), are considerably lower than the bias of −8.79 ppmv and the RMSE (17.5 ppmv) over the complete simulation period.
With increasing altitude, the differences decrease between simulated and observed nocturnal CO 2 contents. Similarly, in Sarrat et al. (2009) CO 2 contents at 60 m agl have been simulated better than at 20 m height. Compared to the surface level, at 121 m height the simulated CO 2 mixing ratios as well as the diurnal amplitudes agree better with the corresponding measurements (RMSE = 12.6 ppmv; bias −6.17 ppmv; Figure 7(d)). Relative to lower levels, the diurnal variation was considerably lower explaining the lower R 2 than at 10 m. The only exception is the night of 8 June showing very high CO 2 contents at all levels at the Jülich tower not occurring in the simulations. For about 2 h, at the tower the wind turned to SSW. Hence, advection of high CO 2 mixing ratios originating from a big coal-fired power plant at 10 km distance disturbs the observations. This temporary wind shift is not simulated, i.e. the CO 2 plume is not captured. This short period (8 June, 0000-0530 UTC) slightly distorts the statistical scores which are significantly better if this period is excluded from the statistical analysis (12.5 m: RMSE = 15.1 ppmv, bias −8.69 ppmv; 102.5 m: RMSE = 9.7 ppmv, bias −5.29 ppmv).
Especially in the daytime of 7 and 9 June, for several hours the simulated CO 2 contents are 10 ppmv (7 June) to 20 ppmv (9 June) lower than observed. On 7 June low simulated CO 2 mixing ratios originating over large arable areas are advected by a gentle to moderate southeasterly wind. As shown in the previous section, the simulated net CO 2 uptake of cereal crops is overestimated (e.g. bias of −5.52 μ mol(CO 2 ) m −2 s −1 for barley in Merzenhausen) causing CO 2 contents which are too low downstream of arable areas. Moreover, the BLH is rather low (< 500 m) in these periods. Similarly, on 9 June the negative deviations result from CO 2 advection from the same region but the simulated synoptic conditions (section 3.1) additionally affect the CO 2 contents, at least in the afternoon. Thus, horizontal advection strongly controls CO 2 mixing ratios measured on a tower in the daytime. This is in accordance to Tolk et al. (2009), who have found that strong CO 2 uptake of crops has a significant effect on CO 2 mixing ratios measured on the Cabauw tower which is located at a grassland site. At night, the Jülich tower is not significantly influenced by arable land in the surroundings because, in the stably stratified nocturnal PBL, locally respired CO 2 is more important than advective CO 2 transport (e.g. Oney et al., 2015).
In order to better understand why the diurnal CO 2 variations are well simulated on some days but not on other days, vertical CO 2 profiles of two contrasting days are analyzed in combination with corresponding temperature profiles. 'Day 1' represents a positive example. On 4 June, between 2100 and 2300 UTC a cold front passed Jülich with enhanced wind speeds which continued over the following 18 h ( Figure 5). Thus, free tropospheric CO 2 could be mixed to the surface explaining the low nocturnal CO 2 increase. The formation of a stable nocturnal PBL was inhibited resulting in similar CO 2 concentrations at all vertical levels (Figure 7). In the daytime, simulated and observed CO 2 mixing ratios continuously decrease and are in very good agreement.
At 2200 and 0000 UTC, the atmosphere is dry-adiabatically stratified (Figure 9(a)). In the well-mixed and deep PBL, a vertically homogeneous CO 2 distribution is simulated which is in accordance to the tower observations. At 0300 UTC, the simulated thermal stratification is no longer a dry adiabat. At the Jülich tower almost isothermal conditions occurred above the forest top, with a weak inversion within the forest canopy. However, consistent with the simulations, a pronounced wind shear (12.5 m: 0.7 m s −1 ; 120 m: 8.2 m s −1 ) allowed efficient dynamically driven turbulent mixing, which also explaines the high BLH (Figure 8). At 1200 UTC, typical measured and simulated CO 2 and temperature profiles of a convective PBL can be seen.
In the night to 7 June ('day 2'), the differences in CO 2 increase at different heights of the tower were most pronounced. In contrast to 'day 1', the simulation strongly underestimates the nocturnal CO 2 contents at 12.5 m (and 32.5 m), whereas at the other levels the CO 2 mixing ratios are better in agreement (Figures 7(c) and (d)). The reasons for these deviations are analyzed in Figure 9(b).
At 2200 UTC, below 30 m a strong temperature inversion (4.2 • C increase with height) was observed whereas TerrSysMP-CO 2 simulates only a slight inversion (0.6 • C below 72 m). Thus, between 12.5 and 52.5 m the observed vertical CO 2 gradient was stronger than simulated. However, for this strong mismatch of the thermal stratification, even a stronger deviating vertical CO 2 gradient could be expected. A rather pronounced wind shear at the Jülich tower, not occurring in the simulation, inhibited a stronger near-surface CO 2 accumulation at this time. Until 0000 UTC, the observed temperature inversion (i.e. stable nocturnal PBL) extended to higher altitudes and the wind shear weakened. Hence, limited vertical turbulence caused an intense near-surface CO 2 accumulation as well as a strong vertical gradient between 12.5 and 52.5 m height. In contrast to the observations, a relatively uniform vertical temperature increase of only 2.8 • C is simulated below 250 m. Thus, the nearsurface CO 2 mixing ratios and respective vertical gradients are significantly underestimated. Between 0000 and 0300 UTC the observed inversion further intensified to an exceptional degree (8.5 • C increase below 120 m) and was most pronounced around the forest top (20.0-32.5 m). This explains the strong near-surface CO 2 accumulation as well as the distinctly lower concentrations at 32.5 and 52.5 m height. Again, the simulated PBL is deeper with strongly underestimated temperature inversion and vertical CO 2 gradients. Finally, at 0500 UTC the measured temperature inversion began to weaken and the wind had a distinct maximum at 82.5 m (not shown), both inducing stronger mixing. Thus, simulated and observed CO 2 profiles are in better agreement than in the previous hours.
This points out that, in situations with efficient turbulent mixing and high BLHs (e.g. 'day 1'), CO 2 mixing ratios and the rather homogeneous vertical CO 2 profiles are well simulated (cf. also 9 June). In contrast, the underestimated CO 2 mixing ratios near the surface result from a deficient representation of stable atmospheric stratifications causing a too strong turbulent vertical transport. This is a common problem of turbulence schemes used in NWP models (overestimating nocturnal BLHs combined with too weak inversions below), resulting in underestimated CO 2 accumulations near the surface (also found in Tolk et al., 2009). The overestimated nocturnal 2 m temperatures ( Figure 5) are a further indication of this problem. However, these mismatches are to some extent caused by the forest surrounding the small clearing where the Jülich tower is located. On clear-sky nights, the measured temperature and wind profiles regularly show a decoupling of the air within the canopy from the air above. This can be identified by almost calm conditions below 20 m reducing horizontal CO 2 transport as well as by a sharp temperature inversion between 20 and 32.5 m (e.g. 'day 2', 0000 and 0300 UTC) inhibiting the transport of CO 2 released by the canopy towards higher altitudes. These canopy effects cannot be considered by the current version of TerrSysMP-CO 2 using the roughness-length approach of a single vegetation surface. With consideration of drag effects, the influence of tall vegetation on the wind field could be accounted for (e.g. Dupont et al., 2004;Aumond et al., 2013).

Interaction with mesoscale meteorology
Canopy fluxes and the spatio-temporal variability of atmospheric CO 2 mixing ratios are analyzed for 26 May 2012. On this day the temperatures were moderate (≈ 20-25 • C in the afternoon) and it was completely cloudless; this allowed a direct comparison of the simulated biosphere-atmosphere exchange of heat, water and CO 2 with observations without consideration of effects by clouds. The wind speeds decreased during the night whereas a moderate easterly wind (≈ 5-8 m s −1 at 10 m agl) occurred over the flat terrain between 1100 and 1400 UTC. Moreover, at this time in some areas the atmospheric humidity decreased considerably, leading to detrimental environmental conditions for plants.
Due to the predominant easterly flow over Central Europe the domain 'East' (Figure 4) was selected for the parent model assuming a background CO 2 content of 396 ppmv for IC and LBCs. The numerical simulation of the inner domain was initialized on 25 May 1800 UTC to be run for 30 h.

Canopy fluxes and NEE
The atmospheric CO 2 content in the PBL and the evolution of the convective PBL itself strongly depend on the canopy fluxes of CO 2 , water and energy (Dolman et al., 2006;Sarrat et al., 2009;Tolk et al., 2009). We performed two simulations, one with the default plant physiological parameters of CLM for all PFTs ('clmcrop'). Since the MODIS land cover data, used for the PFT classification (Figure 3(b)), assigns 81% of all arable land to cereal crops, in a second simulation ('cereal') we replaced the default crop parameters (PFT = 15) by the winter wheat parameters of Sulis et al. (2015), as also used in the simulations of section 3. Figures 10(a)-(d) show NEE for both simulations. At night (0000 UTC) the vegetated canopy is a CO 2 source caused by soil and leaf respiration. In the 'clmcrop' simulation ( Figure 10(a)) the fluxes are stronger for forests than for crops and grassland because (additionally to R auto and R hetero ) the decomposition of above-ground litter and within the forest floor (R O ) contribute to the CO 2 flux. The highest respiration rates are simulated for broadleaf forests (8-9 μ mol(CO 2 ) m −2 s −1 ); respiration of needleleaf forests is lower (6.5-7.5 μ mol(CO 2 ) m −2 s −1 ). Due to stronger autotrophic respiration and considerably stronger leaf respiration of crops in the 'cereal' simulation, the contrast between forests and arable areas is less evident (Figure 10(b)).
With increasing photosynthesis in the early morning, the canopy becomes a net CO 2 sink. At noon (1200 UTC), in most regions the net CO 2 uptake (i.e. negative NEE) reaches its maximum. Figure 10(c) depicts NEE of the 'clmcrop' simulation. The strongest CO 2 uptake is simulated for broadleaf forests (11-14 μ mol(CO 2 ) m −2 s −1 ), whereas the corresponding rates for grassland and crops are lower (6-8 μ mol(CO 2 ) m −2 s −1 ). For crops these low values are surprising because measurements and several modelling studies in European regions (e.g. Ahmadov et al., 2007;Sarrat et al., 2009;Tolk et al., 2009) indicate a much higher CO 2 uptake of cereal crops and vegetables. In contrast, in the 'cereal' simulation agricultural areas show by far the highest NEE (33-38 μ mol(CO 2 ) m −2 s −1 ) caused by 2.5 to 3-fold higher photosynthesis rates (Figure 10(d)). Moreover, canopy transpiration of crops is 2.25 to 2.5-fold higher than in the 'clmcrop' simulation, leading to higher latent heat and lower sensible heat fluxes (not shown). The more efficient photosynthesis mainly results from improved parameters controlling RuBisCO enzyme kinematics (Sulis et al., 2015, Table 1). For winter wheat, amongst others, a lower leaf C:N ratio of 14.0 ('clmcrop' 25.0) and a triplication of N in RuBisCO increases V c,max and thus photosynthesis. This is probably related to fertilization which is not accounted for in the 'clmcrop' parameters (Sulis et al., 2015). Finally, a lower m in Eq. (1) allows for a higher photosynthesis rate for the same r st .
Moreover, Figures 10(c) and (d) show rather low NEE rates of needleleaf forest in the central Eifel region (3-6 μ mol(CO 2 ) m −2 s −1 ) being lower than in the late morning. This reduction is correlated with the above-described decrease of atmospheric humidity in this region (simulated dew point temperatures of 1-3 • C ) in combination with increasing temperatures and wind speeds. An analysis of the simulated stomatal resistance (r st ) indicates values of about 600-1200 s m −1 in the entire domain, which is considerably higher than for optimum conditions (r st ≈ 150 s m −1 ). The highest resistances are simulated for needleleaf forests with 1800-2100 s m −1 (not shown). This means that on this particular day, due to limited soil moisture, the atmospheric conditions require a partial closure of leaf stomata at noon to avoid desiccation, to which needleleaf forests are most sensitive, consistent with the findings of Sarrat et al. (2009). Finally, we have to mention that for urban areas R hetero is neglected and a low LAI of 0.6 only allows for weak fluxes of R auto , R can L and A can , i.e. weak NEE rates both at night and in the daytime.
Figures 10(e) and (f) show the variation of heterotrophic respiration in the mineral soil (R hetero ) and in the forest floor (R O ). The lowest rates are simulated at 0700 UTC with 1.6 μ mol(CO 2 ) m −2 s −1 for arable soils and 3.6-4.6 μ mol(CO 2 ) m −2 s −1 in forests with the highest rates in the Rhine valley (Figure 10(e)). At 1500 UTC, heterotrophic respiration reaches its maximum (Figure 10(f)) with slightly higher values for arable soils and 4.4-5.8 μ mol(CO 2 ) m −2 s −1 in forests, depending on topographic height. The rates of needleleaf and broadleaf forests are similar. The lower respiration rates of arable soils result from lower TOC contents in the entire soil column and from a negligible organic matter layer (R O = 0 in TerrSysMP-CO 2 ). The stronger diurnal variation in forests than in agricultural regions can be explained by higher TOC contents in the topsoil where T soil shows both a pronounced diurnal variation and higher values than in the subsoil in the afternoon. Moreover, respiration of litter and respiration in the first few centimetres of soil (O horizon) are strongly temperature dependent and contribute significantly to soil respiration. Despite rather dry conditions prior to the simulation, soil moisture limitation is weak with a maximum reduction of about 10% in the lee of the Eifel. Hence, in the humid climate moisture reduction of soil respiration (Eq. (4)) in loamy soils is rather low, except after heavy rain events causing (super-)saturated soils (not shown). As in section 3.2, the quality of simulated fluxes is evaluated with EC measurements of NEE, LH and SH fluxes at the crop site Merzenhausen (ME, winter wheat in 2012) and at the grassland sites Rollesbroich (RO) and Niederzier (NI) (Figure 11). For this day no data are available for the forest site Wüstebach.
Obviously, in the daytime the 'clmcrop' simulation (dashed lines in Figure 11) cannot capture the strong CO 2 uptake of winter wheat (Merzenhausen). This is consistent with a significant underestimation (overestimation) of LH (SH) fluxes. A possible reason of these deviations may be a too low LAI (2.2 in TerrSysMP-CO 2 ), compared to LAI measurements performed at this site and other cereal crop fields in this region. Moreover, the N limitation factor of 0.61 (i.e. 39% reduction to maximum photosynthesis) used for crops in CLM by default can be adapted for fertilized fields. However, the most important reason for the too weak CO 2 uptake is the inappropriate plant physiological parameter set of PFT = 15 in CLM. Noticeable in all simulated NEE curves of crops in this and other case-studies is a flattening between 0900 and 1700 UTC not occurring in the observations. Without any other restrictions, the crop parameters in the 'clmcrop' simulation lead to a V c,max of 20-25 μ mol(CO 2 ) m −2 s −1 for the simulated afternoon temperatures. Thus, the capacity utilization limitation (section 2.1.2) as upper limit of photosynthesis is already reached for leaf photosynthesis rates of 10-12.5 μ mol(CO 2 ) m −2 s −1 . In combination with respiration, this limitation explains the saturation of NEE at about −8 μ mol(CO 2 ) m −2 s −1 . In contrast, the winter wheat parameters of Sulis et al. (2015) ('cereal' simulation) clearly represent NEE and LH fluxes better (solid lines) due to a much higher V c,max (Table 1). SH fluxes are still slightly overestimated by the model. Finally, chamber measurements made on the previous day (25 May 2012) show that the simulated soil respiration agrees with observations at a further crop field (Selhausen; Figure 12).
CO 2 assimilation of the meadow in Rollesbroich between 0700 and 1600 UTC is underestimated. The saturation again indicates a capacity utilization limitation of simulated grassland at lower photosynthesis rates than for midlatitude meadows, but in contrast to agriculture, grassland clearly plays a minor role for this domain (5% land cover). Soil respiration is well estimated, proven by the agreement between simulated and measured NEE at night. Consistent with the underestimation of NEE, the simulated LH (SH) fluxes are too low (high). The measured daytime CO 2 fluxes at the meadow near Niederzier better coincide with the simulations. Thus, the partitioning of LH and SH is better than at Rollesbroich with only slightly lower (higher) simulated LH (SH) fluxes. One reason for the different measured fluxes of grassland is the lower canopy height (≈ 10 cm) in Niederzier compared to the two meadows influencing the EC fluxes in Rollesbroich (56 and 20 cm). The reduction of simulated NEE between 1100 and 1500 UTC is caused by a partly stomatal closure as response to low simulated dew-point temperatures and limited soil moisture in this region in the afternoon. A similar but weaker reduction can also be seen at Rollesbroich.

Mesoscale atmospheric CO 2 distribution
The relative influence of synoptic and mesoscale transport, diverse land cover and orography on the spatio-temporal variability of near-surface CO 2 mixing ratios is investigated over the rather complex and heterogeneous region described in section 2.2.1.
During the night, it is mainly respiration which causes a continuous increase of atmospheric CO 2 contents due to a low BLH. In the early morning (0400 UTC), the CO 2 mixing ratios are distributed very heterogeneously (Figure 13(a)). High CO 2 contents (420-440 ppmv) occur in the flat terrain as well as in the Moselle valley and some other narrow valleys in the Eifel (e.g. 6.3-6.5 • E, 50.5-50.75 • N). In contrast, on mountain ridges in the Eifel, the CO 2 concentrations increase only slightly (400-410 ppmv) although these areas are dominated by forests with intense respiration rates. The enhanced CO 2 mixing ratios in narrow valleys seem to be resolved only with the high horizontal resolution of TerrSysMP-CO 2 because smallscale mountain-valley effects are responsible for the horizontal CO 2 gradients. The highest CO 2 mixing ratios are simulated in the windward north of the Eifel (435-450 ppmv). In this region, the synoptic-scale easterly wind in the flat terrain is modified, resulting in a confluent flow (persisting over several hours) and a lowering of the BLH (not shown). This emphasises that in this case-study both the synoptic flow and mesoscale changes of the wind field, induced by orographic effects (e.g. windward of the Eifel, valley breeze), are mainly responsible for the CO 2 distribution in the early morning. The heterogeneous patterns of respiration due to different PFTs are less relevant.
With the onset of photosynthesis after sunrise (≈ 0400 UTC), in rural areas the CO 2 mixing ratios rapidly decrease within 2 h although the net CO 2 uptake is rather low in the early morning ( Figure 13(b)). The transition from a shallow nocturnal PBL to a deep convective PBL occurs later in the morning than the change from a positive to a negative NEE, explaining the strong CO 2 decrease. In contrast, road traffic rapidly increases (morning rush hour at 0600-0700 UTC). Windless conditions (< 1 m s −1 ) in the Rhine valley allow for local CO 2 accumulations around the cities of Bonn, Cologne and Dusseldorf. A similar effect can be seen for Liège. Compared to a simulation without anthropogenic emissions, the CO 2 mixing ratios in the flat terrain are about 5-10 ppmv higher with local deviations > 30 ppmv in urban regions along the river Rhine (not shown). Due to the contrasting behavior of rural and urban areas, at 0600 UTC the greatest heterogeneity is simulated with about 395 ppmv in eastern Belgium and northeast of Cologne and more than 440 ppmv in the metropolitan areas. The local effects of fossil fuel emissions in cities have been also shown (e.g. in Pérez-Landa et al., 2007), but with a very simplified representation of anthropogenic emissions in their study.
Between 0700 and 1200 UTC, in most regions CO 2 continuously decreases to 385-395 ppmv (i.e. below the background CO 2 content). In the flat terrain west of the Rhine, at 0900 UTC (Figure 13(c)) downstream of urban areas and big power plants the CO 2 mixing ratios are still about 3-10 ppmv higher (395-410 ppmv) than without fossil fuel emissions; a rather low BLH (400-500 m) also contributes to this. In general, with growing BLH the horizontal CO 2 gradients begin to diminish. Compared to the morning hours, at 1600 UTC (Figure 13 convective PBL. Vertical cross-sections indicate almost constant CO 2 mixing ratios below the BLH (not shown).
In the afternoon, the CO 2 concentrations within the PBL are below the background CO 2 content which remains rather constant in the free atmosphere. Thus, on this clearsky day, the land surface is a sink of atmospheric CO 2 . Table 3 lists the contributions of the surface CO 2 budget averaged over the entire model domain from 25 May 1800 UTC to 26 May 1800 UTC. It indicates a biogenic sink of −4.297 μ mol(CO 2 ) m −2 s −1 and a total sink (including anthropogenic emissions) of −1.385 μ mol(CO 2 ) m −2 s −1 . In the considered region, anthropogenic emissions are an intense CO 2 source (28% of total CO 2 source). The net surface CO 2 sink in the 'cereal' simulation seems to be more realistic than the gain (+1.757 μ mol(CO 2 ) m −2 s −1 ) in the 'clmcrop' simulation. It is consistent with a measured CO 2 decrease at the Jülich tower in spring and summer (not shown) and with the maps of monthly CO 2 fluxes (CarbonTracker Europe, 2014).

Effect of modified canopy fluxes on PBL conditions
As expected, the strong influence of different plant physiological parameters of crops (PFT = 15) on the canopy exchange of CO 2 (i.e. NEE), water and energy (LH and SH fluxes) has an effect on the atmospheric conditions within the PBL. Higher rates of above-ground and below-ground autotrophic respiration (R can L + R auto ) in the 'cereal' simulation during the night cause 5-15 ppmv higher near-surface CO 2 mixing ratios in the flat terrain where agriculture is the dominant land use (not shown). In contrast, in the daytime the stronger CO 2 uptake results in lower CO 2 mixing ratios near the surface. This means that a higher diurnal CO 2 variation is simulated in the 'cereal' simulation than in the 'clmcrop' simulation, combined with a faster CO 2 decrease in the morning. Figure 14(a) depicts the difference of the nearsurface CO 2 mixing ratios of both simulations at 1200 UTC. Above and downstream of agricultural areas in the flat terrain west of the river Rhine, 5-10 ppmv (locally 15 ppmv) lower CO 2 mixing ratios are simulated in the 'cereal' simulation whereas, in the Bergisches Land east of the Rhine, the difference is less because forests dominate in and upstream of this region. Distinctly higher transpiration rates in the 'cereal' simulation increase the humidity in the PBL. Figure 14(b) shows on average 1-2 • C (locally 4 • C ) higher 2 m T dew at 1200 UTC. In contrast, the 2 m temperature is on average 0.5-0.75 • C (locally > 1 • C ) lower than in the 'clmcrop' simulation due to stronger evaporative cooling (Figure 14(c)). Vertical cross-sections indicate similar differences up to the BLH. Thus, the effects of CO 2 , T dew and T, depicted in Figure 14, are representative for the entire convective PBL.

Summary and conclusions
We have described and applied a new coupled biosphere-atmosphere model (TerrSysMP-CO 2 ) for the simulation of mesoscale and diurnal variations of atmospheric CO 2 mixing ratios in a region with heterogeneous land use, densely populated areas and complex terrain. The process-based canopy parametrization uses the prognostically varying atmospheric CO 2 contents at the surface to calculate canopy transpiration and photosynthesis. Moreover, the model includes an advanced representation of soil respiration, distinguishing between heterotrophic and autotrophic respiration and explicitly considering the effect of varying soil moisture. Another novel concept is the use of measured TOC profiles for the calculation of heterotrophic respiration. In the numerical simulations, the strong heterogeneity of the land surface is accounted for by using a fine horizontal grid resolution (500 m) for the land surface component, but also the atmospheric conditions are simulated with a fine grid spacing (≈ 1 km) to capture relevant mesoscale flow patterns. High-resolution anthropogenic emissions, separated into several emission classes, take into account the human-induced influence on the spatio-temporal CO 2 distribution within the PBL.
The model performance has been tested for eight consecutive days with different weather conditions. A comparison of temperature, radiation, wind and humidity with corresponding measurements on a meteorological tower shows high correlations with rather low biases. Thus, the synoptic situation of this period is captured by the model. A significant weakness is the slightly overestimated 2 m temperatures on cloudless nights caused by simulated turbulent mixing which is too strong. NEE (and energy fluxes) of needleleaf forests and grassland are simulated fairly well, but the model performance strongly depends on the grass height. For agriculture the evaluation is more difficult. With only one PFT used for the entire arable area, the simulated NEE at each agricultural grid point is comparable. In contrast, real NEE of arable land in this region is very heterogeneously dependent on the different physiological behaviour and phenological stages of different crops. Thus, the model fails to simulate NEE for a barley field in 2014 having a rather early phenological cycle. However, in the clear-sky case in 2012, the CO 2 uptake of winter wheat in the mature stage is well simulated. This shows that the physiological parameters are appropriate for cereal crops when the phenological stage coincides with the observed field. We conclude that only one PFT representing agriculture is insufficient for a region where very diverse crops (e.g. cereals, sugar beet, maize, vegetables) are cultivated. In a later version of TerrSysMP-CO 2 , this variation will be considered by integrating additional agricultural PFTs. Finally, we want to mention that only a limited time period has been used for verification. For an extensive validation of the performance of the new model, longer periods and a comparison with several synoptic stations would be necessary.
Both CO 2 measurements on a meteorological tower and the simulations show very different nocturnal CO 2 increases, depending on height and on atmospheric conditions. In situations with efficient turbulent mixing (e.g. convective PBL in daytime, frontal passage with strong winds) the rather homogeneous vertical CO 2 profiles are well predicted. However, especially on cloudless nights with weak winds, the simulated near-surface CO 2 accumulation, vertical CO 2 gradients and vertical stability are underestimated on several days. We conclude that pronounced temperature inversions, i.e. very stably stratified nocturnal PBLs, are not sufficiently reproduced by the model, as supported by slightly overestimated 2 m temperatures on these nights. Thus, too efficient vertical CO 2 transport inhibits realistic CO 2 accumulations in the simulations. This is an important message for the meteorological community indicating that parametrization currently used in most NWP models show deficiencies in simulating the transport of CO 2 (or air pollutants) in such situations. However, we have to mention that some of these mismatches are also affected by the forest surrounding the tower. This often results in a decoupling of the air within the forest canopy from the air above; this cannot be addressed with a single vegetation surface in the model. In a future version of TerrSysMP-CO 2 , the dynamical influence of tall vegetation will be accounted for by including drag terms in the dynamic equations.
The simulated exchange of water, heat and CO 2 between the canopy and the atmosphere has been analyzed for two different crop physiological parameter sets. Compared to EC fluxes on a winter wheat field on a clear-sky day, the simulation using the default parameters of CLM significantly underestimates the CO 2 uptake of cereal crops in the daytime. Similarly, canopy transpiration (and latent heat fluxes) are underestimated, resulting in too high sensible heat fluxes. In contrast, the simulation using plant physiological parameters calibrated for winter wheat in Central Europe indicate a much better agreement of NEE, latent and sensible heat fluxes to the observations. Adopted parameters controlling photosynthetic efficiency, being appropriate for fertilized cereal crops, explain this improvement. The modified canopy fluxes strongly influence the simulated atmospheric conditions showing lower CO 2 contents, higher humidity and lower temperature in the entire convective PBL above and downstream of large agricultural areas. This emphasises that realistic canopy fluxes are essential for an accurate simulation of the atmospheric conditions within the PBL.
Soil respiration is simulated realistically by the advanced parametrization in TerrSysMP-CO 2 . Higher respiration rates and stronger temporal variations are simulated for forests than for arable soils. This is caused by higher TOC contents in total and especially in the topsoil (where distinct variations of T soil occur) as well as by respiration of litter and within the forest floor which are highly sensitive to temperature changes.
The mesoscale application of the model shows no clear correlation between simulated CO 2 mixing ratios and the predominant land use. Hence, the convective boundary-layer budget method (Lloyd et al., 2001), deriving NEE by means of CO 2 profiles in the well-mixed PBL at the landscape scale (≈ 1−10 km), cannot be applied in regions with complex terrain. Instead, distinctly higher near-surface CO 2 concentrations are simulated in valleys or to windward of the Eifel than above elevated mountain ridges. Thus, we conclude that in regions with hilly terrain, the CO 2 patterns are strongly influenced by terrain-induced local circulations, especially at night and in the morning. In the afternoon, efficient turbulent mixing and a high BLH significantly reduce horizontal CO 2 gradients within the convective PBL. In a supplementary study, TerrSysMP-CO 2 will be used to investigate in detail the influence of complex terrain on the thermodynamic behaviour of the atmosphere controlling the near-surface CO 2 distribution. Additionally, the appropriate grid resolution necessary to resolve these processes will be investigated. This is a very important question for the use of CO 2 concentrations from continental measurement sites in regional-scale inverse modelling. Regional models using grid sizes of several kilometres may be insufficient to realistically simulate the spatial CO 2 distribution in regions with complex terrain. This may be one reason for the mismatch between local flux measurements and CO 2 fluxes estimated from inverse modelling. Thus, TerrSysMP-CO 2 provides important information to close the gap between these scales. Moreover, significantly higher CO 2 contents are simulated in and around big cities, indicating intense CO 2 sources from fossil fuel emissions from industrial or urban regions. The very-high-resolution emission dataset used in the TerrSysMP-CO 2 simulations will allow us to quantify the anthropogenic effect for different weather situations in a future study.