Increasing temperatures over an 18-year period shortens growing season length in a beech (Fagus sylvatica L.)-dominated forest

Abstract

Key message

Using long-term measurements in a mature beech (Fagus sylvatica L.)-dominated forest located in east Belgium, this paper showed that spring and autumn temperature increases during the last two decades led to an earlier end and a shortening of the growing season. These phenological shifts impact negatively but not significantly the forest annual net ecosystem productivity.

Context

The mechanisms controlling temperate forest phenology are not fully understood nor are the impacts of climate change and the consequences for forest productivity.

Aims

The aim of this paper is to contribute to the understanding of how temperate forest phenology and net ecosystem productivity (NEP) interplay and respond to temperature and its evolution.

Methods

Indicators of leaf phenology and productivity dynamics at the start and the end of the growing season, as well as combinations of these indicators (length of the growing season), were derived from a long-term (1997–2014) dataset of eddy covariance and light transmission measurements taken over a mature beech-dominated temperate forest.

Results

The start and the end of the growing season were correlated to spring (and autumn, for the end) temperatures. Despite no trends in annual average temperatures being detected during the observation period, April and November temperatures significantly increased. As a result, an earlier but slower start and an earlier end, inducing a shorter length of the growing season, were observed over the studied period. The first shift positively impacts the mixed forest NEP but is mainly related to the presence of conifers in the subplot. The earlier end of the growing season, more related to beech phenology, negatively impacts the forest NEP. Overall, these two effects partially compensate each other, leading to a non-significant impact on NEP.

Conclusion

Increasing temperatures over the 18-year studied period shortened the growing season length, without affecting significantly the mixed forest NEP. However, as beeches are only affected by the earlier end of the growing season, this suggests a phenologically driven beech productivity reduction in the forest.

1 Introduction

The climate on Earth is changing. Between 1850 and 2013, the global surface temperature increased by 0.8 °C. This increase was more pronounced in Europe (1.3 °C) and even more for Belgium (2.4 °C) (Brouwers et al. 2015). This rise occurred in all seasons but is especially higher in spring (Zohner and Renner 2019), e.g., + 3.0 °C for Belgium (Brouwers et al. 2015). Not only is the absolute temperature increasing, but also the temperature variability across Europe (Huntingford et al. 2013) especially in spring (Matiu et al. 2016). In response, temperate forest phenology is also changing. The last decades have globally seen an advance in the start of the growing season and a corresponding lengthening of the growing season (Liu et al. 2018). At the end of the growing season, either earlier (Fu et al. 2014) or later (Schieber et al. 2017), senescence has been reported (Chen et al. 2019).

In this context, understanding the variability of temperate forest net ecosystem productivity (NEP) requires not only the monitoring of weather conditions but also of forest phenology through the years (Richardson et al. 2012). Phenology can be tracked by using phenological indicators, such as leaf out or fall dates, derived from various datasets (terrestrial survey, phenological camera, satellite pictures, meteorological data, etc.). Thanks to all the data collected and the analysis realized, huge breakthroughs have been made in the comprehension of the mechanisms controlling temperate forest phenology.

At the start of the growing season, phenology is mainly controlled by three factors: temperature, day-length, and winter chilling (Vitasse and Basler 2013; Pope et al. 2014). Different tree species can be more or less sensitive to each of the three factors (Körner et al. 2010). For example, beech (Fagus sylvatica L.) is known to have a lower sensitivity to temperature than species such as Quercus sp. and Picea sp. (Vitasse et al. 2009; Cole and Sheldon 2017), to have higher chilling requirements, to have a higher sensitivity to day-length (Basler and Körner 2014; Zohner and Renner 2015), and thus to have a lower leaf out interannual variability (Zohner and Renner 2014). Other factors may be involved such as a genetic differentiation between tree individuals (Vitasse et al. 2010), tree maturity (McMillan and Goulden 2008), previous year’s weather conditions (Fu et al. 2014), spectral composition of incoming light (Brelsford and Robson 2018), and soil properties (Arend et al. 2015).

