Net carbon ecosystem exchange during 24 years in the Sorø Beech Forest – relations to phenology and climate

Abstract

The carbon sequestration of plants through photosynthesis is responsible for removal of a substantial amount of the man-made CO₂ emissions to the atmosphere. In recent years this so-called land-sink has removed about 30% of the man-made emissions to the atmosphere, with forests being the most important sinks. The land-sink is, however, vulnerable to changes in the environment, such as the atmospheric composition, climate change, and extreme events like storms and droughts. It is therefore important to study the effects of such change on terrestrial ecosystems to provide the basis for predicting the future of the sink.

We here report the results of continuous CO₂ flux measurements over a Danish beech forest during the years 1996–2019. Over the years the forest acted as a sink of CO₂ with a net carbon sequestration ranging from about zero to 400 g C m^–2 yr^−1. We found significant trends in net ecosystem exchange (NEE) (increasing in absolute terms with 15 g C m^–2 yr²), gross ecosystem exchange (GEE) (increasing with 25 g C m^–2 yr^–2), and ecosystem respiration (RE) (increasing with 10 g C m^–2 yr^–2). A prolonged growing season explained 73% of the increase in NEE. The increasing CO₂ concentration in the atmosphere and a subsequent increase in photosynthetic capacity together with warming are the most likely main causes of the increased carbon uptake. The severe drought in the summer of 2018 resulted in a reduction of the annual NEE of 25%.

Keywords:

• deciduous temperate forest
• Fagus sylvatica L.
• carbon sequestration
• long-term study
• phenology
• climate

1. Introduction

Terrestrial ecosystems remove a substantial amount (about 30% in 2018) of the anthropogenically emitted CO₂ from the atmosphere, this is the so-called land sink (Friedlingstein et al., 2019). Since the gross flux of carbon between the terrestrial ecosystems and the atmosphere is more than 10 times that of the emissions from man-made sources (fossil fuels CO₂ emissions and land use change emissions), a small change in the terrestrial sink can be extremely important for the amount of carbon in the atmosphere. Therefore, it is crucial to enhance the knowledge of the drivers of the terrestrial sink and its vulnerability towards man-made and natural changes, such as clear-cutting of forests and extreme events like droughts.

CO₂ is taken up from the atmosphere by plants through photosynthesis, but CO₂ is simultaneously lost by respiration processes, both autotrophic respiration (from plants) and heterotrophic respiration (from organisms like animals, fungi, and microbes). The net removal of carbon from the atmosphere (carbon sequestration) is thus a result of the balance between gross photosynthesis and respiration from all elements of the ecosystem (Keenan and Williams, 2018).

Forests are quantitatively the most important terrestrial ecosystems for carbon sequestration. Pan et al. (2011) found that tropical forest constituted the largest sink, followed by temperate and boreal forests. A recent study by Pugh et al. (2019) found that the predominant sink was regrowing forests in mid-high latitude, rather than tropical forests.

Long-term measurements of leaf area index (LAI) from satellites during 1982 to 2019 in combination with ecosystem models have shown a widespread increase of the growing season LAI over 25% to 50% of the global vegetated area, whereas only 4% showed a decrease (Zhu et al., 2016). It was further found that globally, the increasing CO₂ concentration in the atmosphere (CO₂ fertilization) contributed most to the observed LAI trend, followed by nitrogen deposition, climate change and land-use change.

Experiments with forest trees artificially exposed to elevated CO₂ have shown a positive effect of CO₂ fertilization on CO₂ uptake by tree species. However, the net carbon sequestration is limited by the availability of nutrients (mainly nitrogen) in the soil (Norby and Zak, 2011; Jiang et al., 2020).

In a recent study, (Baldocchi et al., 2018) analysed published data-sets of long-term CO₂ flux measurements (5–18 years) over terrestrial ecosystems. They found that for temperate deciduous forests, the length of the growing season was the determining factor for net ecosystem exchange.

Baldocchi et al. (2018) identifies three main areas of the value of long-term flux carbon flux measurements: (1) Long–term carbon flux measurements can provide an understanding of recovery after man-made (e.g. clear-cutting) or natural (e.g. drought) disturbances. (2) Long-term flux measurements are needed to established the rate at which ecosystems respond to climate change such as warming, increased CO₂ concentration in the atmosphere and other man-made perturbations. (3) Long-term flux measurements are essential to identify trends and distinguish them from random noise. In consequence already established sites with long-term measurements need to continue operation (Baldocchi, 2020). Moreover, long-term flux measurements can provide a data base to distinguish between normal and extreme ecosystem responses to environmental drivers.

In a previous paper describing flux measurements in the Sorø beech forest from 1996 to 2009, Pilegaard et al. (2011) observed a significant increase in gross ecosystem exchange (GEE: 29 g C m^–2 yr^–2) and increasing net ecosystem exchange (NEE: 23 g C m^–2 yr^–2), and a positive but not significant increase in ecosystem respiration (RE: 5 g C m^–2 yr^–2). The length of the carbon uptake period (CUP) increased significantly by 1.9 d yr⁻¹, but there was no significant increase in the period with leaves (LP). Part of the increase in NEE could be explained by the increase in CUP, but the remainder was caused by an increase in net photosynthetic capacity, amounting to 15% over the measurement period.

