Simulating the atmospheric CO₂ concentration across the heterogeneous landscape of Denmark using a coupled atmosphere–biosphere mesoscale model system

Although coastal regions only amount to 7 % of the global oceans, their contribution to the global oceanic air–sea CO₂ exchange is proportionally larger, with fluxes in some estuaries being similar in magnitude to terrestrial surface fluxes of CO₂.

Across a heterogeneous surface consisting of a coastal marginal sea with estuarine properties and varied land mosaics, the surface fluxes of CO₂ from both marine areas and terrestrial surfaces were investigated in this study together with their impact in atmospheric CO₂ 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⁻¹. While the marine fluxes from the North Sea and the Danish inner waters were smaller annually, with about −1800 and 1300 GgC yr⁻¹, 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₂ at the only tall tower site in the study region. Although, the tower is positioned next to Roskilde Fjord, the local marine impact was not distinguishable in the simulated concentrations. But the regional impact from the Danish inner waters and the Baltic Sea increased the atmospheric concentration by up to 0.5 ppm during the winter months.

1 Introduction

Understanding the natural processes responsible for absorbing just over half of the anthropogenic carbon emitted to the atmosphere will help decipher future climatic pathways. During the last decade, the ocean and the biosphere are estimated to take up to 2.4±0.5 and 3.0±0.8 PgC yr⁻¹ of the 9.4±0.5 PgC yr⁻¹ of anthropogenic carbon emitted to the atmosphere (Le Quéré et al., 2018). The heterogeneity and the dynamics of the surface complicates such estimates.

Biosphere models of various complexity have been developed to spatially simulate surface fluxes of CO₂, but future estimates of the land uptake are bound with large uncertainties (Friedlingstein et al., 2014) that can be attributed to model structural uncertainties, uncertain observations and lack of model benchmarking (Cox et al., 2013; 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₂ 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., 2007, 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₂ 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₂ 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₂ 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 (Regnier et al., 2013). The air–sea CO₂ exchange is in general numerically larger for estuaries than shelf seas (Chen et al., 2013; Laruelle et al., 2010, 2014) and can, for estuaries, be as large as 1958 gC m⁻² yr⁻¹, while continental shelf seas have fluxes in the range of −154 to 180 gC m⁻² yr⁻¹ (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₂ exchange (Regnier et al., 2013). 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₂ (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₂ (pCO₂) 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₂ greatly increase the air–sea CO₂ exchange (Kuss et al., 2006; Rutgersson et al., 2009; Norman et al., 2013).

In this study we aim to simulate surface exchanges of CO₂ 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₂ fluxes. The Danish Eulerian Hemispheric Model (DEHM) forms the basis of the framework, while the mechanistic biospheric soil–plant–atmosphere 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₂ exchange is simulated at a high temporal resolution with the best applicable surface fields of pCO₂ for the Danish marine areas. Tall tower observations are used to evaluate the simulated atmospheric concentrations of CO₂.

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.

2 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² of land, many land–sea 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₂ and allowing for uptake of atmospheric CO₂. In winter, mineralisation increases pCO₂ (Wesslander et al., 2010), and outgassing of CO₂ to the atmosphere takes place. The North Sea is a persistent sink of atmospheric CO₂, where a continental shelf-sea pump efficiently removes pCO₂ 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₂ 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².

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.

A tiling approach with the seven most common biospheric land cover classifications were selected for the current study, including deciduous forest (3348 km²), evergreen forest (1870 km²), winter barley (1211 km²), winter wheat and other winter crops (9269 km²), spring barley and other spring crops (5368 km²), grassland (6924 km²), and agricultural other (3909 km²), 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).

2.1 Observations of atmospheric CO₂

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², has a mean depth of 3 m and is found in the sector 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₂ 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₂ 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.

3 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₂ 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.

3.1 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 (Christensen et al., 2004), 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₂ (Geels et al., 2002, 2004, 2007; Lansø et al., 2015). The CO₂ 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.

3.1.1 Surface fluxes in DEHM

Anthropogenic emissions of CO₂, wildfire emissions and optimised biospheric fluxes from the NOAA ESRL CarbonTracker system (Peters et al., 2007) version CT2015 were used as inputs to DEHM. Their resolution is 1^∘ × 1^∘, with updated values every third hour. Similarly, CT2015 3-hourly mole fractions of CO₂ were read in as boundary conditions at the lateral boundaries of the main domain.

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 Gyldenkærne, 2011). As the European and Danish emission inventories were from 2005 and 2011, respectively, the emissions were scaled to annual national total CO₂ 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₂ emissions.

3.1.2 Air–sea CO₂ exchange

The exchange of CO₂ between the atmosphere and the ocean, $F_{{\mathrm{CO}}_{\mathrm{2}}}$, was calculated by $F_{{\mathrm{CO}}_{\mathrm{2}}}=K{k}_{\mathrm{660}}\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}$, where K is solubility of CO₂ calculated as in Weiss (1974), k₆₆₀ is the transfer velocity of CO₂ normalised to a Schmidt number of 660 at 20 ^∘C, and ΔpCO₂ is the difference in partial pressure of CO₂ between the surface water and the overlying atmosphere. The transfer velocity parameterisation $k=\mathrm{0.266}{u}_{\mathrm{10}}^{\mathrm{2}}$, where $u_{\mathrm{10}}^{\mathrm{2}}$ 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₂ were described by a combination of the open-ocean surface water climatology of pCO₂ by Takahashi et al. (2014) and the climatology developed by Lansø et al. (2015, 2017) for the Baltic Sea and Danish waters. Furthermore, short-term temporal variability was accounted for in the surface water pCO₂ by imposing monthly mean diurnal cycles onto the monthly climatologies following the method described in Lansø et al. (2017).