At the end of the growing season, contradicting reports are found in the literature concerning the control of leaf senescence and fall (Chen et al. 2019). According to the location and the period considered, leaf senescence has been related to environmental variables (temperature, precipitation, photoperiod) during spring, summer, autumn, or a combination of these periods (Körner et al. 2010; Vitasse et al. 2010; Fu et al. 2014; Estiarte and Peñuelas 2015). Antagonistic effects of growing season and autumn temperature have been reported (Archetti et al. 2013; Liu et al. 2018; Xie et al. 2018; Zohner and Renner 2019). The relative and qualitative roles of the different periods and environmental variables remain unclear, especially for in situ measurements on mature trees. The impacts of climate change on phenological events is therefore not fully understood (Fu et al. 2014; Gallinat et al. 2015; Liu et al. 2016, 2018) nor are the impacts of the phenological events on forest NEP (Richardson et al. 2010; Crabbe et al. 2016; Zohner and Renner 2019).

Temperate forest NEP is known to be related to the incoming solar radiation at daily and seasonal scales. However, beyond these scales, this single meteorological variable does not explain variability of temperate forest NEP (Desai 2014) and other variables need to be considered such as phenology, photosynthetic capacity, previous year’s weather conditions, and biotic stress (Aubinet et al. 2018).

NEP is well-known to interplay with leaf phenology at annual and interannual scales (Morecroft et al. 2003; Keenan et al. 2014), notably in European beech forests (Pilegaard et al. 2011), even if leaf phenology and NEP dynamics may differ (D’Odorico et al. 2015; Fu et al. 2017) due to their different responses to environmental cues (Bauerle et al. 2012; Zohner and Renner 2019). The length of the growing season is one determinant of annual NEP in deciduous broadleaf temperate forests (Barr et al. 2009; Fu et al. 2017). This was shown specifically for pure beech forests in Europe (Granier et al. 2008; Pilegaard et al. 2011). Influences of spring (Richardson et al. 2009) or autumn (Wu et al. 2012; Wu et al. 2013a; Wu et al. 2013b) phenological events on annual NEP have also been reported. It seems clear that patterns differ between sites (Crabbe et al. 2016) and between years, even for a given site, making it difficult to explain the interactions between temperature, phenology, and NEP. Continuous in situ monitoring is still necessary to improve the understanding of current temperate forest phenological shifts (Xie et al. 2018) and their consequences for NEP.

This work uses the phenological and NEP data available from an eddy covariance site (Vielsalm site, east Belgium), a long-term (more than 20 years of continuous monitoring) observation station located inside a mixed temperate forest dominated by beeches, to address the following questions:

• Can variability of leaf phenology (leaf out and fall) and net ecosystem productivity dynamics (carbon uptake increase and decrease) be explained by temperature?

• How variable and interrelated were leaf dynamics and carbon uptake dynamics during the last two decades?

• How did the variability of these dynamics impact the net ecosystem productivity?

2 Materials and methods

2.1 Site description

The Vielsalm site is located in the Ardennes region in eastern Belgium (50° 18′ 18″ N, 5° 59′ 53″ E) at an altitude around 470 m, inside a mature mixed forest. Beech (Fagus sylvatica L.) dominates the forest between 180 and 330° N of azimuthal direction (59% of the trees in a 1-ha wind sector) in the main footprint of the CO₂ flux measurements (Aubinet et al. 2018). Other species present are silver fir (16%), Abies alba Miller; Norway spruce (10%), Picea abies (L.) Karst.; Douglas fir (5%), Pseudotsuga menziesii (Mirb.) Franco; Western hemlock (5%), i.e., Tsuga hetrophylla (Raf.) Sarg.; and Scots pine (5%), Pinus sylvestris L. The maximum leaf area index (LAI_max) was equal to 5 in 2010 and 2011 (Soubie et al. 2016). Beech characteristics are given in Table 1. They were planted in 1908. The mean height for beeches did not change significantly between 1996 and 2014. The understory vegetation is very sparse, if not absent, in the beech subplot. Soils in the region are deep (> 100 cm) dystric cambisol and a slight slope of 3% is present in the north westerly direction. The climate is temperate with an annual mean temperature of 8.4 °C and an annual rainfall of 1000 mm. More site characteristics can be found in previous papers (Aubinet et al. 2001, 2018).