The purpose of the the present paper is to document and interpret the extended long-term (1996–2019) CO₂ flux data-set from the Sorø beech forest. We investigate the data with respect to the influence of weather and climate on the fluxes and for possible trends over the years. In addition, we demonstrate the value of the data-set for investigating effects of periods with extreme events, such as the severe drought in the summer of 2018. We use the longer time-series to test the following hypotheses:

1. The trend of increasing NEE observed during the first 14 years of measurements has continued.

2. The increasing NEE is dominated by the increase in GEE.

3. The NEE increase is mainly due to a prolonged growing season.

4. Summer drought has a substantial influence on NEE.

2. Materials and methods

2.1. Site

The measurements were carried out at the Danish ICOS station, DK-Sor. The station is located in the middle of the forest “Lille Bøgeskov” near Sorø on the island of Zealand at 55°29${}^{\prime }$13${}^{″}$ N, 11°38${}^{\prime }$45${}^{″}$ E, 40 m above mean sea level.

The climate is temperate maritime and is dominated by westerly winds and frequent passes of frontal systems. The climate is characterized by cool summers and mild windy winters. When easterly winds dominate, they result in cold winters and hot summers. The average annual temperature during the period 1996–2010 was 9 °C and the average annual precipitation during the same period was 640 mm (Reyer et al., 2019).

The soils are brown soils classified as either Alfisols or Mollisols (depending on a base saturation under or over 50%) with a 10–40 cm deep organic layer. The carbon pool in the soil (down to 1 m depth) is 9.3 ± 2.8 kg m^–2 (Wu et al., 2013). The C/N ratio is about 20 in the upper organic soil layers falling to about 10 in the lower mineral layers. The parent material is relatively rich in lime (25–50%). However most of this is leached from the upper horizons of the forest soil, resulting in a low pH (4–5) and a lower base-saturation (Østergård, 2000). The ground water table fluctuates between 0.2 m in winter and at least 2 m below the surface during summer (Ladekarl, 2001).

The forest is dominated by European beech (Fagus sylvatica L.) and immediately around the station the trees were planted in 1921. The forest is managed and the different beech sections are thinned about 20% every 10th year, resulting in an average thinning of 2% per year. The average tree height in the forest part where the mast is located was 24 m in 1996 and increased to 28 m in 2018. The roughness length was 1.8 ± 0.7 m and the displacement height 20.6 ± 4 m in 2001 (Dellwik and Jensen, 2005).

The average tree diameter at breast height (DBH) ranged from 4.7 cm in 1944 to 38.8 cm in 2010, and the stand density from 1767 to 288 stems ha⁻¹ over the same time–span (Reyer et al., 2019). The wood increment calculated on the basis of yield tables (Møller, 1933) was approximately 11 m³ ha⁻¹yr^−1. The peak leaf area index of the canopy was 4–5 m² m^–2 at mid summer.