3.1.3 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 interpolation 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).

3.1.4 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⁻¹, 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.

Table 1 Location of the sites used for evaluation of the meteorological drivers together with the time period from which measurements are used and the meteorological variables (Met var.) included in the analysis. The measurements were obtained from Danish Hydrological Observatory (HOBE), FlUXNET and Department of Environmental Science at Aarhus University (AU).

Figure 2 Scatter plots of measured versus modelled 10 m wind velocity for the nine sites used for evaluation of the meteorological drivers. Hourly average values are used for both simulated and measured wind velocities. Observed average wind velocity (Obs), simulated average wind velocity (Mod), correlation squared (R²), root-mean-square error (RMSE) and bias are shown for each site.

Figure 3 Scatter plots of measured versus modelled 2 m temperatures for the nine sites used for evaluation of the meteorological drivers. Hourly average values are used for both simulated and measured temperatures. Observed average 2 m temperature (Obs), simulated 2 m temperature (Mod), correlation squared (R²), root-mean-square error (RMSE) and bias are shown for each site.

Only one site had available surface pressure measurements, and high correlation of R²=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₂ (Lansø, 2016). Moreover, long-range transport and boundary conditions of atmospheric CO₂ 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² 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).

Table 2 Location, species and land cover classification in the model system for the five Danish eddy covariance (EC) sites used for calibration and validation of SPA.

3.2 SPA