Table 1 Mean canopy height (m) and mean diameter at breast height (cm) in 1996, 2009, and 2014 for 24 beeches higher than 20 m in a 200 × 60 m west-east-orientated rectangle around the tower. Standard deviations are given in parentheses

2.2 Site instrumentation and CO₂ flux computation

In addition to classical eddy covariance (CO₂ concentration and wind speed) and climatic variable measurements (air pressure, rainfall, global and net radiation, moisture and temperature at different levels in the air, etc.), photosynthetically active radiation (PAR) measurements were made using PAR sensors (Delta T QS, Cambridge, UK) located above and under the canopy (Aubinet et al. 2002). From 1997 to 2009, eight PAR sensors were placed below the canopy (PAR_b) at 0.5 m above ground level and placed at 5-m intervals along a north-south transect in the beech subplot. The dynamics of the PAR above and under the cover is illustrated in the Appendix section (Fig. 6). In 2009, the eddy covariance tower was replaced and PAR sensors were removed during the operation. After 2009, two sensors were added and all the sensors were moved and grouped by five in a cross shape, still in the beech-dominated subplot. For the full measurement period, a PAR sensor was installed at the top of the eddy covariance tower (PAR_a). In 2009, all the sensors were recalibrated before being reinstalled. Data were collected every 30 s but were stored as half-hourly means. A more detailed site instrumentation description can be found in the abovementioned papers. In 2009, the tower was raised from 40 to 52 m. The impact of this change on the measurements has been described in detail by Hurdebise et al. (2017). The procedure used for data acquisition, CO₂ flux computation, and correction is described in a recent paper (Aubinet et al. 2018). This analysis focuses on the period from 1997 to 2014. Raw data are available on the FLUXNET website (De Ligne et al. 2015), except the PAR data which are available on Mendeley data (Hurdebise et al. 2019).

Only the data from the wind sector dominated by beech (from 180 to 330° N) were selected in order to limit the influence of tree species heterogeneity on NEP. A site-dedicated data gap filling was developed (Aubinet et al. 2018). Besides heterogeneity, this gap filling takes into account the tower height change and the climate bias due to the weather disturbance regime associated to south-west winds in Western Europe.

2.3 Leaf phenology indicators derived from PAR measurements

Based on the PAR measurements, the relative light transmission (RLT) (Pilegaard et al. 2011) through the forest was computed for each sensor and for each selected half-hour as the ratio $\frac{{\mathrm{PAR}}_{\mathrm{b}}}{{\mathrm{PAR}}_{\mathrm{a}}}$. Data corresponding to low radiation (PAR_b < 5 W/m² or PAR_a < 100 W/m²) were removed. The threshold used corresponds to 10% of the value typically observed during a bright day. By doing so, around 40% of the data was removed. Daily RLT was estimated for each sensor as the median (less sensitive to extreme values (due to sunfleck) than the mean) of the half-hourly $\frac{{\mathrm{PAR}}_{\mathrm{b}}}{{\mathrm{PAR}}_{\mathrm{a}}}$ ratio in a 20-day window centered on the considered day. In order to avoid artificial interannual variability due to sensor weathering or setup changes, these values were normalized (RLT_corr) by subtracting the minimal ratio observed on the corresponding year and dividing the result by the ratio range of this year. Finally, the days of the year (DoY) when 80% and 20% of the maximal RLT_corr value were reached were identified for each sensor, and the mean DoY between sensors were defined as the leaf out start and end (LO_S and LO_E) and the leaf fall end and start (LF_E and LF_S) days, respectively (Fig. 1). The differences between these phenological indicators were then computed (LO_L = LO_E − LO_S and LF_L = LF_E − LF_S). The LO_L and LE_L were referred to as the length of the leaf out and fall periods, respectively. The length of the leafed period (L_L) was defined as LF_S − LO_E. Furthermore, the relationships between leaf phenology indicators were investigated using correlation matrices.

Fig. 1

Schematic view of the dynamics of RNEP. The relative net ecosystem productivity (a) and of RLT. The relative light transmissivity (b). The different leaf phenology indicators derived from RLT are presented in the lower part, while the different productivity dynamics indicators derived from RNEP are presented in the upper part. The abbreviations “LF,” “LO,” “CUI,” and “CUD” refer to “leaf fall,” “leaf out,” “carbon uptake increase,” and “carbon uptake decrease,” respectively. The subscripts “_S,” “_E,” and “_L” refer to the “start,” the “end,” or the “length” of the phenological phases mentioned above