The main part of the forest surrounding the measurement mast consists of beech trees of varying age but there are also scattered stands of conifers (mainly Norway spruce (Picea abies (L.) Karst.) as well as single trees of other conifers such as European larch (Larix decidua Mill.). In total conifers constitute about 20% of the forest. The terrain is flat and the distance to the forest edge from the measurement site is 0.5–1.6 km depending on direction (Pilegaard et al., 2011).

2.2. Instrumentation

2.2.1. Eddy covariance system

During the course of the measurement, the main concept of the eddy covariance (EC) system was kept. The sampling rate was 10 Hz. The gas analyser (GA) was kept in an hut at the base of the tower. The air is carried from 43 m height through a 50 long tube (8 mm inner diameter) at 28 stdl min⁻¹ to the hut, where it is sub-sampled initially at a flow rate of 1–2 stdl min⁻¹ from 2007 at 3–4 stdl min⁻¹ through the GA. The advantage of this system is that the GA is kept in controlled and stable conditions and maintenance is easy. A disadvantage is that the long tube leads to comparably high spectral attenuation (Ibrom, Dellwik, Flyvbjerg, et al., 2007). For long-term measurements the stability of the system and robustness is an important factor. During the long measurement period, several sensors have been used. The GA model changed in 2006 from LI-6262 (LI-COR Lincoln, USA) to LI-7000 (LI-COR), the sonic anemometer changed from Solent 1012 R2, (Gill Instruments, Lymington, UK) to Gill HS-50 in 2014. While the new sensors had improved performance, the principle properties of the EC system did not change considerably apart from the increase in mass flow through the sensor in 2007, which increased cut-off frequency from 0.14 Hz (first years) to now 0.34 Hz. In 2014 the data acquisition moved from a MS-DOS based IBM AT PC (Pilegaard et al., 2001) to a Campbell CR3000 logger (Campbell, Logan, USA). In none of the cases, the parameters of the spectral correction changed substantially after a sensor change.

2.2.2. Meteorology

Throughout the observation period, wind speed and direction, incoming shortwave radiation, air temperature and relative humidity, precipitation and atmospheric pressure were continuously measured and later augmented with a number of additional sensors of which incoming and transmitted, i.e. below canopy, photon flux density of photosynthetically active radiation is of relevance here to determine the leaved period. The placements and the types of sensors, as well as the period of observation are mentioned in Table 1.

Table 1. Meteorological variables measured at the site.

Variable	Sensor	Unit	Height (m)	Start	End
Wind speed	cup anemometer (P2244, Risø National Laboratory, Risø, Denmark)	m s^–1	57	1996	2013
	sonic anemometer (Gill HS50, Gill, Lymington, UK)	m s^–1	43	2014	2019
Wind direction	wind vane (Risø)	^∘	57	1996	2013
	sonic anemometer (Gill HS50)	^∘	43	2014	2019
Air temperature	Pt-100 (Risø)	^∘ C	57	1996	2013
	HMP155a (Vaisala, Helsinki, Finland)	^∘ C	43	2014	2019
Relative humidity	HMP 45 A, (Vaisala)	%	43	1996	2013
	HMP155a (Vaisala)	%	43	2014	2019
Precipitation	tipping bucket rain guage (Semi Pro, Pronamic, Silkeborg, Denmark)	mm (30 min)^–1	25	1996	2013
	tipping bucket (RAIN-O-MATIC Meteorological Snow & Rain, Pronamic)	mm (30 min)^–1	30	2014	2019
Short wave radiation	pyranometer (CM11, Kipp & Zonen, Delft, The Netherlands)	W m^–2	57	1996	2006
	pyranometer CM14b, (Kipp & Zonen)	W m^–2	57	2006	2013
	pyranometer CM14b, (Kipp & Zonen)	W m^–2	43	2014	2019
Incoming PAR	quantum sensor, LI190 (LI-COR, Lincoln, Nebraska, USA)	µmol m^–2 s^–1	57	1996	2013
	quantum sensor, (LI190)	µmol m^–2 s^–1	45	2014	2019
Transmitted PAR	quantum sensor (LI190)	µmol m^–2 s^–1	4	1998	2007
	quantum sensor (LI190)	µmol m^–2 s^–1	1	2007	2019
	10 quantum sensors (SQ-110-L-15, Apogee Instr. Inc, Logan, USA)	µmol m^–2 s^–1	1	2014	2019
Atmospheric pressure	barometric pressure transmitter PTA 427 (A), (Vaisala)	Pa	2	1996	2013
	barometric pressure transducer (Model 278, Setra Inc., Boxborough, USA)	Pa	2	2014	2019

2.3. Flux calculations

The turbulent fluxes were calculated as described in Pilegaard et al. (2011). After the change in data acquisition, the only change in the data processing was to transform the raw data from the CR3000 ASCII format to the Risø ASCII format, i.e. the format that the software rcpm (Morgenstern, 2000) which was needed to calculate the raw statistics and spectra.

The flow of data processing followed mainly the recently released ICOS protocol (Sabbatini et al., 2018). To summarise, the data were read into rcpm and despiked (Vickers and Mahrt, 1997). For the Gill R2 anemometer the wind components were corrected for flow distortion (Nakai et al., 2006). Then the wind data were rotated (2D rotation, (Moore, 1986)), linearly detrended and, for covariance calculation, shifted to the time lag that yielded maximum covariance. The time lag differed between CO₂ and H₂O, see (Ibrom, Dellwik, Larsen, et al., 2007). The software calculates also the values of the stationarity test (Foken and Wichura, 1996) and binned power and co-spectra (n = 75) based on the FFT algorithm, (Press et al., 1992). In a next step the data quality criteria were determined, where it was distinguished between meteorological data quality constraints (e.g. positive momentum flows) and technical (e.g. flow rate through the gas tubing). To define technical thresholds for sensor malfunctioning, absolute values and ranges were found most effective. The values of the thresholds were identified from frequency distributions and subsequent raw data inspection.

With R-scripts, the final transformations such as cross-wind correction (Gill R2 only), unit transformations, and power- and co- spectral correction were performed. The spectral properties of the GA were analysed for monthly and half annual time periods for CO₂ and H₂O, respectively, after Ibrom, Dellwik, Flyvbjerg, et al. (2007). Other than in Ibrom, Dellwik, Flyvbjerg, et al. (2007), the spectral correction factor was calculated as the ratio of an attenuated model spectrum (Horst, 1997) to the un-attenuated spectrum, where the spectral transfer functions considered both high-pass (Rannik and Vesala, 1999) and low-pass filtering.

The calculated fluxes were then converted into NEE using the concentration at the measurement height, because CO₂ concentration profile data were not available over the entire measurement period. For data gap-filling and flux partitioning, we used the R-package REddyProc (Wutzler et al., 2018) with, to stay consistent with earlier procedures, a constant u* threshold of 0.1 m s⁻¹ and the nighttime look-up table approach. In total 20% of the data were gap-filled with REddyProc, with the exception of one longer data gap from January to April 2014. Here we decided to fill this gap with averages from the previous four years’ data.

In the period from November 2017 to April 2018, the operating LI-7000 sensor” forgot” its factory calibration, which made it ”more sensitive”. The raw data from this period were then converted by inverting the measured concentration values with the wrong polynomial and recalculating the corrected values with the correct polynomial. The user calibration parameters were reconstructed with some remaining uncertainty. However, as the LI-7000 sensor is generally very stable, the span parameter did never change more than 2%, with practically no effects on the covariances. After the correction both variances and covariances reached their expected ranges.

2.4. Derived parameters

2.4.1. Leafed period (LP)

The duration and location of the vegetation period can be defined as the time where beech trees carry leaves (leafed period, LP). This was estimated by the time course of the difference between photosynthetic photon flux density above (Q_a) and below canopy (Q_b). Relative light transmission (LT_r) was calculated as the ratio of Q_b and Q_a from daily averages of half hourly values with $Q_a>=$ 50 μmol m^–2 s⁻¹ throughout the years 1999–2019. Q_a was measured at a height of 57 m and Q_b was measured at a point below the canopy. This location changed over time. A typical LT_r time series has two periods one with high and one with low LT_r, winter and summer, respectively, and transitions between them. Large winter LT_r fluctuations complicate automatic detection of phenological phases.

The breakpoints between foliated and non-foliated periods were calculated by using the R package” segmented” (Muggeo, 2008) retrieving piecewise regressions for the time before bud-break, the time from bud-break to fully un-folded leaves, the peak leaf area during summer, the transition to complete fall-off, and the time thereafter. LP was defined as the time between the day of year when 50% of the peak leaf area was obtained in the spring (LP_start(50)) and the day of year when 50% of the leaf area was lost in the autumn (LP_end(50)). The exact days where these occur were found by applying the relationship of light extinction and leaf area. An example result of this procedure is given in Fig. 1, showing results from 2018, with LP_start(50)) at day 121 and LP_end(50) at day 303, resulting in an LP of 182 days.

