Edinburgh Research Explorer Simulating the atmospheric CO2 concentration across the heterogeneous landscape of Denmark using a coupled atmosphere–biospheremesoscale model system

. Although coastal regions only amount to 7 % of the global oceans, their contribution to the global oceanic air–sea CO 2 exchange is proportionally larger, with ﬂuxes in some estuaries being similar in magnitude to terrestrial surface ﬂuxes of CO 2 . Across a heterogeneous surface consisting of a coastal marginal sea with estuarine properties and varied land mo-saics, the surface ﬂuxes of CO 2 from both marine areas and

Abstract. Although coastal regions only amount to 7 % of the global oceans, their contribution to the global oceanic air-sea CO 2 exchange is proportionally larger, with fluxes in some estuaries being similar in magnitude to terrestrial surface fluxes of CO 2 .
Across a heterogeneous surface consisting of a coastal marginal sea with estuarine properties and varied land mosaics, the surface fluxes of CO 2 from both marine areas and terrestrial surfaces were investigated in this study together with their impact in atmospheric CO 2 concentrations by the usage of a high-resolution modelling framework. The simulated terrestrial fluxes across the study region of Denmark experienced an east-west gradient corresponding to the distribution of the land cover classification, their biological activity and the urbanised areas. Annually, the Danish terrestrial surface had an uptake of approximately −7000 GgC yr −1 . While the marine fluxes from the North Sea and the Danish inner waters were smaller annually, with about −1800 and 1300 GgC yr −1 , their sizes are comparable to annual terrestrial fluxes from individual land cover classifications in the study region and hence are not negligible. The contribution of terrestrial surfaces fluxes was easily detectable in both simulated and measured concentrations of atmospheric CO 2 at the only tall tower site in the study region. Although, the Luo et al., 2012;Lovenduski and Bonan, 2017). To have the best chance of accurately predicting the future evolution of the carbon cycle, and its implications for our climate, it is important to minimise the uncertainties that exist presently (Carslaw et al., 2018). Enhanced knowledge and a better process understanding in ecological theory and modelling could potentially reduce the model structural uncertainties (Lovenduski and Bonan, 2017), which together with improvements in spatial surface representation could minimise the current uncertainties.
Studying surface exchanges of CO 2 on regional to local scale can be accomplished with mesoscale atmospheric transport models. Their resolution is in the range from 2 to 20 km, and their advantage is their capability to get a better process understanding of both atmospheric and surface exchange mechanisms in order to improve the link between observations and models at all scales, i.e. for both mesoscale, regional and global models (Ahmadov et al., 2007). The higher spatial resolution of mesoscale models allows for a better representation of atmospheric flows and for a more detailed surface description, which is necessary in particular for heterogeneous areas. In previous mesoscale model studies, biosphere models have been coupled to the mesoscale atmospheric models ranging in their complexity, from a simple diagnostic (Sarrat et al., 2007b;Ahmadov et al., 2007Ahmadov et al., , 2009) to mechanistic-process-based biosphere models (Tolk et al., 2009;Ter Maat et al., 2010;Smallman et al., 2014;Uebel et al., 2017). The modelled CO 2 concentrations and surface fluxes from mesoscale model systems compare better with observations than global model systems (Ahmadov et al., 2009). The atmospheric impact on surface processes related to the ecosystem's sensitivity and CO 2 exchange can be examined in greater detail (Tolk et al., 2009), and tall tower footprints can be studied more concisely (Smallman et al., 2014).
Heterogeneity can also be considerable in coastal oceans, and like terrestrial surface fluxes, the high spatio-temporal variability leads to large uncertainties in estimates of coastal air-sea CO 2 fluxes (Cai, 2011;Laruelle et al., 2013). Coastal seas play an important role in the carbon cycle facilitating lateral transport of carbon from land to open oceans, but almost 20 % of the carbon entering estuaries is released into the atmosphere, while 17 % of the carbon inputs to coastal shelves come from atmospheric exchange . The air-sea CO 2 exchange is in general numerically larger for estuaries than shelf seas (Chen et al., 2013;Laruelle et al., 2010;Laruelle et al., 2014) and can, for estuaries, be as large as 1958 gC m −2 yr −1 , while continental shelf seas have fluxes in the range of −154 to 180 gC m −2 yr −1 (Chen et al., 2013). The large spatial and temporal heterogeneity of the coastal ocean adds to the large uncertainty related to the annual estimates of the air-sea CO 2 exchange . The observed high spatial and temporal variability (Kuss et al., 2006;Leinweber et al., 2009;Vandemark et al., 2011;Norman et al., 2013;Mørk et al., 2016) are not always included in marine models (Omstedt et al., 2009;Gypens et al., 2011;Kuznetsov and Neumann, 2013;Gustafsson et al., 2015;Valsala and Murtugudde, 2015), let alone being taken into account in atmospheric mesoscale systems simulating CO 2 (Sarrat et al., 2007a;Geels et al., 2007;Law et al., 2008;Tolk et al., 2009;Broquet et al., 2011;Kretschmer et al., 2014). But a recent study found that short-term variability in the partial pressure of surface water CO 2 (pCO 2 ) can substantially affect simulated annual fluxes in certain coastal areas -in their case the Baltic Sea (Lansø et al., 2017). Moreover, direct eddy covariance (EC) measurements in the Baltic Sea have shown that upwelling events with rapid changes in pCO 2 greatly increase the air-sea CO 2 exchange (Kuss et al., 2006;Rutgersson et al., 2009;Norman et al., 2013).
In this study we aim to simulate surface exchanges of CO 2 at a high spatio-temporal resolution across a region neighbouring the Baltic Sea, alternating between land and coastal sea together with mesoscale atmospheric transport. A newly developed mesoscale modelling system is used to assess and understand the dynamics and relative importance of the marine and terrestrial CO 2 fluxes. The Danish Eulerian Hemispheric Model (DEHM) forms the basis of the framework, while the mechanistic biospheric soil-plantatmosphere model (SPA) is dynamically coupled to the atmospheric model. Both models are driven by methodological data from the Weather Research and Forecasting model (WRF). The air-sea CO 2 exchange is simulated at a high temporal resolution with the best applicable surface fields of pCO 2 for the Danish marine areas. Tall tower observations are used to evaluate the simulated atmospheric concentrations of CO 2 .
Section 2 is dedicated to describing the study region, which is followed by a detailed description and evaluation of the atmospheric and biospheric model components of the developed model system in Sect. 3. Section 4 contains results, while the Discussion and Conclusions follow in Sects. 5 and 6.