2.4 Productivity dynamics indicators

Productivity dynamics indicators were derived from the dynamics of daily net ecosystem productivity (NEP). Similar to RLT measurements, daily NEP (NEP_day) values were smoothed using a 1-month mobile window and normalized to obtain a relative signal (RNEP) varying between 0 and 1. The length of the smoothing window was larger for RNEP than for RLT_corr (30 instead of 20 days), due to the higher variability observed in daily NEP. The days of the year when RNEP reached 20% and 80% of the maximum of RNEP value for the first and last time were defined as carbon uptake increase start and end (CUI_S and CUI_E), and carbon uptake decrease start and end (CUD_S and CUD_E), respectively (Fig. 1). The length of the carbon uptake period (CU_L) was defined as CUD_S − CUI_E. Similarly, as for RLT measurements, the use of relative instead of fixed thresholds (Richardson et al. 2010) provides indicators that are insensitive to setup changes (in this particular case, the 2009 measurement height change). Furthermore, the relationships between productivity dynamics and between these dynamics and leaf phenology indicators were investigated using correlation matrices.

2.5 Temperature dependencies and temporal trends

Temporal trends of the leaf phenology and productivity dynamics indicators at an interannual scale were analyzed using linear regression.

The relationships between the abovementioned indicators and the mean air temperature during different periods were analyzed using a stepwise regression. In order to do so, a moving and expandable temperature averaging window was used. The start of the averaging window varied between January 1 and the mean value of the phenological indicators considered, with a step of one week (7 days). The length of the averaging window varied between 3 weeks and 5 weeks, with a 1-week step. For each indicator, the correlation coefficients of all the regressions between the considered indicator and temperatures averaged on all possible averaging windows were stored in a 2D matrix (where the columns correspond to the start of the averaging window and the rows correspond to the length of the averaging window) and displayed in a correlogram using a color scale.

Among all the regressions, the best was kept, its residuals were computed, and their correlation with temperatures averaged on all possible temperature averaging windows was again tested and so on until no more significant (p value < 0.05) relationship was found.

Finally, the significance of temperature temporal trends was tested for all the possible temperature averaging windows and the correlation coefficients of these trends were presented in a correlogram.

In addition, chilling and forcing requirements were evaluated by computing the number of chilling and forcing days (Vitasse and Basler 2013), i.e., the days with average temperature below 5 °C from the November 1 to leaf out and above 5 °C from January 1 to leaf out, respectively.

3 Results

3.1 Year to year variability

The mean value and the range of the different phenological indicators are given in Table 2 and synthesized in Figs. 1 and 2.

Table 2 Mean, range, and standard deviation (SD) of the leaf phenology and productivity dynamics indicators observed between 1997 and 2014. Slopes and regression coefficients of the indicators vs. year are also given. In all tables, italicized values correspond to significant (p value < 0.05) correlations, and bold values to highly significant (p value < 0.01) correlations

Fig. 2

Interannual variability of the dynamics of RNEP. The relative net ecosystem productivity (a) and of RLT. The relative light transmissivity (b). The continuous black lines connect the different phenological indicators averaged from 1997 to 2005, while the discontinuous grey lines connect the different phenological indicators averaged from 2006 to 2014. The significant trends observed (indicated by a cross) are described in Table 2

Productivity dynamics indicators are generally significantly earlier (CUI_S, CUD_S) or longer (CUI_L, CUD_L) than the corresponding leaf phenology indicators (LO_S, LF_S, LO_L, LF_L). The difference between these indicators will be discussed in “Section 4.3.” Significant long-term trends were observed. The characteristics of these trends are given in Table 2 and synthesized in Fig. 2, while Fig. 3 presents the time evolution of different phenological indicators.

Fig. 3

Temporal dynamics of the productivity dynamics (a) and leaf phenology (b) indicators. See the caption for Fig. 1 for the definition of the abbreviations. The dotted grey lines indicate the significant temporal trends observed and are described in Table 2