Fig. 1. Time course of relative light transmission (LTr, dots), with piece-wise regressions (blue lines). The LTr at peak LAI and at 50% of peak LAI is indicated by dashed green lines. The day of year for 50% peak leaf area in the spring and autumn, respectively, are indicated with vertical lines.

2.4.2. Carbon uptake period (CUP)

The start and end of the carbon uptake period (CUP), i.e. the day of year when the NEE turns from positive to negative in spring, respectively the day of year, when the NEE turns from negative to positive in the autumn were calculated from the time course of daily values of NEE applying a variable span scatterplot smoother (function ” supsmu” in R (R Core Team, 2019)). The start and end of CUP were identified as the day of year, when the smoothed curve crossed zero g m^–2 d^–1. The carbon uptake period (CUP) was calculated as the number of days between the start and end of CUP. Because of a difficulty in determining the start of CUP in 2004 due to a lot of scatter in the fluxes during the spring period, we substituted the value for start of CUP with an estimate based on the relationship between CUP and bud-break. The start of CUP occurs on average only 1 day after bud-break.

2.4.3. Flux footprint

Flux footprints were calculated according to Kljun et al. (2015) using the R implementation. The measurement height (z_m) was calculated as: (1) $$z_m=z_{\mathit{receptor}}-z_d,$$(1) where z_receptor was calculated as the height of the sonic (43 m) minus the height of top of the canopy. The height of the canopy was calculated according to Babst et al. (2014). The displacement height (z_d) was 20.6 m in 2001 according to Dellwik and Jensen (2005). The displacement height was adjusted by an annual height increment of 0.17 m. The boundary layer height (h) was not measured but estimated by the interpolation formula proposed by Nieuwstadt (1981) for neutral to stable conditions: (2) $$h=\frac{L}{3.8}\left[-1+{\left(1+2.28\frac{u_{*}}{fL}\right)}^{1/2}\right],$$(2) where L is the Monin–Obukhov length and f the Coriolis parameter. For unstable situations h was set to 1500 m.

Flux footprints were calculated annually for half-hourly measured data meeting the conditions h > 40 m and $u_{*}$ > 0.15 m^−1. The relative contribution to the total footprint area (2 × 2 km centred at the mast) was calculated for the land cover types: beech forest, coniferous forest, clear-cut (including Christmas tree plantations), agriculture, semi-urban (built-up area), and woodland (scattered woody patches in the agricultural area surrounding the forest).

2.4.4. Photosynthetic capacity

The net photosynthetic capacity was calculated from NEE by an empirical model (Hollinger and Richardson, 2005) based on Michaelis-Menten kinetics; the exact method is documented in Pilegaard et al. (2011). From the model results we derive the photosynthetic capacity at maximum incident light, Q = 1800 μmol m^–2 s^−1, measured above the canopy at the top of the mast.

2.5. Trend analysis

The effect of the summer drought in 2018 was analysed by means of linear trend estimation based on monthly trends during 1996-2017. The observed monthly NEE in 2018 was compared to the predicted values from the monthly time series.