SPA is a mechanistic terrestrial biosphere model (Williams et al., 1996, 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.

Ecosystem carbon cycling and phenology is determined by a simple carbon cycle model (DALEC; 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, 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₂, H₂O and heat, and atmospheric observations of CO₂ from aircraft and a tall tower (Smallman et al., 2013, 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₂ fluxes. With this change, the spatial resolution for the biosphere fluxes was increased from 1^∘ × 1^∘ to 5.6 km × 5.6 km, allowing for a better representation of the Danish surface and hence also the biospheric fluxes.

On an hourly basis, DEHM provides atmospheric CO₂ concentrations and the meteorological drivers obtained from WRF to SPA, while SPA returns net ecosystem exchange (NEE) to DEHM each hour.

3.2.1 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.

3.2.2 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.

Figure 4 Monthly averaged values of measured (black dashed, calibration period; black solid, validation period) and simulated (red) net ecosystem exchange (NEE) for the Danish EC sites with measurements in the simulation period. The shaded areas show the standard deviations for the modelled and measured NEE calculated using hourly fluxes. The model mean (Mean_model), observational mean (Mean_obs), correlation squared (R²) and root-mean-square error (RMSE) for the validation period (2013 and 2014) are shown for each site.

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² 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.

Table 3 Statistic metrics for the validation period (2013–2014) for the fives sites that have measurements of NEE during the validation period for hourly, daily and monthly values. Measured mean (mean_obs), modelled mean (mean_model), correlation squared (R²) and root-mean-square error (RMSE) are shown for each site and temporal resolution.

4 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₂ 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₂ concentrations.

4.1 Surface fluxes

4.1.1 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 5 Net ecosystem exchange (NEE) for January (a) and July (b) 2011.

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₂ 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 $-\mathrm{4982}\pm \mathrm{385}$ GgC. The total surface exchange of CO₂ between the atmosphere and Danish biosphere is $-\mathrm{7337}\pm \mathrm{1468}$ GgC yr⁻¹, where winter wheat has the largest contribution with $-\mathrm{2342}\pm \mathrm{1045}$ GgC yr⁻¹.

Figure 6 Total monthly average of NEE for each land cover classification for the simulation period of 2011–2014 together with the monthly air–sea CO₂ from the Danish marine areas that have been divided into the North Sea and Skagerak (NSSK) and Kattegat and the Danish straits (inner waters).

4.1.2 Marine fluxes

The air–sea CO₂ 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₂ in the Danish inner waters, resulting in outgassing of CO₂ to the atmosphere, while uptake occurs during spring and summer months following the decrease in surface water pCO₂ due to biological activities.

Figure 7 Simulated air–sea CO₂ exchange within the model framework for January (a) and July (b) 2011. The spatial resolution follows those of nest 4 from the DEHM model (i.e. 5.6 km × 5.6 km). The formulation by Ho et al. (2006) was used to calculate the air–sea CO₂ exchange.

The simulated annual air–sea CO₂ exchange in the 105 000 km² covered by the Danish EEZ amounts to −422 GgC yr⁻¹. 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⁻¹, the fluxes from the near-coastal Danish inner waters varied in the range −46 to 540 GgC per month, releasing 1343 GgC yr⁻¹ to the atmosphere.

4.2 Atmospheric CO₂ concentrations

The time series of measured and simulated CO₂ show good agreement (Fig. 8) with R²=0.77 and RMSE = 4.87 ppm for daily averaged time series, demonstrating that the model is capable of capturing the synoptic scale variability. Also, good statistical measures are obtained for the hourly time series with R²=0.71 and RMSE = 5.95 ppm, but the short-term variability was not always fully captured by the model. All in all, the evaluation shows that the model can capture the overall variability in the atmospheric CO₂ 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₂ concentrations (Fig. S10).

Figure 8 Hourly averages and daily averages of modelled and continuously measured atmospheric CO₂ at the Risø site for 2013–2014. The trends have been removed from the time series.

To investigate the origin of the CO₂ simulated at the Risø site, concentration rose plots of simulated atmospheric CO₂ have been made (Fig. 9). The concentration rose shows the wind direction and associated CO₂ concentrations. Division has been made between seasons and daytime and nighttime values, both showing distinct seasonal and diurnal patterns. The highest values of CO₂ are obtained during winter, where very little diurnal variation is seen. During summer the lowest values are obtained in particular during daylight, when photosynthesis occurs.

Figure 9 Concentration roses of modelled atmospheric CO₂ (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₂ concentrations with mean removed that have been transported to the site from the given wind directions.

The individual contribution from fossil fuel emissions, marine and biospheric exchanges to the atmospheric CO₂ (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₂, 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₂ 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₂ during spring while release of CO₂ 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₂ 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₂ 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₂ towards the Risø site (Fig. 10a), while at the same time a detectable increase in the oceanic contribution to the CO₂ concentration at the Roskilde Fjord system is seen (Fig. 10b). The model system simulates the small peak in the observed atmospheric CO₂ 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₂ concentration at Risø (Fig. 11b) reveals no oceanic impact, and hence there is no apparent influence from Roskilde Fjord during the storm event.

Figure 10 (a) Hourly averages of atmospheric CO₂ concentrations including the annual background across Denmark on 28 October 2013 06:00 UTC during the October storm. (b) The contribution from the marine exchange alone to the hourly averaged atmospheric CO₂ concentration on 28 October 2013 06:00 UTC. The less negative values at the Roskilde Fjord system indicate release of CO₂ to the atmosphere.

Figure 11 (a) Hourly averages of modelled and continuously measured atmospheric CO₂ at the Risø site for 19–29 October 2013 with annual means removed. (b) Contributions from fossil fuel emission, oceanic surface exchange and biospheric surface exchange to the atmospheric CO₂ concentration, shown as 1 h averages of modelled concentrations at the Risø site for the same period.

5 Discussion

5.1 Surface fluxes

The simulated annual uptake by deciduous forest of $-\mathrm{284}\pm \mathrm{21}$ gC m⁻² yr⁻¹ 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⁻² yr⁻¹ (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 $-\mathrm{355}\pm \mathrm{41}$ gC m⁻² yr⁻¹ is in the low range of previous estimates of temperate evergreen forests with −402 gC m⁻² yr⁻¹ (Luyssaert et al., 2007) and Danish evergreen plantations of −503 gC m⁻² yr⁻¹ (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⁻² yr⁻¹ estimated from a mixed agricultural landscape (Soegaard et al., 2003) and −245 gC m⁻² yr⁻¹ at a winter barley site (Herbst et al., 2011). The SPA–DEHM model system simulated annual uptakes for winter wheat of $-\mathrm{252}\pm \mathrm{113}$ gC m⁻² yr⁻¹, and spring crops of $-\mathrm{179}\pm \mathrm{28}$ gC m⁻² yr⁻¹, while winter barley had a smaller uptake of $-\mathrm{82}\pm \mathrm{91}$ gC m⁻² yr⁻¹ 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, particularly 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 intra-seasonal crops (Smallman et al., 2014).

The current study estimated the Danish grasslands to be a sink of CO₂ with $-\mathrm{210}\pm \mathrm{43}$ gC m⁻² yr⁻¹, which is similar, albeit slightly smaller, than the −267 gC m⁻² yr⁻¹ observed at the Skjern Enge grassland site during 2009–2011 (Herbst et al., 2013) and the −312 gC m⁻² yr⁻¹ 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⁻² yr⁻¹, 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⁻² yr⁻¹, which agrees with most previous studies. Wesslander et al. (2010) estimated Kattegat to act as a small sink of −14 gC m⁻² yr⁻¹ based on measurements of water chemistry, while Norman et al. (2013), on the contrary, found a release of 19 gC m⁻² yr⁻¹ 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₂, with values in the range of 41 to 104 gC m⁻² yr⁻¹ (Gazeau et al., 2005; Mørk et al., 2016). The current study estimates the North Sea to be a sink of −29 gC m⁻² yr⁻¹, which is very close to previous estimates, both measured and modelled, of −20 and −25 m⁻² yr⁻¹ (Thomas et al., 2004; Prowe et al., 2009).

5.2 Atmospheric CO₂ 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₂ 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₂. The signal from Roskilde Fjord is difficult to detect in the simulated CO₂ 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₂ 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₂ increased during the storm event on 28 October (Fig. 10a), but no concurrent increase was seen in the oceanic contribution to atmospheric CO₂ at the Risø site (Fig. 10b). This might be explained by the southerly winds that transported the CO₂ 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₂ (the diurnal amplitude for October was approximately 10 µatm), while measurements suggested that an increase in surface water pCO₂ of approximately 300 µatm also sustained the observed CO₂ flux (Mørk et al., 2016). The lack of such an increase in surface water pCO₂ in the current modelling study could explain why no impact on the simulated atmospheric CO₂ 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₂ 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₂ concentrations at Risø.

However, the air–sea CO₂ 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₂ 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⁻¹ from the North Sea and 1343 gC yr⁻¹ 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⁻¹, evergreen at −665 gC yr⁻¹ and grasslands at −1467 gC yr⁻¹), the air–sea CO₂ fluxes are 1 order of magnitude smaller than the biospheric fluxes, with 30 gC m⁻² yr⁻¹ for the Danish inner waters and −29 gC m⁻² yr⁻¹ for the North Sea and Skagerak.

5.3 Uncertainties in relation to surface exchanges of CO₂

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₂ 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 current 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₂ concentrations (Fig. 4) with R²=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₂ maps. Sensitivity analysis of global transfer velocity parameterisation based on ¹⁴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₂ map could change the study region from an annual sink to source of atmospheric CO₂ (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₂ 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₂ but reducing it.

6 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₂ 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⁻¹. 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₂ and the Danish inner waters experiencing small uptake in summer and release of CO₂ during winter resulting in a total marine annual uptake of −422 GgC yr⁻¹.

Good accordance between simulated and observed concentrations was found between modelled and observed atmospheric CO₂ concentrations for 2013 and 2014 at the Risø campus tower. The origin of the modelled CO₂ 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₂ 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.