The mean temperatures at the start (from DoY 92 to 120) and at the end (from DoY 295 to DoY 328) of the growing season have both significantly increased by 0.18 and 0.14 °C year⁻¹ (3.2 and 2.5 °C over the studied period), respectively (Fig. 4). At the same time, significant advances were observed for the start of the carbon uptake increase and decrease (1.61 and 2.07 days year⁻¹ (29 and 37 days over the studied period) for CUI_S and CUD_S, respectively), as well as for leaf fall start (0.65 days year⁻¹ for LF_S (12 days over the studied period)). As a result, a significant decrease over time was observed for the length of the carbon uptake and of the leafed period (2.64 and 0.62 days year⁻¹ for CU_L and L_L, respectively). Finally, carbon uptake increase and leaf out lengths both increased significantly (2.18 and 0.50 days year⁻¹ (39 and 9 days over the studied period) for CUI_L and LO_L, respectively) (Fig. 2). No other significant long-term trend was observed, except a clear increase of interannual variability in the start and the end of the leaf out (LO_S and LO_E) after 2006 (Fig. 3).

Fig. 4

Correlograms of the yearly trends observed for the temperature averaged on a mobile window of 3, 4, or 5 weeks. Each pixel (x, y) corresponds to a temperature averaging window spanning y weeks (vertical axis) after the date x (horizontal axis). The color of a pixel (x, y) corresponds to the correlation coefficient between the years and the mean temperature on the corresponding temperature averaging window (starting at the date x and spanning y weeks after). Significant (p value < 0.05) correlation coefficients are indicated by a cross. A red square corresponds to a positive correlation and a blue square to a negative correlation. Three periods were characterized by a significant increase: ~April (e.g., from DoY 92–120, 0.18 °C year⁻¹), in ~June (e.g., from DoY 141–169, 0.10 °C year⁻¹), and in ~November (e.g., from DoY 295–323, 0.14 °C year⁻¹)

3.2 Relations of indicators to temperature

Most of the phenological indicator variability described in the preceding section can be explained by the variability in the mean temperature at some specific periods of about 1 month length which are presented below. For each phenological indicator, correlograms of the first step of the stepwise regression are presented (Fig. 5). The successive significant relationships found using the stepwise regression are presented in Table 3.

Fig. 5

Correlograms for the relationship between the mean temperature on a mobile window of 3, 4, or 5 weeks and the different leaf phenology and productivity dynamics indicators (annotated on the top of the graph). Each pixel (x, y) on each plot corresponds to a temperature averaging window spanning y weeks (vertical axis) after the date x (horizontal axis). The color of a pixel (x, y) corresponds to the correlation coefficient between the indicator considered in a given graph and the mean temperature on the corresponding temperature averaging window (starting at the date x and spanning y weeks after). Significant (p value < 0.05) and highly significant (p value < 0.01) correlation coefficients are indicated by a cross and a target, respectively. A red square corresponds to a significant positive correlation and a blue square to a significant negative correlation. The vertical dotted line indicates the mean value of the considered phenological(s) indicator(s)

Table 3 Results of the stepwise regression between leaf phenology or productivity dynamics indicators (Ind.) and temperature (T_m) averaged over different periods. In each case, the averaging period (day of year (DoY) and corresponding months (shaded)) corresponding to the best correlation is shown for each step of the stepwise regression. Slopes and coefficients of determination are given, as well as the coefficient of determination of the multivariate regression (R_m²). Italicized values correspond to significant (p value < 0.05) correlations, and bold values to highly significant (p-value < 0.01) correlations

The mean temperature in April was the best predictor for LO_S but also, more surprisingly, of CU_L, CUD_S, and LF_S. It was also significantly related to LO_L, LF_E, L_L, and CUD_L.

LF_S was also highly significantly related to the mean temperature in May which was also the best predictor of LO_L variability. Temperatures from April to May affected both LO_E and CUI_E.

In addition to April and May temperatures, LF_S was also positively related to August temperature (also observed in the opposite direction for CUD_L, LF_L, and LF_E) and negatively to October temperatures (also observed for L_L and, in the opposite direction, for LF_L). For CUD_S, a positive relationship to August temperatures was also found but only in the first step of the stepwise regression, the negative relationship to June temperatures being the most significant in the second step. A similar relationship was found for L_L and CU_L. July temperatures were positively correlated to LF_L and LF_E.