3 Results

3.1. Fluxes

The flux measurements were carried out continuously since June 1^th, 1996. The overall data coverage is very high (85%) after removal of values due to quality control. The main part of the missing values is caused by a major cap from December 2013 to April 2014 due to the collapse of the mast in the severe storm” Bodil”, December 4th - 7th, 2013.

An overview of the daily NEE values are given in Fig. 2. NEE shows a release during winter-time (December–February) of on average 1.9 g C m^–2 d⁻¹ (corresponding to RE) and peak uptake values in June of on average of 6.3 g C m^–2 d^−1. The annual course of daily GEE shows values around zero during wintertime with a change to negative, i. e. uptake starting in March and peaking in June with an average of −14.4 g C m^–2 d^−1, followed by a decrease until the loss of chlorophyll in the leaves in the autumn. The pattern for RE is similar but with opposite sign, with winter values down to a minimum in February of on average 1.8 g C m^–2 d⁻¹ and a peak in June–July of on average 8.2 g C m^–2 d^−1.

Fig. 2. Net ecosystem exchange (g C m^{– 2} d⁻¹⁾ as a function of day of year.

Differences between the years are seen in the absolute values, but also in the details of the seasonal pattern. Some years show a summer minimum of NEE in late July followed by a smaller peak in September. This is seen for several years, but especially pronounced in 2008 and 2018 (Fig. 2). The same pattern is found for GEE.

Figure 3 shows an image plot of daily values of NEE over the years 1996–2019. The period of release and uptake is clarly seen from the reddish and greenish colours, respectively. Superimposed on the image are the start and end of the carbon uptake period (CUP). There is a significantly earlier start of CUP over the years (0.24 ± 0.13 d yr⁻¹), and also a significantly later end of CUP (0.56 ± 0.23 d yr⁻¹). The trend of CUP over the years is highly significant with an annual increase of 0.95 ± 0.22 d yr${}^{-1}$ (Fig. 4). In contrast to CUP, there was no significant increase in the leaf period (LP) (Fig. 4). The leaf out showed a non-significant tendency of a 0.4 d yr⁻¹ earlier start. No tendency was found for leaf fall.

Fig. 3. Imageplot of daily values of NEE (g C m^{– 2} d⁻¹ over the years 1996–2019. Values are missing in the period January 1^st through May 31^st, 1996. Green dots shows beginning and end of carbon uptake period (CUP). Lines show regressions.

Fig. 4. Length (days) of leaf period (LP) and carbon uptake period (CUP) vs. year.

The GEE, NEE, and RE all show significant increasing trends over the years (Fig. 5). GEE increased by 25.4 ± 3.6 g C m^–2 yr⁻¹, NEE by 15.4 ± 2.5 g C m^–2 yr⁻¹, and RE by 9.9 ± 4.0 g C m^–2 yr⁻¹. There is a considerable inter-annual variation shown as scatter around the lines. According to the regressions, GEE has increased from around −1540 g C m^–2 yr⁻¹ to −2020 g C m^–2 yr⁻¹ over the period 1996–2019, NEE from around −55 g C m^–2 yr⁻¹ to −380 g C m^–2 yr^−1, and RE from 1480 g C m^–2 yr⁻¹ to 1730 g C m^–2 yr^−1.

Fig. 5. Annual RE, NEE, GEE g C m^{– 2} yr⁻¹⁾ vs. year. Full line = regression line, dotted lines = 95% significance limits.

There is a highly significant correlation between annual NEE and CUP with an an increase of 11.7 ± 1.7 g C m^–2 yr⁻¹ d⁻¹ (Fig. 6).

Fig. 6. Annual NEE (g C m^{– 2} yr⁻¹⁾ vs. carbon uptake period (CUP (d)).

The long time-series allow for investigating trends within each month. Such analysis is shown for NEE in Fig. 7. For many months there is a significant trend of increasing NEE. The exceptions are April, July, October, November, and December. The strongest trend is seen in June. Also emphasised in the figure are years with a strong summer minimum in late July (possible drought years: 2008 and 2018). It is seen that in 2018, July and August had significantly lower NEE than expected from the trend. For 2008 this is only the case in July. It is also seen, that the NEE in May and June in these two years does not seem to differ significantly from the overall trend. In 2019, the year following the severe summer drought of 2018, there was an exceptionally high RE in June, August and September (Fig. A2).

Fig. 7. NEE by month and year with regression lines and 99.9% confidence limits (shaded areas). Significance of slope indicated for each month. n.s = not significant, * = p < 0.05, ** = p < 0.01, *** = p < 0.001. Years with significant signs of summer drought indicated with different colours.

Similar results are found for GEE with May and June showing equally strong trends (Fig. A1). The GEE in July and August 2018 are much lower than expected from the overall trend, similar to NEE. In 2008, the GEE in July was as expected, but somewhat lower in August. For RE significant trends are found within many of the months (Fig. A2). The RE of July 2018 is much lower than the average, and the RE in August 2018 also comparably low.

The results of the trend analysis of the effects of the summer drought in 2018 (Fig. 8) showed an increased NEE in May and June, however, not significantly different from the predicted values. The NEE in July and August was strongly reduced and significantly different from the predicted values. Overall, the NEE was reduced 25% compared to the value predicted from the trend in 1996–2017. The low NEE values in July and August were somewhat compensated for by higher than average values in May and June and the low RE in July.

Fig. 8. Comparison between observed NEE in 2018 and NEE predicted by the monthly trends during 1996 to 2017.

3.2. Flux footprint

The results of the calculations of source area contribution to the measured flux is shown in Fig. 9. On average over the years, the dominating source was beech forest (84.2%), followed by agriculture (10.1%), coniferous forest (4.7%), clear-cut (0.79%), semi-urban (0.17%), and woodland (0.02%). Over the years there was a slight increase in the source contribution from beech forest (82–86%), a slight decrease in the contribution from agriculture (12–9%), and from coniferous forest (5.5–4%). Figure 10 shows the flux footprint for 2018 superimposed on a land-use map of a 16 km² area centered at the mast. Footprint contour lines are shown in steps of 10% from 50% to 90%.

Fig. 9. Source contribution for the years 1996–2019.

Fig. 10. Flux foot print for 2018. The red dot depicts the tower location. Footprint contour lines are shown in steps of 10% from 50% to 90%. The background shows land use: beech forest (light green), coniferous forest (green), clear-cut (dark green), agriculture (pale yellow), built-up area (blue), and woodland patches (grey).

3.3. Photosynthetic capacity

The maximum rate of net photosynthetic assimilation (Fig. 11) had a range of −20.3 to −34.7 μmol m^–2 s⁻¹ and shows a significant increase over the years from around −24 μmol m^–2 s⁻¹ to around −30 μmol m^–2 s^−1. The average annual increase was 0.24 ± 0.07 μmol m^–2 s^−1, corresponding to a total increase of 23% over the 24 year period.

Fig. 11. Maximum rate of photosynthetic assimilation at Q = 1800 mol m^{– 2} s⁻¹ (F₁₈₀₀) in June each year from 1996 to 2019. Dotted lines are 95% confidence limits.

4. Discussion

4.1. Trends

The trends in GEE and NEE over the period 1996–2019 were less strong than reported for the shorter period 1996–2009 (Pilegaard et al., 2011). Whether this is a result of a slowdown in the later years or just a difference due to increase in the number of observations can not be concluded. GEE shows no signs of a slowdown in increase, but NEE show some signs of less increase in the later years. However, a statistical test (Davies test (Muggeo, 2008)) does not suggest a change in slope. The signs of slowdown in NEE might rather be due to a relative larger increase in RE in the later years.

From an analysis of several temperate deciduous forests with long–term CO₂ flux measurements, Baldocchi et al. (2018) concluded that the length of growing season is the main determinant factor affecting NEE. In our study, we also found a high correlation between NEE and CUP, with CUP explaining 78% of the increase in NEE. Fernandez-Martinez et al. (2017) examined the trends in 23 temperate and boreal forests (including our forest) and found that annual NEE and GPP increased by 8.4 and 11.2 g C m^–2 yr^−1. The increase in RE of 2.9 g C m^–2 yr⁻¹ was not statistically significant.

The annual increase of 12 g C m^–2 d⁻¹ of CUP is somewhat higher than the global average increase (6 g C m^–2 d⁻¹) estimated by Baldocchi (2020).

The measurement error of the annual NEE at our site was estimated to 40 g C m^–2 yr⁻¹ (Wu et al., 2013). The trend observed in the present study (15 g C m^–2 yr^–2) is well above the threshold for detecting a trend exceeding the detectable limit due to random causes for a record of 24 years ($\approx$ 2 g C m^–2 yr^–2), according to Baldocchi et al. (2018).

The trend in NEE is mainly driven by GEE, but reduced by a simultaneous, but smaller increase in RE (Fig. 5).

4.2. Interannual variation

The interannual variation is probably best evaluated by the plots of monthly values and their trends. For NEE (Fig. 7), the largest variations around the regression lines were found in the months May (variation in timing of bud-burst), July (incidence of summer drought), and October (timing of leaf fall). For GEE (Fig. A1, similar to NEE, the largest variations were found in May and July). RE (Fig. A2) behaved somewhat differently with the large interannual variations found in all the months from May through October. Most of the interannual variation is believed to be weather driven, e.g. temperature controls the exact timing of the bud-burst and light and water availability determines the GEE during summertime. RE is driven by both GEE and temperature and soil moisture.

4.3. Global change

In the study of 23 temperate and boreal forests mentioned above, Fernandez-Martinez et al. (2017) found that annual trends of NEE and GEE were mostly correlated with the increasing trend of CO₂. They also point out that the CO₂ fertilization has increased NEE more than RE, which is similar to our findings. The estimated increase in NEE was 4.8 g C m^–2 d⁻¹ppm⁻¹ CO₂. Estimates of global sensitivity to CO₂ ranged from 0.23 to 0.45 C m^–2 d⁻¹ppm⁻¹ CO₂ in a study based on satellite data Fernández-Martínez et al. (2019). The global increase in atmospheric CO₂ over the period of this study (1996-2019) has been around 48 ppm according to the Mouna Loa data (Dr. Pieter Tans, NOAA/ESRL (www.esrl.noaa.gov/gmd/ccgg/trends/) and Dr. Ralph Keeling, Scripps Institution of Oceanography (scrippsco2.ucsd.edu/)). In order to explain the trend in our data, the sensitivity of NEE to CO₂ should thus have been 8 g C m^–2 d⁻¹ppm⁻¹ CO₂. We therefore suspect that the CO₂ fertilization alone cannot explain the observed increase over the 24 years, but that other factors are contributing to the increased uptake, among which the most likely candidates being temperature and precipitation.

Warming has been found to affect phenology with earlier bud-burst (0.46 d yr⁻¹) and later autumn senescence (0.83 d yr⁻¹) (Keenan et al., 2014). While we found both a significantly earlier start of CUP over the years (0.27 d yr⁻¹) and a significantly later end of CUP (0.52 d yr⁻¹), neither the leafed period (LP), nor the timing of leaf-out and leaf-fall showed any significant trends. The start of CUP normally coincides with the bud burst, so CUP start can be used as a proxy for leaf out.

The bud-burst of beech is regulated by a complex combination of light and temperature and is difficult to predict (Vitasse and Basler, 2013). The loss of leaves from the trees occurs after browning and often during a storm. CUP is a result of both uptake and release of CO₂. Thus it involves the processes of photosynthesis, autotrophic respiration and heterotrophic respiration. CUP end is therefore not necessarily a good proxy for autumn senescence. However, there is a chance to gain more insight into the phenology of the beech forest from analysis of images taken from a phenocam above the forest. This study has yet to be completed.

4.4. Drought

The summer drought of 2018 was the most severe drought experienced through the measurement period (Fig. 12). One other year, 2008, also showed a severe drought, but in contrast to the 2018 drought that started in early May and lasted more or less for the rest of the year, the 2008 drought was confined to May–July and was less severe.

Fig. 12. Monthly SPEI index for grid cell covering DK-Sor 1996–2019 (Vicente-Serrano et al., 2010).

Based on our trend analysis we found that NEE was somewhat higher than expected in May and June, but significantly lower in July and August. The effect is similar for GEE with the highest value for May observed throughout the whole time-series and the lowest value for July (Fig. A1). The drought effect is less strong for RE (Fig. A2), but the value for July is much lower than expected. This finding is in accordance with Xu et al. (2020), based on a study of effects of heat waves and drought on North American forests. An early warming and thus high carbon uptake in the spring might compensate for the loss by a drought later in the summer (Wolf et al., 2016). In our study, however, the drought continued throughout the summer with only a little rain in August and September, so a full compensation was not possible.

At the start of the drought in May and in the first part of June, there was a sufficient amount of water available in the soil and the GEE was high due to high solar radiation. From late June to beginning of August, water was clearly a limiting factor resulting in a much lower GEE and a corresponding low NEE.

Liu et al. (2020) found that the 2018 heatwave was different from the European heat-waves in 2003 and 2010 because the coupling between soil moisture and temperature was relatively week in the northern European center of the heat-waves and that forest sites with deeper root profiles were less susceptible to surface drying.

There seems to be a carry-over effect to 2019 with a generally higher RE during the summer months, probably because of accumulated soil organic matter in 2018 due to reduced heterotrophic respiration during the drought that extended throughout autumn. This in turn resulted in a comparably low annual NEE of 2019. Such effects have also previously been reported, e.g. by Thomas et al. (2009), who in a study of the effect of drought in a mature ponderosa pine stand, found that decomposition of additional litter may affect the post-drought carbon balance. Another such lagged effect of drought can be reduced growth in the following year due to lower carbohydrate storage in the year of the drought (Frank et al., 2015).

The summer drought in 2008 did not result in a reduced NEE on annual basis. First of all, the drought was not as severe as in 2018 and because of a sufficient rainfall in August, the total effect was much less.

4.5 N and S deposition

In a study based on long–term flux measurements in 23 forest sites in Europe and North America, Fernandez-Martinez et al. (2017) found that the ecosystem respiration and gross photosynthesis was influenced by a decrease in atmospheric deposition of sulphur. In contrast, they did not find any influence of the simultaneous decrease in nitrogen deposition. Over the period of the measurements at our site (1996–2019), a substantial reduction of deposition of S and N has occurred (Ellermann et al., 2018). The two stations in the Danish air pollution monitoring network closest to our site is Risoe/Lille Valby (37 km NE) and Keldsnor (103 km SW). Both stations are in rural areas, with the station Risoe/Lille Valby somewhat influenced from nearby urban areas, and the station Keldsnor more distant to rural areas and thus less influenced. The amount of deposition of air pollutants at our site is believed to be somewhere in between these two stations. Available data from these stations are shown for NO_x (Fig. 13) and for S (in particles) (Fig. 14). The NO_x concentration in the air was reduced to about 1/3 over the period from 1990 to 2017 and the deposition of S to about 1/3 from 1990 to 2010.

Fig. 13. Concentration of NO_x in the atmosphere (μg m^–3) at two air pollution background stations in Denmark (Ellermann et al., 2018).

Fig. 14. Concentration of S in PM10 (μg m^–3) at two air pollution background stations in Denmark (Ellermann et al., 2018).

The reduction of S and N deposition is expected to reduce soil acidity and thus increase microbial activity. On the other hand, less N input might lead to a limitation for plant growth. Fernandez-Martinez et al. (2017) found that the reduction of sulphur generally increased RE and GEE, but since the increase in RE was higher, it resulted in a reduction of NEE. The decrease in nitrogen deposition had a small but not significant reducing effect on both RE and GEE.

At our site, we believe that the system N is abundant due to the relative high atmospheric input of N. We therefore do not expect an immediate direct effect of a reduction of N deposition. The reduction of S deposition might be important, but due to the relative high limestone content in the soil, the acidifying effect of S has probably not been significant. We can, however, not rule out that there could have been an effect in the upper soil layers. We did not observe a significant increase in RE at the site over the years, so we cannot at present verify an effect of the reduction of S deposition.

4.6. Outlook

Our data show significant increases in both GEE and NEE over the 24 years of measurements. We have observed that although the carbon uptake period is increasing, the period with leaves is not. Of course, the carbon uptake period must be somewhat shorter than the leafed period due to the browning of the leaves and the subsequent slow down in CO₂ uptake. We have, however, observed an increase in the canopy exchange capacity that enhances the GEE. Climate change leads to a general warming and can therefore extend the growing season; the higher concentrations of CO₂ in the atmosphere stimulates growth and reduces the water loss from vegetation. The most likely explanation of our observation is therefore climate change, maybe combined with a recovery from a disturbance some time before we started our measurements.

In the future, we expect a further increase in GEE, mainly caused by CO₂ fertilization. However, it is not likely that the carbon uptake period can increase substantially in the future, because leaf-out of beech trees is also controlled by light. Therefore, the increase in annual GEE might slow down. If we, as predicted, experience more severe and frequent summer droughts, we might even see a reduction in GEE. The evolution of NEE is difficult to predict, because it is the balance between GEE and RE.

5. Conclusion

1. Over the course of the measurements (1996–2019), the GEE of the forest increased from around −1550 g C m^–2 yr⁻¹ to around −2070 g C m^–2 yr⁻¹ and the NEE increased from around zero to about -400 g C m^–2 yr^−1, while the RE remained unchanged at around −1600 g C m^–2 yr^−1.

2. We found significant trends of increases for NEE (15 g C m^–2 yr^–2), GEE (25 g C m^–2 yr^–2), and RE (10 g C m^–2 yr^–2). (Hypotheses 1 and 2).

3. The prolonged growing season explains 73% of the increase in NEE. The prolongation is most likely caused by a combination of different aspects of atmospheric change (CO₂ fertilisation, temperature increase, reduced N and S deposition, and changed precipitation pattern). (Hypothesis 3).

4. The increase in photosynthetic capacity is mainly explained by the increased CO₂ concentration in the atmosphere, but changes in the soil due to reduction in S and N deposition might play a role.

5. A considerable year-on-year variation is due to the actual weather in different seasons and years.

6. The most severe summer drought occurred in 2018, but other years show clear signs of decreased C-uptake in the summer months; especially 2008.

7. The summer drought in 2018 resulted in a much lower GEE than normally. The NEE was moderately affected because of a very low RE. The annual NEE in 2018 was about 100 g C lower than expected from the trend over previous years. (Hypothesis 4).

8. The strength of having a long continuous time-series of carbon flux measurements are clearly demonstrated in the identification of a trend in carbon uptake and the ability to analyse the effects of extreme events like the 2018 summer drought.

Acknowledgements

The work has been funded by multiple sources: National Funding (Danish Environmental Research Programme, Danish Council for Strategic Research, DTU Climate Centre, Danish national project ECOCLIM, Danish Agency for Science and Higher Education), EU Funding (Euroflux, CarboEuroFlux, CarboEurope-IP, CORE, IMECC, Nitro-Europe, and CARBO-Extreme). For allowing us access to the site, we thank the land-owner, Sorø Akademi (Directors Jens Thomsen and Jens Kristian Poulsen and forest manager Anders Grube).

Appendix

Fig. A1. GEE by month and year with regression lines and 99.9% confidence limits (shaded areas). Significance of slope indicated for each month. n.s = not significant, * = p < 0.05, ** = p < 0.01, *** = p < 0.001. Years with significant signs of summer drought indicated with different colours.

Fig. A2. RE by month and year with regression lines and 99.9% confidence limits (shaded areas). Significance of slope indicated for each month. n.s = not significant, * = p < 0.05, ** = p < 0.01, *** = p < 0.001. Years with significant signs of summer drought indicated with different colours.