Study area
The study area is comprised of Denmark, a country that is characterised by a mainland (Jutland) and many smaller islands, all containing a varied land mosaic of urban, forest and agricultural areas. With more than 7300 km of coastline encircling approximately 43 000 km 2 of land, many landsea borders are found throughout the country, adding to the complexity (Fig. 1). Denmark is positioned in the transition zone between the Baltic Sea, a marginal coastal sea with low salinity, and the North Sea, a continental shelf sea. Bordering the Baltic Sea, the Danish inner waters are rich in nutrients and organic material (Kuliński and Pempkowiak, 2011). This fosters high biological activity in spring and summer, lowering surface water pCO 2 and allowing for uptake of at- Figure 1. The study region of Denmark (land masses in grey), with the location of the five EC sites shown in black and the Risø campus tall tower site indicated in red. mospheric CO 2 . In winter, mineralisation increases pCO 2 (Wesslander et al., 2010), and outgassing of CO 2 to the atmosphere takes place. The North Sea is a persistent sink of atmospheric CO 2 , where a continental shelf-sea pump efficiently removes pCO 2 from the surface water and transports it to the North Atlantic Ocean (Thomas et al., 2004). This study uses the definition of the Danish exclusive economic zone (EEZ) to estimate the Danish air-sea CO 2 exchange, as the coastal state (in this case Denmark) has the right to explore, exploit and manage all resources found within its EEZ (United Nations Chapter XXI, 1984). The Danish EEZ is approximately 105 000 km 2 .
A tiling approach with the seven most common biospheric land cover classifications were selected for the current study, including deciduous forest (3348 km 2 ), evergreen forest (1870 km 2 ), winter barley (1211 km 2 ), winter wheat and other winter crops (9269 km 2 ), spring barley and other spring crops (5368 km 2 ), grassland (6924 km 2 ), and agricultural other (3909 km 2 ), but excluding urbanised areas. The agricultural other land cover classification includes all agriculture that does not classify as cereals, and as such contains root crops, fruits, corn, hedgerows and "undefined" agriculture. This classification corresponds to the actual crop distribution of 2011 (Jepsen and Levin, 2013).

Observations of atmospheric CO 2
One tall tower is found within the study area on the eastern inner shore of Roskilde Fjord. Here atmospheric continuous measurements have been conducted at the Risø campus tower site (55 • 42 N, 12 • 05 E) during 2013 and 2014. The tower is located on small hill 6.5 m above sea level (Sogachev and Dellwik, 2017). Roskilde Fjord is a narrow micro-tidal estuary that is 40 km long with a surface area of 123 km 2 , has a mean depth of 3 m and is found in the sec-tor 200-360 • relative to the Risø campus tower (Mørk et al., 2016). The city of Roskilde, with around 50 000 inhabitants, is positioned approximately 5 km to the southwest of the site, while Copenhagen lies 20 km to the east.
The tall tower continuous measurements of atmospheric CO 2 concentrations at the Risø campus tower were carried out by the use of a Picarro G1301 placed in a heated building. The inlet was 118 m above the surface, and the tube flow rate was 5 slpm. The Picarro was new and calibrated by the factory. The calibration was checked by a standard gas of 1000 ppm CO 2 in atmospheric air (Air Liquide). During the measurement period from the middle of 2013 to the end of 2014, the instrument showed no other drift than the general increase in the global atmospheric concentration.

Model set-up
The model framework used in the present study consists of two models: DEHM and SPA. A coupling between the two was made for the innermost nest of DEHM in order to simulate the exchange of CO 2 between the atmosphere and terrestrial biosphere at a high temporal (1 h) and spatial resolution (5.6 km × 5.6 km) for the area of Denmark.

DEHM
DEHM is an atmospheric chemical transport model covering the Northern Hemisphere, with a polar stereographic projection true at 60 • N. Originally developed to study sulfur and sulfate (Christensen, 1997), the DEHM model now contains 58 chemical species and nine groups of particular matter (Brandt et al., 2012). This adaptable model has been used to study atmospheric mercury , persistent organic pollutants (Hansen et al., 2004), biogenic volatile organic compounds' influence on air quality (Zare et al., 2014), emission and transport of pollen (Skjøth et al., 2007), ammonia and nitrogen deposition (Geels et al., 2012a, b), and atmospheric CO 2 Geels et al., 2004;Geels et al., 2007;Lansø et al., 2015). The CO 2 version of DEHM was used in the present study. DEHM has 29 vertical levels distributed from the surface to the 100 hPa surface with approximately 10 levels in the boundary layer. Horizontally, DEHM has 96 × 96 grid points, which through its nesting capabilities increase in resolution from 150 km × 150 km in the main domain to 50 km × 50 km, 16.7 km × 16.7 km and 5.6 km × 5.6 km in the three nests. The two-way nesting replaces the concentrations in the coarser grids by the values from the finer grids.