CUI_S was only significantly related to January temperatures. No significant relationships were found for CUI_L and CUD_E.

Finally, the number of chilling days always exceeded 100 days, except for one year (2007) and no relationship between the numbers of chilling and forcing days was found.

3.3 Relations between indicators

Table 4 gives the correlations between indicators. All leaf phenology indicators were significantly correlated to their equivalent productivity dynamics indicators, except for the transition period lengths (CUI_L vs. LO_L, CUD_L vs. LF_L). Correlations were also found between the phenological indicators derived from the same dataset. For example, the leaf out start (LO_S) was highly significantly (p value < 0.01) related to the end of the leaf out (LO_E), which itself was significantly related to the leaf fall start (LF_S), while neither LO_S nor LO_E was correlated to LO_L. The significant correlations appearing in Table 4 can be explained by a dependence on temperature for common periods (see also Fig. 5).

Table 4 Correlation coefficients between leaf phenology and productivity dynamics indicators. Italicized values correspond to significant (p value < 0.05) correlations, and bold values to highly significant (p value < 0.01) correlations. See the first column for the definition of the abbreviations of the first line

4 Discussion

4.1 Leaf out

LO_S occurred around DoY 126 (range, 112–137), which is a little bit later than was found by Vitasse and Basler (2013) (DoY 117) but in the same range as observed across Western Europe by Kramer et al. (2017) (DoY 109 to 137). The LO_S interannual variability (standard deviation, 7 days) is also in agreement with the results of Vitasse and Basler (2013) (8 days) who highlight its low value compared with other deciduous species.

Some variability of LO_S has been found in relationship to April temperatures. A strong relationship between these two variables has been found with a temperature sensitivity of 3.2 days °C⁻¹, which is slightly larger than the 2 days °C⁻¹ predicted by Vitasse and Basler (2013).

Chilling requirements appeared to be filled at the site: except for one year (2007), the number of chilling days always exceeded 100 days. In addition, no relationship between chilling day number and temperature forcing, as suggested by Vitasse and Basler (2013), was found. This appears to be in agreement with Fu et al. (2012) who suggest that chilling requirements are generally fulfilled at Belgian sites.

Concerning the time evolution, no significant trends in LO_S were found for the period investigated. However, three years (2007, 2011, and 2014) were characterized by higher April temperatures (12.2, 12.7, and 9.7 °C, respectively, against 7.2 °C on average), and for these years, LO was advanced by about 15 days. As a significant increase of April temperatures (0.18 °C year⁻¹), mainly driven by the exceptional temperatures in the three abovementioned years, was observed at our site and confirmed by independent measurements made by the Royal Meteorological Institute (Brouwers et al. 2015), a further increase in these temperatures is likely which would lead to more precocious leaf out starts in the future. This would be in agreement with long-term observations (from 1982 to 2011) on beech in Western Europe by Fu et al. (2014).

The LO_L mainly depended on the temperatures underwent during leaf out, i.e., in May, being slowed (longer LO_L) under lower temperatures. It was also subject to a significantly positive trend, increasing with time (0.50 days year⁻¹), mainly driven by three years (2010, 2013, and 2014). Despite these three years corresponding to lower than average temperatures, the trend in May temperatures was not significant and we cannot ensure that the observed leaf out slowing can be explained by cooling.

The LO_E interannual variability mainly resulted from that of the LO_S to which it was highly significantly correlated, probably due to the fast leaf unfolding of the beech (13 days on average) compared with other species (Duputié et al. 2015). No significant trend was observed for this variable but a larger variability was observed after 2006 mainly due to the greater April and May temperature variability.

4.2 Leaf fall

LF_S appears around DoY 304 (range, 291–318) and is negatively correlated to April and October temperatures and positively correlated to August temperatures. The first correlation is likely to be indirect and results from the dependencies of LF_S to LO_E that is itself correlated to April temperatures, as discussed above. This would suggest that the impacts of leaf out earliness are propagated until the end of the season, leading to an earlier leaf fall. Such dependency between the beginning and the end of leaf life span in beeches has been emphasized, in particular by Fu et al. (2014) for beech saplings. However, leaf fall timing does not only depend on LO_E as some variability in leaf life span is also observed. Indeed, both L_L and LF_S are negatively correlated to October temperatures, while LF_S (as LF_E) is positively correlated to August temperatures and L_L is negatively correlated to June temperatures.

Different mechanisms may be invoked to explain these dependencies. The negative relationship to June temperatures observed for L_L could be explained by an increase of metabolism rates under higher growing season temperatures and thus to an earlier leaf senescence as suggested by Liu et al. (2018). The positive relationship to August temperatures (i.e., the advance of leaf fall under colder temperatures) may be associated to the low-temperature leaf senescence trigger hypothesis as presented by Delpierre et al. (2009): senescence occurs only when the photoperiod goes under a given threshold and then depends on the accumulation of cold degree days. In contrast, the negative relationship to October temperatures suggests that, once the senescence is sufficiently advanced, there is an inversion of the relationship to temperature. This suggests that different, and sometimes antagonistic, effects of temperature may impact leaf senescence. In large-scale study (53 meteorological stations, several hundreds of phenological stations in Germany) re-examining the possible impacts of meteorological factors on leaf coloring in beech, Estrella and Menzel (2006) confirmed that the two first effects were also observed but only the second was significant.

4.3 Seasonal carbon uptake increase

The differences in the leaf phenology indicators are representative of the beech stand only, whereas the productivity dynamics indicators are representative of the whole subplot captured by the eddy covariance system. Consequently, they are not only representative of beech phenology but also depend on radiation and, despite the precautions taken to limit this effect, on the assimilation of conifers spread across the subplot. As a first consequence, CUI_S (range, 48–118) is on average 53 days earlier and subject to a three times greater variability compared with LO_S. The CUI_S is correlated with January temperatures, which is probably explained by the activation of conifer assimilation as soon as the temperature rises enough at the end of winter. January temperatures are not subject to significant trends, and the trend observed for CUI_S is mainly driven by a cold winter in the beginning of the measurement period (1997) and a warm winter at the end of this period (2014).

CUI_E is well-correlated to LO_E and both dates coincide (DoY 143 and 139, respectively), which appears plausible as the full leaf development in the beeches coincides with the maximal carbon uptake capacity in the whole subplot. Consequently, the characteristics and dynamics of CUI_E and LO_E are close together and, similar to LO_E, CUI_E is most correlated to temperature averages for April and May and is not subject to any significant trend. The sensitivities to April temperatures found for CUI_E (− 3.1 day °C⁻¹) were close to the mean value reported for temperate deciduous forests (− 2.6 day °C⁻¹) (Fu et al. 2017).

4.4 Seasonal carbon uptake decline

The CUD_S occurs quite early, around DoY 209 (range, 171–256). The possible causes of carbon uptake decrease are multiple, including not only leaf senescence in the beech and in the conifers but also a decrease in incoming radiation that becomes perceptible from late July. However, the high correlation of CUD_S to LF_S and to April and June temperatures, similar to LF_S, suggests that despite these two events being separated by 95 days, the senescence of leaves of beeches is highly significantly correlated to CUD_S. This dependence suggests that the propagation of the impact of leaf out earliness affects not only the leaf phenology but also the productivity dynamics, as supported by the dependence between April temperatures and CUD_L or CU_L. This indicates that spring temperature is the dominant explanatory factor for the variability in leaf phenology and productivity dynamics at the site, which is in agreement with a recent large-scale study on phenology in the Northern Hemisphere (Liu et al. 2016).

The other important feature is the positive correlation of CUD_S (and the associated negative correlation of CUD_L) with August temperatures, with higher temperatures inducing shorter CUD_L. Once again, this could be related to the low-temperature leaf senescence trigger hypothesis.

4.5 Seasonal trends and impact on NEP