Surface fluxes in DEHM
Anthropogenic emissions of CO 2 , wildfire emissions and optimised biospheric fluxes from the NOAA ESRL Carbon-Tracker system (Peters et al., 2007)  Hourly anthropogenic emissions on a 10 km × 10 km grid from the Institute of Energy Economics and Rational Energy Use (IER, Pregger et al., 2007) were applied for Europe instead of emissions from CT2015. Furthermore, these are, for the area of Denmark, substituted by hourly anthropogenic emissions with an even higher spatial resolution of 1 km × 1 km (Plejdrup and Gyldenkaerne, 2011). As the European and Danish emission inventories were from 2005 and 2011, respectively, the emissions were scaled to annual national total CO 2 emissions of fossil fuel and cement production conducted by EDGAR (Olivier et al., 2014) in order to include the yearly variability in national anthropogenic CO 2 emissions.

Air-sea CO 2 exchange
The exchange of CO 2 between the atmosphere and the ocean, F CO 2 , was calculated by F CO 2 = Kk 660 pCO 2 , where K is solubility of CO 2 calculated as in Weiss (1974), k 660 is the transfer velocity of CO 2 normalised to a Schmidt number of 660 at 20 • C, and pCO 2 is the difference in partial pressure of CO 2 between the surface water and the overlying atmosphere. The transfer velocity parameterisation k = 0.266u 2 10 , where u 2 10 is the wind speed at 10 m, determined by Ho et al. (2006), has been found to match Danish fjord systems (Mørk et al., 2016) and was applied in the current study. Surface values of marine pCO 2 were described by a combination of the open-ocean surface water climatology of pCO 2 by Takahashi et al. (2014) and the climatology developed by Lansø et al. (2015Lansø et al. ( , 2017 for the Baltic Sea and Danish waters. Furthermore, short-term temporal variability was accounted for in the surface water pCO 2 by imposing monthly mean diurnal cycles onto the monthly climatologies following the method described in Lansø et al. (2017).

Meteorological drivers
The necessary meteorological parameters for DEHM were simulated by the WRF (Skamarock et al., 2008), nudged by 6-hourly ERA-Interim meteorology (Dee et al., 2011) for the period 2008 to 2014, and were also used as initial and boundary conditions. In WRF the Noah land surface model, the Eta similarity surface layer and the Mellor-Yamada-Janjic boundary layer scheme were chosen to simulate surface and boundary layer dynamics. The CAM scheme was used for long-and short-wave radiation, the WRF single-moment five-class microphysics scheme was applied for microphysical processes, and the Kain-Fritsch scheme was used for cumulus parameterisation (Skamarock et al., 2008). In WRF the same nests were chosen as in DEHM, and the meteorological outputs were saved every hour. To get the sub-hourly values that match the time step in DEHM, a temporal inter-polation was conducted between the hourly time steps when DEHM was reading the hourly meteorological data. Furthermore, a correction of the horizontal wind speed was conducted in DEHM to ensure mass conservation and compliance with surface pressure (Bregman et al., 2003).

Evaluation of meteorological drivers
We did not have access to measurements from official meteorological observational sites. Thus, for the evaluation of the meteorological drivers, measurements from different types of monitoring sites were used, comprised of three air pollution monitoring sites, three FLUXNET sites, and three sites from the Danish Hydrological Observatory (HOBE; Table 1).
Wind directions, investigated by comparing wind roses made from WRF outputs and measurements, were at most sites reasonably captured by WRF (see Figs. S1-S3 in the Supplement). At several sites the frequency of wind directions from the west were overestimated by WRF, mainly at the expense of southern winds. However, the opposite was the case at Aarhus, where the effect of street canyons was likely causing higher occurrences from due west in the observed wind directions. The wind velocities were in general overestimated by WRF with an average of 1.1 m s −1 , with the greatest differences at the same sites experiencing most problems in reproducing the observed wind direction patterns (Fig. 2). Moreover, at the Risø campus tower site the wind velocities were underestimated.
Only one site had available surface pressure measurements, and high correlation of R 2 = 0.99 was obtained with the simulation surface pressures, indicating that WRF was capable of reproducing the actual pressure system across the study region (Fig. S4). Comparisons of wind velocities and wind rose likewise indicate that WRF captured the general atmospheric flow patterns, though the overestimated wind velocities might induce atmospheric mixing too quickly. The simulated mixing layer heights have previously been evaluated, and although the diurnal boundary layer dynamics were reproduced together with the rectifier effect, problems with accurately modelling the nighttime boundary layer were observed, possibly overestimating nighttime surface concentrations of CO 2 (Lansø, 2016). Moreover, long-range transport and boundary conditions of atmospheric CO 2 concentrations have previously been shown to be captured by the model system across northern Europe, where the current study area is positioned, using observations from Mace Head, Pallas, Westerland, the oil and gas platform F3, Lutjewad, and Östergarnsholm (Lansø et al., 2015).
When evaluating the meteorological variable also important for the biospheric model component, the surface temperature showed high correlation with R 2 above 0.93 for all sites (Fig. 3). The total short-wave incoming radiation (R in ) mirrors the measured R in from the three HOBE sites but is overestimated during summer by WRF (Fig. S5). The values of photosynthetic active radiation (PAR) passed to SPA   from DEHM might thus be overestimated, as PAR is proportional to R in . However, in SPA there is a cap on the limiting carboxylation rate in the calculation of the photosynthesis, and the effect of the overestimated PAR can thus be limited. Precipitation was lacking from all sites, but the annual accumulated modelled precipitation at the nine sites follows the countrywide annual estimates (Cappelan et al., 2018), however, with higher values for the westernmost sites, since many frontal systems enter Denmark from the west (Fig. S6).

SPA
SPA is a mechanistic terrestrial biosphere model (Williams et al., 1996(Williams et al., , 2001. SPA has a high vertical resolution with up to 10 canopy layers (Williams et al., 1996), allowing for variation in the vertical profile of photosynthetic parameters and multi-layer turbulence (Smallman et al., 2013). Within the soil, up to 20 soil layers can be simulated (Williams et al., 2001). The radiative transfer scheme estimates the distribution of direct and diffuse radiation and sunlit and shaded leaf areas (Williams et al., 1998). SPA uses the mechanistic Farquhar model (Farquhar and von Caemmerer, 1982) of leaf-level photosynthesis and the Penman-Monteith model to represent leaf-level transpiration (Jones, 1992). Photosynthesis and transpiration are coupled via a mechanistic model of stomatal conductance, where stomatal opening is adjusted to maximise carbon uptake per unit of nitrogen within hydraulic limitations, determined by a minimum leaf water potential tolerance, to prevent cavitation.  Williams et al., 2005), which is directly coupled to SPA. DALEC simulates carbon stocks in foliage, fine roots, wood (branches, stems and coarse roots), litter (foliage and fine roots) and soil organic matter (including coarse woody debris). Photosynthate is allocated to autotrophic respiration and living biomass via fixed fractions, while turnover of carbon pools is governed by first order kinetics. In addition, when simulating crops, a storage organ (i.e. the crop yield) and foliage pools that are dead but still standing are added, influencing both radiative transfer and turbulent exchange (Sus et al., 2010).
SPA has been extensively validated against site observations from temperate forests (Williams et al., 1996(Williams et al., , 2001, temperate arable agriculture (Sus et al., 2010) and the Arctic tundra (Williams et al., 2000). SPA has more recently been coupled to the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008), and the resulting WRF-SPA model was used in multi-annual simulations over the United Kingdom and assessed against surface fluxes of CO 2 , H 2 O and heat, and atmospheric observations of CO 2 from aircraft and a tall tower (Smallman et al., 2013(Smallman et al., , 2014. SPA needs vegetation and soil input parameters. Initial soil carbon stock estimates were obtained from the Regridded Harmonized World Soil Database (Wieder, 2014). The vegetation inputs and plant traits for SPA were partly taken from previous parameter sets used in SPA but also from the Plant Trait Database (TRY; Kattge et al., 2011) and from literature (Penning de Vries et al., 1989;Wullschleger, 1993). As these parameters and plant traits were determined at various sites that do not necessarily correspond to Danish conditions, a calibration of the vegetation inputs to SPA was conducted for Danish eddy covariance (EC) flux sites (Table 2). Only data from five sites were available, and these were divided in two sets -one for calibration (all available observations before 2013) and the other for validation (all available observations from 2013 and 2014).
In the innermost nest of DEHM for the area of Denmark, a coupling was made between DEHM and SPA. Thus, the coarser optimised biospheric fluxes from CT2015 were, for Denmark, replaced by hourly SPA simulated CO 2 fluxes. With this change, the spatial resolution for the biosphere fluxes was increased from 1 • × 1 • to 5.6 km × 5.6 km, al-lowing for a better representation of the Danish surface and hence also the biospheric fluxes.
On an hourly basis, DEHM provides atmospheric CO 2 concentrations and the meteorological drivers obtained from WRF to SPA, while SPA returns net ecosystem exchange (NEE) to DEHM each hour.

SPA calibration
The calibration was conducted by selecting a set of inputs parameters (plant traits, carbon stocks, etc.), and for each parameter, five values within a realistic range were chosen. Next, 200 SPA simulations with randomly chosen parameter values were conducted. These results were statistically evaluated against observations of NEE from the different flux sites, with the aim of selecting the parameter combination with the lowest root-mean-square error (RMSE) in combination with highest correlation that captured the observed variability and onset of the growing season. However, it was not always possible to have all these conditions satisfied (e.g. Fig. S7). Based on this random parameter testing, it was possible to choose the best set of realistic vegetation input parameters that could improve the model performance at the Danish sites. The best found vegetation parameters values corresponded in some cases to the values already applied in SPA for the given land cover.

SPA evaluation
Comparing to observations of NEE, SPA was, in general, able to capture the phenology and seasonal cycle throughout the entire simulation period (Fig. 4). Correlations and RMSEs between the model and the independent data from the validation period (Fig. 4) likewise indicate a good model performance. At Sorø, variability, as inferred from the standard deviations, the amplitude and the onset of the growing season were reproduced well by SPA. However, difficulties with simulating the evergreen forest at Gludsted are evident, with more variation modelled than given by the observations and a lag of the start of the growing season when compared to the observations. The evergreen plant functional type in SPA lacks a labile or non-structural carbohydrate store needed to drive rapid leaf expansion with the onset of spring; instead leaf expansion is dependent on available photosynthate in a given time step. Therefore, SPA's leaf area index (LAI) is lower early in the growing season, resulting in a biased slow photosynthetic activity and an underestimation in the magnitude of NEE, as seen at Gludsted (Fig. 4). Voulund alternates between winter and spring barley for the calibration period starting with winter barley in 2009. Note that the whole observed time series of NEE at Voulund is shown together with model NEE of both winter and spring barley (Fig. 4c  and d). While the phenology and amplitude are captured well for spring barley at Voulund, SPA is not able to capture the seasonal amplitude of the winter barley that seems to be more sensitive to the meteorological drives, and seasons with harder winters had lower NEE peaks in summer (winter 2010-2011 and 2012-2013). At the grassland site Skjern Enge, NEE is reasonably modelled for winter, spring and the first part of the summer. The difficulties for late summer and autumn arise from the management practices at the site, where both grazing and grass cutting are conducted, limiting NEE (Herbst et al., 2013). Although grazing is included in SPA, it does not simulate the same reduction in NEE.
Examining the performance of SPA at a higher temporal resolution (Table 3), the correlations are better for hourly values than daily for the land cover classification having problems with the phenology (evergreen forest and winter barley) because SPA is capable of reproducing the diurnal variability. For the remaining land cover classifications, R 2 and RMSE are improved when going from hourly to monthly averages of NEE. Zooming in on shorter time windows, the timing of the diurnal cycle is in accordance with measured NEE (see e.g. Fig. S9), but the amplitude is underestimated by SPA.

Results
The model system was run from 2008 to 2014, with the first 3 years regarded as a spin-up period. In the following sections the terrestrial and marine surface fluxes will be presented first, followed by measurements of atmospheric CO 2 from the Risø campus tower that will be used to assess the performance of the DEHM-SPA system and evaluate local impacts from fjord systems on atmospheric CO 2 concentrations.

Biospheric fluxes
As shown in Fig. 5, SPA simulates an east-west gradient in NEE for both January and July in 2011. Larger values of NEE are found in the western part of Denmark, while the islands and eastern Jutland have lower biosphere fluxes. This gradient follows the distribution of the individual land cover classifications (Fig. S8), their phenology and productivity but also reflects the urbanisation, which is denser in the eastern part of the country. During January, total ecosystem respiration dominates NEE. Evergreen forests and grasslands are represented well in western Jutland, and even though these land cover classes have gross primary production (GPP), they are still dominated by total ecosystem respiration, but their total ecosystem respiration can be higher than the other land cover classifications because of the contribution from the autotrophic respiration that depends on GPP in SPA. During July, the productivity is at its highest for all land cover classes dominating total ecosystem respiration, resulting in negative NEE (Fig. 6), and the gradient across the country is more likely to be a result of the urbanisation. Figure 6 shows the average monthly contribution from each land cover classification to the countrywide NEE, which inherently follows their productivity but also reflects the area covered by each land cover type, with highest peaks for winter wheat and grasslands during June. During winter, the spread amongst the land cover classifications is smaller but still has numerically larger monthly fluxes for the land cover classifications with largest area. Integrating over all land cover classifications, the Danish terrestrial land surfaces are a net source of CO 2 to the atmosphere in the months from October to April, with the highest release of 1063±154 GgC per month in December. From May to September, the biosphere is a net sink with a maximum uptake in June of −4982 ± 385 GgC. The total surface exchange of CO 2 between the atmosphere and Danish biosphere is −7337 ±  1468 GgC yr −1 , where winter wheat has the largest contribution with −2342 ± 1045 GgC yr −1 .

Marine fluxes
The air-sea CO 2 exchange in the Danish inner waters experiences large seasonal variations, while the variations in the North Sea are less pronounced, as illustrated by Fig. 7. The mineralisation in winter increases the surface water pCO 2 in the Danish inner waters, resulting in outgassing of CO 2 to the atmosphere, while uptake occurs during spring and summer months following the decrease in surface water pCO 2 due to biological activities. The simulated annual air-sea CO 2 exchange in the 105 000 km 2 covered by the Danish EEZ amounts to −422 GgC yr −1 . However, this number masks large spatial differences and monthly numerical larger fluxes (Fig. 6).
While the North Sea area contained within the EEZ continuously had uptakes in the range −73 to −191 GgC per month with an accumulation of 1765 GgC yr −1 , the fluxes from the near-coastal Danish inner waters varied in the range −46 to 540 GgC per month, releasing 1343 GgC yr −1 to the atmosphere.

Atmospheric CO 2 concentrations
The time series of measured and simulated CO 2 show good agreement (Fig. 8)    the overall variability in the atmospheric CO 2 concentrations and fluxes. Moreover, the higher resolution in both the transport model and surface fluxes results in a better model performance in simulating atmospheric CO 2 concentrations (Fig. S10).
To investigate the origin of the CO 2 simulated at the Risø site, concentration rose plots of simulated atmospheric CO 2 have been made (Fig. 9). The concentration rose shows the wind direction and associated CO 2 concentrations. Division has been made between seasons and daytime and nighttime values, both showing distinct seasonal and diurnal patterns. The highest values of CO 2 are obtained during winter, where very little diurnal variation is seen. During summer the low-est values are obtained in particular during daylight, when photosynthesis occurs.
The individual contribution from fossil fuel emissions, marine and biospheric exchanges to the atmospheric CO 2 (see Figs. S11-S13) indicate that the biosphere contributes most to the variations simulated at Risø (Fig. S12) -both seasonally and daily. Emissions of fossil fuel experience little diurnal variability, but seasonally having the greatest contribution during autumn and winter (Fig. S11). The highest values are seen originating from the sectors encapsulating the city of Roskilde and the capital region. In all seasons, the simulated oceanic contribution is negative, i.e. indicating uptake of atmospheric CO 2 , but the marine contribution is small, with little variation (Fig. S13). The less negative values in autumn and winter may be a result of the simulated outgassing of CO 2 from the Baltic Sea and Danish inner waters during the winter season (Lansø et al., 2015), which, however, is still dominated by the uptake by global open oceans.
The local impact from Roskilde Fjord is difficult to detect in the marine concentration plots. Flux measurements at Roskilde Fjord have shown uptake of CO 2 during spring while release of CO 2 in the remaining seasons (Mørk et al., 2016), which is accurately captured by the modelling system (Lansø et al., 2017). A footprint analysis of the Risø tower has shown that the fluxes from Roskilde Fjord have a contribution to the total CO 2 flux measured at the top of the 118 m high tower, but it is only minor, since fluxes over water are typically an order of magnitude smaller than fluxes over land (Sogachev and Dellwik, 2017). Therefore, we investigated a period with observed large outgassing from Roskilde Fjord -a storm event in October 2013 that was observed to increase the monthly release of CO 2 in the fjord by 66 % (Mørk et al., 2016). The storm event passed Denmark on 28 October 2013, and at 06:00 UTC southerly winds transport air  masses with higher CO 2 towards the Risø site (Fig. 10a), while at the same time a detectable increase in the oceanic contribution to the CO 2 concentration at the Roskilde Fjord system is seen (Fig. 10b). The model system simulates the small peak in the observed atmospheric CO 2 concentrations for 28 October (Fig. 11a) at the Risø site, but distinguishing between contributions from fossil fuel emissions, the biosphere and the ocean to the atmospheric CO 2 concentration at Risø (Fig. 11b) reveals no oceanic impact, and hence there is no apparent influence from Roskilde Fjord during the storm event.

Surface fluxes
The simulated annual uptake by deciduous forest of −284 ± 21 gC m −2 yr −1 for the period 2011-2014 is within the observed range of annual estimated NEE at Sorø from 1996 to 2009, spanning from 32 to −331 gC m −2 yr −1 (Pilegaard et al., 2011). Improvements to the evergreen plant functional type in SPA are needed, and an addition of a labile pool to the evergreen carbon assimilation would omit the seasonal lag (Williams et al., 2005). Such adjustments have already been made to the DALEC carbon assimilation system utilised by SPA (Smallman et al., 2017), substantially improving the representation of terrestrial phenology but not yet being incorporated into SPA. The estimated uptake of −355 ± 41 gC m −2 yr −1 is in the low range of previous estimates of temperate evergreen forests with −402 gC m −2 yr −1 (Luyssaert et al., 2007) and Danish evergreen plantations of −503 gC m −2 yr −1 (Herbst et al., 2011). This could be caused by the slow leaf onset in spring, inhibiting the productivity at the beginning of the growing season. Previous annual estimates at Danish agricultural field sites found carbon uptake of −31 gC m −2 yr −1 estimated from a mixed agricultural landscape (Soegaard et al., 2003) and −245 gC m −2 yr −1 at a winter barley site (Herbst et al., 2011). The SPA-DEHM model system simulated annual uptakes for winter wheat of −252 ± 113 gC m −2 yr −1 , and spring crops of −179 ± 28 gC m −2 yr −1 , while winter barley had a smaller uptake of −82 ± 91 gC m −2 yr −1 with large standard deviation potentially resulting in small annual releases. The calibration and validation (Fig. 4c) show difficulties in simulating the observed NEE during growing seasons for winter barley, par- Figure 9. Concentration roses of modelled atmospheric CO 2 (ppm) at the Risø site for 2011-2014. The wind direction is split into 10 • intervals, and the frequency is indicated by the concentric circles. The colours indicate the CO 2 concentrations with mean removed that have been transported to the site from the given wind directions. ticularly after cold and snow-covered winters. As pointed out in previous studies, the crop modelling component in SPA could likewise be improved, e.g. by inclusion of intraseasonal crops (Smallman et al., 2014).
The current study estimated the Danish grasslands to be a sink of CO 2 with −210 ± 43 gC m −2 yr −1 , which is similar, albeit slightly smaller, than the −267 gC m −2 yr −1 observed at the Skjern Enge grassland site during 2009-2011 (Herbst et al., 2013) and the −312 gC m −2 yr −1 observed at the Lille Valby grassland site, Denmark (Gilmanov et al., 2007). The European grassland study by Gilmanov et al. (2007) found large variation in annual fluxes from grassland driven by environmental conditions and management practices at the sites varying from 171 to −707 gC m −2 yr −1 , but with most sites having an annual uptake of carbon. As seen in Fig. 4e, more work on grassland calibration could have been done, but the conditions and management regimes at Skjern Enge do not necessarily fit the rest of the Danish grasslands. With the chosen parameters, very comparable results were obtained, indicating that such an additional calibration might not be advantageous.
A tilling approach has been used for the land cover classification in the SPA-DEHM modelling framework, including sub-grid heterogeneity in the model system. However, the seven land cover classes do not fully encompass the ecosystem variability in Denmark. Both grassland and agricultural other cover a broad range of subcategories, with both heather and meadow included in the grassland class, while agricultural other contains, for example, vegetables fields, hedgerows, woodland patches and uncultivated land, highlighting the need to adopt approaches allowing for generating novel spatially varying parameter sets (Bloom et al., 2016). Moreover, large urbanised areas are not accounted for in the current classes either. Adding more land cover classifications could give a better and more realistic surface description, if data for both calibration and validation for the lacking land cover classes, preferably from a similar climatic region to Denmark, were available. Compatible marine fluxes to previous estimates are obtained for the study region. On an annual basis, the Danish inner waters were found to be a source of 30 gC m −2 yr −1 , which agrees with most previous studies. Wesslander et al. (2010) estimated Kattegat to act as a small sink of −14 gC m −2 yr −1 based on measurements of water chemistry, while Norman et al. (2013), on the contrary, found a release of 19 gC m −2 yr −1 using a biogeochemical model of the Baltic Sea. Measurements from Danish fjords, on the other hand, consistently point towards these marine areas being annual sources of CO 2 , with values in the range of 41 to 104 gC m −2 yr −1 (Gazeau et al., 2005;Mørk et al., 2016). The current study estimates the North Sea to be a sink of −29 gC m −2 yr −1 , which is very close to previous estimates, both measured and modelled, of −20 and −25 m −2 yr −1 (Thomas et al., 2004;Prowe et al., 2009).

Atmospheric CO 2 and land-sea signals
WRF is in general capable of simulating the observed wind patterns, while the overestimation of the wind velocity could lead to an overestimation of the atmospheric mixing. However, the SPA-DEHM modelling system resembles the synoptic and diurnal variability in the atmospheric CO 2 concentrations measured at Risø campus tower site. The variability at the Risø site is dominated by the biospheric impact and fossil fuel emissions of CO 2 . The signal from Roskilde Fjord is difficult to detect in the simulated CO 2 concentrations. Even when the marine contribution to the atmospheric concentration alone is examined, the Roskilde Fjord signal is hard to distinguish at the Risø campus tower. Moreover, sea breezes from the narrow Roskilde Fjord might be difficult to detect by the model system with its 5.6 km horizontal resolution.
As Roskilde Fjord previously was found by a footprint analysis to have an impact on the atmospheric CO 2 concentration at the top of the tower (Sogachev and Dellwik, 2017), a period with observations of large outgassing from Roskilde Fjord was examined to more clearly envision its impact in the simulated concentration fields. Both the simulated and observed atmospheric CO 2 increased during the storm event on 28 October (Fig. 10a), but no concurrent increase was seen in the oceanic contribution to atmospheric CO 2 at the Risø site (Fig. 10b). This might be explained by the southerly winds that transported the CO 2 released from the fjord northward and away from the Risø campus tower, which is positioned in the southern part of the fjord. Moreover, in this study the increased flux from Roskilde Fjord was only caused by increased wind speed together with the imposed diurnal cycle of marine pCO 2 (the diurnal amplitude for October was approximately 10 µatm), while measurements suggested that an increase in surface water pCO 2 of approximately 300 µatm also sustained the observed CO 2 flux (Mørk et al., 2016). The lack of such an increase in surface water pCO 2 in the current modelling study could explain why no impact on the simulated atmospheric CO 2 is seen from the marine component during the storm event. Thus, the results could indicate that (i) the narrow Roskilde Fjord was not sufficiently resolved in the current model framework, where the horizontal grid resolution is 5.6 km × 5.6 km, (ii) the surface water pCO 2 was not described in enough detail in the model system, (iii) Roskilde Fjord is not in the footprint of the tower during the storm event or (iv) the fjord only has a minor impact on the atmospheric CO 2 concentrations at Risø.
However, the air-sea CO 2 exchange from the Danish inner waters (including all fjord, inner straits and Kattegat) has an impact during winter. Between November and February, the air-sea fluxes from the Danish inner water correspond to 23 %-60 % of the monthly NEE (see Fig. 6). Moreover, the higher values of about 0.5 ppm in the concentration roses of the marine contribution to the atmospheric CO 2 concentrations at the Risø campus tower site in winter likewise emphasise the marine impact, though the outgassing from the neighbouring Baltic Sea also has a contribution. Although the annual total numerical marine fluxes of 1765 gC yr −1 from the North Sea and 1343 gC yr −1 from the Danish inner waters are comparable to the sizes of annual NEE for individual land cover classifications (e.g. deciduous at −987 gC yr −1 , evergreen at −665 gC yr −1 and grasslands at −1467 gC yr −1 ), the air-sea CO 2 fluxes are 1 order of magnitude smaller than the biospheric fluxes, with 30 gC m −2 yr −1 for the Danish inner waters and −29 gC m −2 yr −1 for the North Sea and Skagerak.

Uncertainties in relation to surface exchanges of CO 2
Some of the largest uncertainties lie in the parameters underlying the terrestrial carbon cycle, in particular those governing allocation to plant tissues and their subsequent turnover.
Most often these are based on maps of land cover or plant functional type, but parameter estimation via data assimilation analysis has shown substantial spatial variation in terrestrial ecosystem parameters within plant functional type groupings, with consequences for carbon cycling predictions (Bloom et al., 2016). Increasing the quantity and type of observations available for data assimilation systems can have a significant impact on reducing uncertainty of model process parameters and simulated fluxes (Smallman et al., 2017). In particular, availability of repeated above-ground biomass estimates was able to half the uncertainty of net biome productivity estimates for temperate forests (Smallman et al., 2017). Above-ground biomass estimates are currently available from remote sensing sources (e.g. Thurner et al., 2014;Avitabile et al., 2016), with future planned missions such as the ESA Biomass mission (Le Toan et al., 2011) and NASA GEDI (https://gedi.umd.edu/, last access: January 2019) providing high-quality observations over the tropics and global scales, respectively. While SPA also uses DALEC to simulate carbon allocation and turnover, it is currently impractical to conduct a similar data assimilation analysis to optimise DALEC (or SPA) parameters based on comparison with observations of atmospheric CO 2 concentrations, as this would require repetition of computationally intense simulations of atmospheric transport. While conducting such an analysis remains a future ambition, we consider it to be out of the scope of the cur-rent study, since the terrestrial surface fluxes in this study are constrained by one data stream consisting of EC measurements. This study has focussed on surface fluxes over a relative short time period, and the model framework was capable of producing such fluxes, including their aggregated impact on atmospheric CO 2 concentrations (Fig. 4) with R 2 = 0.77 and RMSE = 4.87 ppm for daily values.
Uncertainties of the marine fluxes can be associated with both the choice of transfer velocity parameterisation, choice of the wind speed product and the used surface water pCO 2 maps. Sensitivity analysis of global transfer velocity parameterisation based on 14 C bomb inventories shows uncertainties of 20 %, while varying the applied wind speed products for these formulation increases the difference in the global annual flux by 40 % (Roobaert et al., 2018). Including empirical formulations of transfer velocity parameterisation in the analysis increased the sensitivity of the wind speed product to nearly 70 %, while the uncertainty of the parameterisation itself rose to more than 200 %. More than a doubling of the annual uptake by the usage of different transfer velocity formulations has likewise been shown for the study region (Lansø et al., 2015), while the choice of the surface water pCO 2 map could change the study region from an annual sink to source of atmospheric CO 2 (Lansø et al., 2017). As shown by Roobaert et al. (2018) the ERA-Interim and the transfer velocity formulation by Ho et al. (2006) used in the present study have a combined uncertainty estimate around 20 %. The improved data-driven near-coastal Danish pCO 2 climatology better reflects the observed spatial dynamics and seasonality in the Danish inner waters (Mørk, 2015), albeit not diminishing the uncertainty related to surface maps of pCO 2 but reducing it.

Conclusions
By usage of the designed mesoscale modelling framework, it was possible to get detailed insight into the spatio-temporal variability in the Danish surface exchanges of CO 2 and the relative contribution from the different surface types. The simulated biospheric fluxes experienced an east-west gradient corresponding to the distribution of the land cover classes, their biological activity and the urbanisation pattern across the country. The relative importance of the seven land cover classes varied throughout the course of the year. Grasslands had a high contribution to the monthly NEE through all seasons, while croplands influence grew from March to July. On an annual basis, winter wheat had the largest impact on the biospheric uptake, with −2342 GgC yr −1 . However, the simulated biospheric uptake could benefit both from model improvement and divisions into more land cover classes. The marine fluxes, being subdivided into the North Sea including Skagerak and the Danish inner waters, had annual fluxes of opposite signs, with the North Sea being a continuous sink of atmospheric CO 2 and the Danish inner waters experiencing Biogeosciences, 16, 1505-1524, 2019 www.biogeosciences.net/16/1505/2019/ small uptake in summer and release of CO 2 during winter resulting in a total marine annual uptake of −422 GgC yr −1 . Good accordance between simulated and observed concentrations was found between modelled and observed atmospheric CO 2 concentrations for 2013 and 2014 at the Risø campus tower. The origin of the modelled CO 2 concentrations at Risø varied, with biospheric fluxes having the largest impact on diurnal variability, while on a seasonal scale fossil fuel emissions also had a dominant role. The local impact from Roskilde Fjord was difficult to detect, while regional impact from the Baltic Sea and Danish inner straits is apparent in winter. The results may indicate that Roskilde Fjord and its localised impact (i.e. at the Risø campus tower site) on atmospheric CO 2 are not adequately resolved in the current model set-up or only have a modest effect. Numerically, the annual fluxes from the North Sea and the Danish inner water were comparable in size to the annual net terrestrial fluxes from the individual land cover classifications.
In order to further examine the air-sea signal at the complex Risø site surrounded by a mosaic of fjord systems, land masses and the Danish inner water, more model experiments could be made, where a larger focus is put on other marine areas than Roskilde Fjord, e.g. the Danish inner straits, Kattegat and the Baltic Sea. Although the total annual marine flux was small, it disguises large monthly variations, and further investigations could help to understand the carbon dynamics in coastal regions. A runoff component in the modelling system would moreover be beneficial for such studies.
Code availability. Scientists with an interest in the atmospheric chemical transport model, DEHM, can contact Jesper H. Christensen (jc@envs.au.dk) with enquiries. Scientists with an interest in the soil-plant-atmosphere model, SPA, can visit its web page (https://www.geos.ed.ac.uk/homes/mwilliam/spa. html, last access: January 2019) or contact Mathew Williams (mat.williams@ed.ac.uk).
Author contributions. ASL and CG designed the experiment. ASL was responsible for coupling the model, which JHC assisted with, and running the experiments. LLS made contributions to the marine set-up, while MW and TSL were responsible for SPA. KP conducted the atmospheric measurements of CO 2 at the Risø campus tall tower. All authors contributed to the discussion of the results and the preparation of the paper.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This study was carried out as part of a PhD study within the Danish ECOCLIM project funded by the Danish Strategic Research Council (grant no. 10-093901). CarbonTracker CT2015 results that were provided by the NOAA ESRL, Boulder, Colorado, USA, from the website http://carbontracker.noaa.gov (last access: January 2016) have contributed to this work. Emission inventories form IER, EDGAR and Aarhus University have likewise made an important contribution. This study has been supported by the TRY initiative on plant traits (http://www.try-db.org, last access: June 2015). The TRY initiative and database is hosted, developed and maintained by Jens Kattge and Gerhard Bönisch (Max Planck Institute for Biogeochemistry, Jena, Germany). TRY is currently supported by DIVERSITAS-Future Earth and the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig. Moreover, this work used eddy covariance data acquired and shared by the Danish Hydrological Observatory, HOBE, founded by the Villum Foundation and by the FLUXNET community. The FLUXNET eddy covariance data processing and harmonisation were carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center and the OzFlux, ChinaFLUX and AsiaFlux offices. We are grateful for many fruitful discussions regarding the atmospheric measurements of CO 2 at the Risø tall tower and its applications in a modelling framework with Ebba Dellwik at DTU Wind Energy.
Review statement. This paper was edited by Trevor Keenan and reviewed by two anonymous referees.