Although no trend was found for the annual mean temperature, some seasonal averages were subject to significant increasing trends (April, June, November), some others were not (January, May, July, August), and none were subject to a significant decreasing trend. Leaf phenology and productivity dynamics indicators correlated with these temperatures generally followed the same trends, which led to significant advances in some cases. The most important impacts were observed on the indicators of the season end (CUD_S, − 2.07 days year⁻¹, p = 0.03; CUD_E, − 0.88 days year⁻¹, p = 0.05; LF_S, − 0.65 days year⁻¹, p = 0.02, LF_E, − 0.71 days year⁻¹, p = 0.06). A less important advance was found for CUI_S (− 1.61 days year⁻¹, p = 0.04). Despite no significant trend being found for LO_S, its strong dependency on April temperatures does not exclude the possibility of more precocious starts in the future.

All our results suggest that the observed temperature increases, if they lead to an advance of leaf out and carbon uptake start, would mainly shorten the growing period by making CUD and LF occur earlier. The resulting shortening would range between 0.04 and 1.20 (mean 0.62) day year⁻¹ for LL and between 1.30 and 3.99 (mean 2.67) day year⁻¹ for CUL. The latter value appears very high but is probably partly compensated by a lengthening of CUD_L which, even if not significant, leads to a lower decreasing rate of CUD_E.

The potential impacts of these shifts on the net ecosystem productivity (NEP) are partly compensated. On the one hand, an advance of CUI_S would increase NEP while an earlier CUD_S would decrease it. However, as the advance of the start is less pronounced than the advance of the end, the total expected impact of phenological shifts on NEP is expected to be negative, decreasing the forest productivity.

The computation of the area under the curves in Fig. 2 suggests an advance of CUI_S and the associated lengthening of CUI_L would lead to an increase of 1.09 ± 0.66% year⁻¹ of NEP at the start of the growing season. Conversely, the advance in CUD_S and the associated shortening of CU_L would lead to a decrease of 2.64 ± 1.35% year⁻¹ of NEP per year at the middle of the growing season. Finally, at the end of the growing season, the lengthening of CUD_L would lead to a non-significant increase of 0.60 ± 1.02% year⁻¹. The whole impact of the phenological shifts on NEP at our mixed site would thus not be significantly negative (0.95 ± 3.03% decrease of NEP per year). Similarly, no significant trend in NEP had been observed between 1997 and 2016 for the same site by Aubinet et al. (2018). However, as the advance in CUI_S that we observed was mainly due to the early assimilation of conifer spread across the site, it is likely that the NEP decrease would be more pronounced in a pure beech stand. This contradicts reports based on satellite measurements that highlighted a recent greening of global vegetation (Crabbe et al. 2016).

Finally, besides their direct influence on NEP, indirect impacts of leaf phenology on NEP are also possible. In particular, potential future earlier leaf out starts may increase the risk of damage due to late frosts or spring droughts that could negatively impact the annual productivity by reducing the forest photosynthetic capacity (Zohner et al. 2019).

5 Conclusion

The variability observed in phenological indicators from the start of the growing season and, more interestingly, from the end of it was related to the mean temperatures at different periods. Both direct and (intraannual) lagged effects were identified. More specifically, the temperature correlation analysis showed April temperature was the dominant explanatory variable for the leaf phenology and for the productivity dynamics as it affected almost all the phenological indicators, not only at the beginning but also at the end of the growing season.

The leaf phenology and the productivity dynamics were related, especially at the end of the leaf out and at the start of the leaf fall, i.e., when the leaves are present and fully developed but they differ strongly at the start of the leaf out and the end of the leaf fall as a result of the site heterogeneity.

For the period considered (from 1997 to 2014), a significant advance was observed in the leaf fall start and in the carbon uptake decrease start, as well as a decrease in the length of the leafed and of the carbon uptake periods. Based on a temperature correlation analysis, these trends were mainly explained by the April temperature dynamics during the study period. The potential negative impact of this reduction on the mixed forest net ecosystem productivity was partly compensated by a lengthening of the start of the growing season so that, on the whole, the mixed forest productivity was not subject to any trends.

Increasing temperatures over the 18-year studied period shortened the growing season length, without affecting significantly the mixed forest NEP. However, as beeches are only affected by the earlier end of the growing season, this suggests a phenologically driven beech productivity reduction in the forest.

Appendix

Fig. 6

PAR dynamics above (in blue) and under (in red) the cover. Half-hourly data are presented for one of the eight working sensors in 1998 at the Vielsalm Station.