Impact of winter warming on CO₂ fluxes in evergreen needleleaf forests

Compared to drought and heat waves, the impact of winter warming on forest CO₂ fluxes has been less studied, despite its significant relevance in colder regions with higher soil carbon content. Our objective was to test the effect of the exceptionally warm winter of 2020 on the winter CO₂ budget of cold-adapted evergreen needleleaf forests across Europe and identify the contribution of climate factors to changes in winter CO₂ fluxes. Our hypothesis was that warming in winter leads to higher emissions across colder sites due to increased ecosystem respiration. To test this hypothesis, we used 98 site-year eddy covariance measurements across 14 evergreen needleleaf forests (ENFs) distributed from the north to the south of Europe (from Sweden to Italy). We used a data-driven approach to quantify the effect of radiation, air temperature, and soil temperature on changes in CO₂ fluxes during the warm winter of 2020. Our results showed that warming in winter decreased forest net ecosystem productivity (NEP) significantly across most sites. The contribution of climate variables to CO₂ fluxes varied across the sites: in southern regions with warmer mean temperatures, radiation had a greater influence on NEP. Conversely, at colder sites, air temperature played a more critical role in affecting NEP. During the warm winter of 2020, colder regions experienced larger air temperature anomalies compared to the other sites; however we did not observe a significantly larger increase at colder sites due to winter warming. The varying responses of NEP across different sites highlight the complex interactions between climate variables such as air temperature, soil temperature, and radiation. These findings underscore the importance of integrating winter warming effects to more accurately predict the impacts of climate change on forest carbon dynamics.

1 Introduction

One of the key challenges in assessing the role of forests in mitigating climate change lies in understanding how forest CO₂ fluxes respond to extreme climatic conditions, particularly increases in air temperature. While forests serve as a significant sink for anthropogenic CO₂ emissions (Friedlingstein et al., 2023), extreme warming events can compromise their ability to sequester carbon effectively (Shekhar et al., 2023; Gharun et al., 2024). Although much research has focused on extreme events during the growing season, the impacts of warming winters remain relatively understudied (Kreyling et al., 2019).

In regions where evergreen conifers predominate, such as northern latitudes or higher altitudes, winter warming events can be particularly pronounced (IPCC, 2014). For instance, in 2020, Europe witnessed its warmest winter on record since 1981, with the most significant deviation from the reference period (1981–2020) observed in northeastern Europe (Copernicus Climate Change Service (C3S), 2020). However, the specific effects of such winter warming on CO₂ fluxes, especially in forested areas covered by snow and rich in soil carbon content, remain unclear.

1.1 Effect of warming on forest carbon fluxes

Forest net ecosystem productivity (NEP) depends on the balance between gross ecosystem CO₂ uptake (gross primary productivity, GPP) and emission (ecosystem respiration, R_eco). Both of these flux components are highly sensitive to climate drivers (e.g., air temperature, soil temperature, solar radiation). When canopy structural changes from one year to another are negligible, the interannual variations can be predominantly explained by changes in the climatic conditions (Hui et al., 2003). Net ecosystem productivity can increase or decrease with changes in temperature (Shekhar et al., 2023). In temperature-limited ecosystems for example, an increase in air temperature increases photosynthesis, which leads to higher gross productivity and potentially increased net CO₂ uptake (if respiration does not increase more) (Lin et al., 2021). However with warming and increased temperatures, respiration (autotrophic and heterotrophic) can also increase, and the balance of this with changes in gross productivity could lead to an increase, no change, or a reduction in net CO₂ uptake (Gharun et al., 2020).

Evergreen forests in the Northern Hemisphere contribute significantly to terrestrial carbon (C) storage and exchange (Beer et al., 2010; Thurner et al., 2014). High-latitude evergreen forests have shown an increase in gross primary productivity (GPP) with increasing temperature largely due to longer growing seasons (Myneni et al., 1997; Randerson et al., 1999; Forkel et al., 2016). Multiple other changes under warming however could counteract such an increase for the overall CO₂ uptake capacity of the forest (e.g., due to an increase in ecosystem respiration). In the absence of soil moisture limitation, respiration increases exponentially with an increase in temperature (Law et al., 1999). Additionally, in the presence of winter warming, despite more favorable conditions for photosynthesis, factors such as water stress or photoinhibition caused by high photon flux densities, in combination with low air temperatures, could downregulate photochemical efficiency and negatively affect net photosynthesis, which could reduce gross primary productivity (Troeng and Linder, 1982).

The temperature sensitivity of ecosystem respiration regulates how the terrestrial CO₂ emissions respond to a warming climate. Within naturally occurring temperature ranges, ecosystem respiration (sum of autotrophic and heterotrophic respirations) typically shows an exponential increase with temperature (Lloyd and Taylor, 1994). Previous studies, such as Chen et al. (2020), have demonstrated that Q₁₀, which reflects the rate of increase in soil respiration with a 10 °C rise in temperature, tends to increase as site mean temperatures decrease. However, the temperature sensitivity of ecosystem respiration is influenced not only by the direct effects of temperature on metabolic activity in plants and microorganisms but also by indirect factors such as moisture levels, the leaf area index, photosynthate input, litter quality, and microbial community composition. For example soil moisture affects the microbial activity and decomposition rates, which in turn influence respiration rates. In moist conditions, microbial activity increases, leading to increased decomposition and respiration rates. Conversely, in dry conditions, microbial activity slows down, reducing the respiration rates. The amount of organic matter produced through photosynthesis affects the availability of substrates for microbial decomposition, and higher photosynthate input results in increased carbon availability, stimulating microbial activity and respiration rates (Reichstein et al., 2002; Fierer et al., 2005; Lindroth et al., 2008; Migliavacca et al., 2011; Karhu et al., 2014; Collalti et al., 2020). The temperature response of net ecosystem productivity is the product of the sensitivity of GPP and ecosystem respiration to temperature (Lloyd and Taylor, 1994; Niu et al., 2011), and the temperature sensitivity of respiration (Q₁₀) changes proportionally with the site mean temperature (e.g., higher Q₁₀ at colder sites; Chen et al., 2020).

In forests, the tree canopy serves as a thermal buffer, moderating the microclimatic below by insulating sub-canopy air. This leads to notable differences between sub-canopy and open-air temperatures, with sub-canopy air temperatures averaging 2 °C higher during winter months (Haesen et al., 2021). Teasing apart the effect of soil versus air temperature on fluxes is important firstly because soil and air temperatures can diverge by as much as 10 °C (Lembrechts et al., 2022) and secondly because air temperature and soil temperature affect CO₂ fluxes differently. Warming of the air increases leaf temperature, predominantly affecting the photosynthetic process in the leaves and CO₂ uptake, while soil temperature more directly influences root and microbial activity and thus influences soil respiration and CO₂ release (Berry and Björkman, 1980; Lloyd and Taylor, 1994).

1.2 Importance of the winter period for evergreen needleleaf forests (ENFs)

Environmental cues such as temperature, photoperiod, and light quality control a network of signalling pathways that coordinate cold acclimation and cold hardiness in trees that ensure survival during long periods of low temperature and freezing (Öquist and Huner, 2003; Ensminger et al., 2006). These signalling pathways include the gating of cold responses by the circadian clock, the interaction of light quality and photoperiod, and the involvement of phytohormones in low-temperature acclimation (Chang et al., 2021). Soluble carbohydrates, including sucrose (most abundant), accumulate in response to low temperatures, starting in late autumn and continuing throughout winter (Strimbeck and Schaberg, 2009; Chang et al., 2015). Persistent uninterrupted cold periods thus play an important role in forming the photosynthetic capacity of the trees as warmer winter temperatures increase the chance of photo-oxidative frost damage during earlier stages of the growing season (Gu et al., 2008; Chamberlain et al., 2019) which would compromise the forest's capacity for CO₂ uptake throughout the year (Desai et al., 2016). The risk of photo-oxidative frost damage increases with winter warming, as warmer winter temperatures can lead to an accumulation of photosynthetically active compounds in plants, and when sudden frost events occur during periods of high radiation, the combination of low temperatures and intense sunlight can induce photo-oxidative stress in plant tissues. This occurs because the photosynthetic machinery is still active but the low temperatures impair the plant's ability to dissipate excess energy, leading to the production (and imbalance) of reactive oxygen species (ROSs) that can damage cells and tissues. Photochemical damage can also happen in the case of high radiation, low water content in the leaf tissue, and low temperature, when photosynthesis and protein turnover become inhibited by low temperatures and when non-photochemical, heat dissipation mechanisms are insufficient to deal with excess excitation (hence the negative effect of freezing temperatures after dehardening) (Öquist and Huner, 2003).

Experimental evidence from temperature-sensitive conifers shows that warm spells in winter can induce premature dehardening of buds and result in stunted shoot development in the following spring (Nørgaard Nielsen and Rasmussen, 2008). Additionally, increased respiration due to warming can deplete stored non-structural carbohydrates (NSCs) and tree hydraulic functioning (if combined with drought) and affect tree functioning in spring (Sperling et al., 2015).

The winter of 2019–2020 was reported as the hottest winter in the last 4 decades (1981–2022) across Europe (Copernicus Climate Change Service, ECMWF). When compared to the average conditions, up to 45 fewer winter ice days were detected in eastern Europe and Russia (Copernicus Climate Change Service, Koninklijk Nederlands Meteorologisch Instituut (KNMI)). In Finland, for example, the average air temperature for January and February was over 6 °C higher than the 1981–2010 mean (Copernicus Climate Change Service, ECMWF). In this study we investigated how the exceptionally warm winter of 2019–2020 affected ENFs in Europe and whether increasing winter temperature increased or decreased the carbon uptake of the forest. Our objectives were to (1) evaluate the relative change in air and soil temperature and incoming radiation during the winter of 2019–2020, compared to a 6-year reference period of 2014–2019; (2) quantify the relative changes in the winter CO₂ fluxes across coniferous sites with available ecosystem-level CO₂ flux measurements; (3) tease apart the contribution of air temperature versus soil temperature versus solar radiation to changes in CO₂ fluxes during the warm winter; (4) test the sensitivity of CO₂ fluxes to each of the climatic drivers; and (5) test if the sensitivity of CO₂ fluxes to temperature changed during the warmer winter compared to previous years. Our hypothesis was that warming in winter will lead to a larger negative effect on net ecosystem productivity (i.e., higher CO₂ emissions) across colder forests due to increased ecosystem respiration. We addressed these objectives and tested our hypothesis by exploring ecosystem-level CO₂ fluxes measured with the eddy covariance method over 98 site years in 14 evergreen needleleaf forests distributed from the boreal to the Mediterranean regions of Europe.

Figure 1 Location of the 14 evergreen needleleaf forest (ENF) sites included in this study. Basemap is the MODIS Land Cover product (MOD12Q1, 500 m spatial resolution) showing the distribution of ENFs in Europe in 2020. Elevation of the sites ranges from 4 m a.s.l. (IT-SR2) to 1735 m a.s.l. (IT-Ren).

2 Material and methods

2.1 Site description

We selected 14 evergreen needleleaf forests where continuous CO₂ fluxes and meteorological measurements were available for at least 6 years until the end of 2020. Selected sites were located from the northern to the southern edge of ENF forest distribution in Europe (Fig. 1). The most northern site studied is located in Sweden at 64.2° N (SE-Svb), and the most southern site is in Italy at 43.7° N (IT-SR2). Mean annual air temperature across the sites varies between 1.8 °C (at SE-Ros and SE-Svb) and 15.4 °C (at IT-SR2). Mean annual total precipitation varies from 527 mm (at SE-Nor) to 1316 mm (at CZ-BK1). Elevation ranges from 4 m a.s.l. (at IT-SR2) to 1730 m a.s.l. (at IT-Ren). CZ-BK1 has the largest leaf area index (LAI; 4.52 ± 0.09 SE), and SE-Ros has the smallest (2.59 ± 0.09). Table 1 summarizes the description of sites including their dominant canopy species.

Table 1 Description of the 14 ENF study sites. Mean annual temperature and total precipitation refer to the 2014–2019 period. Mean number of days with snow cover for each site is based on the MODIS satellite observations. Sites are listed in decreasing order of mean annual temperature.

2.2 Dataset

We used the eddy covariance dataset for the warm winter of 2020 processed with the FLUXNET pipeline (compatible with the FLUXNET2015 collection) in this study (Warm Winter 2020 Team and ICOS Ecosystem Thematic Centre, 2022; https://www.icos-cp.eu/data-products/2G60-ZHAK, last access: 10 July 2024) (Pastorello et al., 2020). We included the analysis of soil and air temperature during the spring season at each site to check for any significant changes in the climate immediately after the winter season. Winter months included December, January, and February, and spring months included March, April, and May. The 6-year reference period was from 2014–2019. This period was selected to have sufficient temporal overlap between the sites. Net ecosystem exchange (NEE) quality-checked with a constant friction velocity (u^*) threshold was used for all sites (NEE_CUT_REF) (Shekhar et al., 2023). For an easier interpretation, we present net ecosystem exchange as net ecosystem productivity (NEP = −NEE), where a negative NEP indicates that the forest is a net source and positive NEP indicates that the forest is a net sink of CO₂ (Chapin et al., 2006).

In terms of climatic variables we selected those that overlapped in data availability across all sites during the study period. These included incoming shortwave radiation (R_g), air temperature (T_air), soil temperature at 5 cm (T_soil), and precipitation and topsoil water content. Given that continuous long-term snow depth measurements were not available at all sites, we used remotely sensed snow depth products to quantify mean snow depth and snow depth anomalies in winter 2020. The snow depth data were derived from the simulation of the Famine Early Warning System Network (FEWS NET) Land Data Assimilation System (FLDAS) (McNally et al., 2017). FLDAS data are produced from the Noah version 3.6.1 land surface model (LSM) at a monthly resolution with a global coverage at a spatial resolution of 0.1°×0.1° (approx. 10 km × 10 km) (Kumar et al., 2013) and have been used in the past to study global spatiotemporal patterns of snow depth and cover (Notarnicola, 2022). For snow cover we used the MODIS/Terra (MOD10A2) and MODIS/Aqua (MYD10A2) (Hall and Riggs, 2021) Snow Cover 8-Day L3 Global 500m SIN Grid, Version 6, dataset, which provides maximum snow cover extent at 8 d temporal resolution and 500 m spatial resolution. As a quality check, we compared the measured snow depth against the remotely sensed snow depth for one site (DE-Tha) for which measurements were available during the study period and found a reasonable agreement between the two datasets (r=0.86, p<0.001). For each forest site, we derived the average (2014–2019) leaf area index (LAI) from the LAI Collection 300 m Version 1.1 product (LAI300) provided by the Copernicus Global Land Service (Fuster et al., 2020). The average LAI was estimated for each site during the mean net CO₂ uptake period. Following Shekhar et al. (2023), the start of the net carbon uptake period was defined as when daily NEP crosses from negative to positive and the end is the inverse.

2.3 Statistical analysis

We compared average daily and daytime (when R_g>10 W m⁻² and 08:00–18:00 local time) means of each variable (v; climate drivers, CO₂ fluxes) during the winter and spring of 2020 to the mean from a 6-year reference period (2014–2019) using a t test (p<0.05). Daily means of each variable were calculated only using the measured and good-quality gap-filled half-hourly data (variable quality control = 0 or 1). To understand the major drivers of winter NEP, R_eco, and GPP for each forest site, we derived the conditional variable importance (CVI_v) of each predictor variable (R_g, T_air, and T_soil) based on a random forest regression model built for the site (Breiman, 2001). For training the random forest model of R_eco, we additionally used GPP as an explanatory variable. In addition to the influence of abiotic drivers, the empirical relationship between photosynthesis (and thus GPP) and ecosystem respiration in forests has been established by a large body of research (Brüggemann et al., 2011; Koerner, 2013; Migliavacca et al., 2011; Shekhar et al., 2024). Soil water content (SWC) was removed from the drivers analysis (1) because of its negligible effect on the overall model (see details below), (2) since not all sites had complete measurements throughout the study period, (3) and because soil water content measurements at freezing soil temperature levels are not reliable and we observed that for several sites soil temperature in winter remained near or below 0 °C (Fig. S1). The effect of soil water content on the random forest (RF) model proved to be negligible after we compared the random forest results once with and once without SWC. The comparison revealed that the difference in the variance explained (R²) was less than 3 %, indicating a negligible improvement in model performance based on the percentage of variance explained (see Fig. S2).

We tuned the random forest model by iterating the “n_tree” parameter (number of trees to grow) from 100 to 500 with steps of 50 and “m_try” parameter (number of variables to try at each split) from 1 to 3 with steps of 1 and chose the parameter (n_tree = 300 and m_try = 2) with the minimum mean square error. CVI_v accounts for the correlation between the predictor variables and was calculated using the party R package (Hothorn et al., 2006). Based on a 7 d moving window (centered on the central value of the window) we calculated the mean daily (and daytime) NEP, T_air, R_g, and T_soil. To compare the CVI_v across sites, for each site we calculated the relative CVI (RCVI) for each variable as per Eq. (1).

$$\begin{array}{}\text{(1)}& {\mathrm{RCVI}}_{\mathrm{v}}\left(\mathit{\%}\right)={\displaystyle \frac{{\mathrm{CVI}}_{\mathrm{v}}}{\sum {\mathrm{CVI}}_{\mathrm{v}}}}\times \mathrm{100},\end{array}$$

where ∑CVI_v is the sum of the CVI_v of all variables used in the model. We expressed changes in the variable during 2020 (v₂₀₂₀) and the reference period (v_reference) based on its relative anomaly (Δv_r) and absolute anomaly (Δv_a) as per Eqs. (2) and (3).

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }{v}_{\mathrm{r}}\left(\mathit{\%}\right)={\displaystyle \frac{v_{\mathrm{2020}}-v_{\mathrm{reference}}}{\left|v_{\mathrm{reference}}\right|}}\times \mathrm{100}\text{(3)}& {\displaystyle }\mathrm{\Delta }{v}_a=v_{\mathrm{2020}}-{\underset{\mathrm{¯}}{v}}_{\mathrm{reference}}\end{array}$$

To better understand how absolute anomalies in different variables (R_g, T_air, and T_soil) contributed to variations in ΔNEP, we applied the RCVI (as described in Eq. 2) derived from a random forest regression model. This model used ΔNEP as the response variable and R_g, T_air, and T_soil as predictors, with hyperparameters set to n_tree = 100 and m_try = 3 (optimized for the lowest mean squared error). The analysis was conducted at each site, with a minimum of 80 d of data points. The percent variance explained of the model (r²) was based on the out-of-bag estimates. The sensitivity of NEP anomalies in winter (ΔNEP) to anomalies of climate drivers (R_g, T_air, and T_soil) was quantified as the slope of ΔR_g, ΔT_air, and ΔT_soil when regressed with ΔNEP using a multivariate linear regression (ΔNEP ∼ $\mathrm{\Delta }{R}_g+\mathrm{\Delta }{T}_{\mathrm{air}}+\mathrm{\Delta }{T}_{\mathrm{soil}}$).

3 Results

3.1 Conditions of the warm winter of 2019–2020 across different sites

According to the in situ data, compared to the reference period (2014–2019), winter 2020 was the warmest winter across 10 sites. At seven sites, the winter was also drier than normal (Fig. S3). Positive air temperature anomalies in winter 2020 were significantly larger at sites with a lower mean (2014–2019) air temperature (p<0.05, $r=-\mathrm{0.53}$), with the largest significant anomaly of 4.79 °C at RU-Fyo and lowest significant positive anomaly of 0.87 °C observed at IT-SR2 (Fig. 2). Incoming shortwave radiation did not change significantly across any of the sites during the warm winter (data not shown here).

Figure 2 Seasonal changes in air temperature (T_air) and soil temperature (T_soil) in 2020 compared to the 6-year reference period (2014–2019). Asterisks mark where means in 2020 were significantly different from the reference period (p<0.05). Anomalies were calculated from daily values. Sites are listed from the top to the bottom in decreasing order of site mean annual air temperature.

The average number of snow cover days per year was highly variable across the study sites (Table 1). The southernmost site studied here (IT-SR2) has no snow cover in winter, while the sub-Alpine forest in Switzerland (CH-Dav) has snow cover on 139 d yr⁻¹ on average (Table 1). At those sites with consistent snow cover in winter (11 out of 14 sites), snow depth declined at 9 out of 11 sites during the warm winter of 2020, and this reduction was considerable at FI-Let, RU-Fyo, SE-Nor, DE-Obe, DE-RuW, and DE-Tha (Fig. 3). At SE-Svb, FI-Let, and DE-Obe soil temperature at 5 cm was continuously above the freezing level in winter 2020 (Fig. S1), unlike the mean conditions at the sites where soil temperature fluctuates around 0 °C in winter. Changes in air and soil temperature were more significant in winter than in spring (Fig. 2), which is the reason why we focus on the effect of winter warming on CO₂ fluxes.

Figure 3 December–May snow depth changes in winter 2020 compared to the average winters during the reference period (2014–2019). Note that only 11 out of the 14 sites have persistent snow cover in winter. Sites are ordered from the top left to the right by increasing the site mean temperature (SE-Svb coldest and DE-Tha warmest).

3.2 Effect of climate drivers on winter CO₂ fluxes

The annual NEP of the ENFs varied from a maximum sink (± SD) of 797 (± 320) g C m⁻² yr⁻¹ (CZ-BK1) to a maximum source of −311 (± 93) g C m⁻² yr⁻¹ (SE-Nor) during the 6-year reference period (2014–2019) (Table 2). Interannual variation in NEP was largest at CZ-BK1 (320 g C m⁻² yr⁻¹) and lowest at SE-Svb (35 g C m⁻² yr⁻¹) (Table 2). The length of the net CO₂ uptake period was on average 178 d but varied between the sites from 105 d (at RU-Fyo) to 315 d (at DE-RuW) (Table 2). Except FR-Bil and DE-RuW, all sites were a CO₂ source in winter under reference conditions (Table S1).

Table 2 Mean total annual net ecosystem productivity (NEP) and the standard deviation (interannual variation) during the reference period (2014 and 2019). The start of the net carbon uptake period (start of season, SOS; day of year, DOY) is when daily NEP changes from negative to positive, and the end (end of season, EOS) is the inverse (following Shekhar et al., 2023). Sites are listed in decreasing order of mean annual air temperature.

During the warm winter of 2020, mean daily NEP (i.e., annual winter CO₂ sink or source strength) changed significantly (p<0.05) at 9 out of the 14 sites (BE-Bra, CZ-BK1, DE-Obe, FI-Let, IT-Ren, IT-SR2, SE-Svb, SE-Nor, and RU-Fyo, grouped together as the “affected” sites) compared to the 2014–2019 reference period, with changes in both positive and negative directions (Fig. 4). For example, at BE-Bra, DE-Obe, IT-Ren, SE-Svb, and FI-Let, the forest became a significantly larger source of CO₂ in winter 2020, while at SE-Nor, CZ-BK1, and RU-Fyo the forest shifted towards being a smaller source for CO₂ and at IT-SR2 it turned into a net sink in winter 2020 (Fig. 4, Table S1). IT-SR2 showed the largest increased daily NEP in winter (331 %), and BE-Bra showed the largest decline in daily NEP (−98 %) (Fig. 4). During the warm winter, ecosystem respiration (approximated by nighttime NEP) increased significantly across 10 out of the 14 sites, indicated by a negative anomaly in nighttime NEP (Fig. 4). Daytime NEP however (dominated by productivity) increased significantly with warming at only five sites (Fig. 4).

Figure 4 Relative change (anomaly, %) in mean daily, mean nighttime, and mean daytime NEP in winter 2020 compared to the 6-year reference winters (2014–2019). Asterisks mark where the mean in 2020 was significantly different from the reference period (p<0.05). Positive NEP change indicates increased net uptake (due to increased uptake or reduced emission), and negative change indicates decreased net uptake (due to reduced uptake or increased emission). Sites are listed from the top to the bottom in decreasing order of mean annual air temperature.

Figures 5–7 show the relationship between air temperature, soil temperature, and incoming shortwave radiation with NEP. Winter CO₂ fluxes in general (for the reference period of 2014–2019) showed a clear decline in NEP in response to an increase in soil temperature across several sites (SE-Svb, SE-Ros, FI-Let, SE-Nor, DE-Obe, BE-Bra, IT-SR2) (Fig. 7).

While the response of NEP to R_g was more consistent across sites, the effect of soil and air temperature on NEP varied largely across sites. The average variation explained by the random forest regression for NEP in winter was 78 % (Fig. S4). The relative importance results of the random forest regression analysis showed that across tested variables, R_g generally had the largest effect on NEP. However, with a decrease in the site baseline (i.e., mean) temperature, the effect of R_g declined (Fig. 8). For example, at the three coldest sites (SE-Svb, CH-Dav, IT-Ren) R_g had a relative importance of 52 %, 23 %, and 41 % for the variations in NEP respectively, while at the three warmest sites (IT-SR2, FR-Bil, BE-Bra), R_g had a relative importance of 73 %, 81 %, and 58 % for NEP respectively (Fig. 8). When looking into partitioned fluxes, radiation dominated the effect on winter GPP and temperature dominated the effect on winter respiration fluxes (Fig. 8). Particularly at the colder sites was the effect of radiation the least important (Fig. 8).

Figure 5 Comparison of the binned response of NEP versus R_g (incoming shortwave radiation) during the winters of the reference period (2014–2019) and winter 2020 across all sites (arranged from the top left to the bottom based on increasing mean air temperature). The daily mean NEP is aggregated (mean ± 95 % CI as error bars, confidence interval) at 10 W m⁻² R_g bins.

Figure 6 Comparison of the binned response of NEP versus T_air (air temperature) during the winters of the reference period (2014–2019) and winter 2020 across all sites (arranged from the top left to the bottom based on increasing mean air temperature). The daily mean NEP is aggregated (mean ± 95 % CI as error bars) at 1 °C T_air bins.

Figure 7 Comparison of the binned response of NEP versus T_soil (soil temperature) during the winters of the reference period (2014–2019) and winter 2020 across all sites (arranged from the top left to the bottom based on increasing mean air temperature). The daily mean NEP is aggregated (mean ± 95 % CI as error bars) at 1 °C T_soil bins.

3.3 Effect of warming on NEP anomalies

At warmer sites (low latitude or altitude < 1000 m a.s.l.) where NEP showed significant changes in winter 2020 (IT-SR2, BE-Bra, DE-Obe), the average NEP anomaly increased by 75 %. Conversely, at colder sites where NEP significantly changed in winter 2020 (SE-Nor, CZ-BK1, RU-Fyo, FI-Let, IT-Ren, SE-Svb), the average NEP anomaly decreased by 8.8 %, indicating reduced net carbon uptake (Fig. 4). Overall, at sites where daily NEP differed significantly during the warm winter, NEP increased at only three sites, while it declined significantly at the remaining six sites (Fig. 4). Changes in NEP are attributed only to changes in climatic factors because the forests did not undergo significant changes in the canopy density, except at FI-Let. While FI-Let was affected by a partial cut in 2016 (Korkiakoski et al., 2019, 2020), winter fluxes remained relatively stable in all pre- and post-harvest years as the partial cut affected mostly the summer fluxes (data not shown here).

Figure 9 shows the sensitivity of NEP anomalies to anomalies of air temperature, soil temperature, and radiation across different sites. Overall, the sensitivity of NEP anomalies to soil temperature anomalies was larger than to anomalies of air temperature and radiation as shown by the test of the slope of change in NEP anomalies (Fig. 9).

Figure 8 Relative conditional variable importance (RCVI, %) of three climatic variables for explaining the variance in daily winter NEP, GPP, and R_eco and the overall variability explained (r²) (marked with red triangles) estimated from the random forest regression (RFR) analysis. The RFR model was trained on winter observations during the reference period (2014–2019). Sites are ordered by increasing site mean annual temperature (from left to right). For modeling R_eco, GPP was used as an additional predictor (see the methods section for more detail).

Figure 9 Sensitivity of NEP anomalies in winter (ΔNEP) to (a) anomalies of incoming solar radiation (ΔR_g), (b) anomalies of air temperature (ΔT_air), and (c) anomalies of soil temperature (ΔT_soil). The sensitivities represent the slope of ΔR_g, ΔT_air, and ΔT_soil when regressed with ΔNEP using a multivariate linear regression (ΔNEP ∼ ΔR_g + ΔT_air+ΔT_soil). The non-significant (p<0.05) sensitivity is shown as a transparent point. Error bar shows the 95 % CI of the slope obtained from the multivariate linear regression. Sites are ordered by increasing site mean air temperature (from left to right).

While the relationship between air temperature and soil temperature was stronger than the relationship between radiation and air temperature during the winter, we observed large variability in the coupling between soil and air temperature across the sites, as is shown in Table 3.

Table 3 Pearson correlation coefficient between mean daily incoming shortwave radiation (R_g), air temperature (T_air), and soil temperature at 5 m (T_soil) at each site during the reference period (2014–2019). Sites are ordered by decreasing mean air temperature. Leaf area index (LAI) values are shown as mean across the study period ± standard error of the mean.

4 Discussion

4.1 Warming of the air and the soil in winter

We tested how climate variables and CO₂ fluxes deviated from a reference period (2014–2019) during the warm winter of 2020 across 14 evergreen needleleaf forest sites distributed from the north to the south of Europe (from Sweden to Italy). The sites where winter 2020 was particularly warm and dry were not clustered in a certain climatic region; however we observed a consistent pattern that the warming of the air was more pronounced at the colder sites (p<0.05, $r=-\mathrm{0.53}$, Figs. 2, S11).

The strength of the coupling between the air and the soil temperature was not similar across all sites. In forests, topsoil temperature is directly affected by changes in air temperature; however, several underlying processes and properties modify the magnitude of decoupling between air and soil temperatures. This decoupling can reach up to 10 °C, depending on the season and the properties of the biome type (Lembrechts et al., 2022). These underlying factors and processes include for example (1) a vertically complex and horizontally continuous forest structure that leads to higher decoupling of the soil temperature from the air temperature; (2) soil moisture content as moisture increases the soil heat storage; (3) insulation by the litter or snow cover; (4) cloud cover, ground surface albedo, and the rate of evapotranspiration which collectively affect the radiation balance and energy exchange between the soil and the air; and (5) microtopography that affects the drainage of air (e.g., cool air drains in low-lying areas) (Guan et al., 2009; Lozano-Parra et al., 2018; De Frenne et al., 2021; Gril et al., 2023). Although the direct effect of canopy closure on snow distribution, accumulation, and melting in different periods was not tested here, it was evident that sites that had a larger LAI also showed a tighter coupling between air temperature and soil temperature (p<0.05, r=0.69, Table 3) as forest canopy structure influences the coupling of air and soil temperature in forest ecosystems, for example by shading the soil and reducing the snow depth beneath denser canopies (Woods et al., 2006; Gao et al., 2022). Snow cover plays a crucial role in regulating the coupling between soil and air temperatures due to its insulating properties (Liu et al., 2023). At the sites selected for this study, the coupling between air and soil temperatures also tended to weaken as snow depth increased, though the relationship was not statistically significant (r=0.46, p>0.05). A more detailed analysis across a larger number of sites is required to test if this may be due to uncertainties in remotely estimating snow depth beneath the forest canopy.

4.2 Winter warming effect on forest CO₂ fluxes

Winter warming had both positive and negative effects on the net ecosystem productivity (NEP) of coniferous forests. However, at most sites, NEP declined during the warm winter. This variation can largely be attributed to the differential effects of warming on soil versus air temperatures (Bond-Lamberty and Thomson, 2010), which influence both soil respiration and tree CO₂ uptake. When soil temperatures rise above freezing, root activity and respiration increase. Meanwhile, tree productivity may respond to higher air temperatures, potentially enhancing CO₂ uptake. However, air warming alone, without corresponding soil warming, might not disrupt dormancy (Bowling et al., 2024) if the soil within the rooting zone remains frozen. For instance, at IT-Ren, where daytime NEP significantly declined during the warm winter, air temperatures rose by over 3.5 °C above normal but soil temperatures stayed at freezing levels (Fig. S1).

Our hypothesis was that warming in winter will lead to a larger negative effect on net ecosystem productivity (i.e., higher CO₂ emissions) across colder forests. While we observed that (1) across most sites winter emissions increased during the warm winter and (2) generally emissions in winter increase in response to an increase in soil temperature (observed at seven sites, Fig. 7), we did not find a link between warming of the air and increased emissions that would confirm this general hypothesis.

CO₂ fluxes are influenced by both temperature and light (Figs. S5–S9). However, evidence suggests that the temperature response of biochemical processes is shaped more by the plant's growth temperature than by the immediate temperature (Fürstenau Togashi et al., 2018). Moreover, the net ecosystem productivity (NEP) response to similar temperatures can vary seasonally, exhibiting a clear hysteresis effect, and is further influenced by environmental factors such as solar radiation and soil moisture (Niu et al., 2011). Across different sites, the sensitivity of NEP to temperature generally increases as site mean temperatures decrease. Once temperature is no longer a limiting factor at higher site mean temperatures, radiation becomes the dominant constraint on NEP (Running et al., 2004; Fig. 9).

Chamber-based observations from boreal forests indicate that snow depth and soil moisture significantly influence the temperature sensitivity of soil CO₂ fluxes, particularly due to the effects of freeze–thaw cycles on soil moisture content (Du et al., 2013). Warmer winters can lead to increased respiration rates and greater nutrient availability for trees, as thawed soils exhibit a higher Q₁₀ than frozen soils, resulting in a more rapid increase in soil respiration in response to warming (Wang et al., 2014). However, this expected increase in respired CO₂ may be moderated by microbial carbon limitation, which can be alleviated through higher inputs of labile carbon from plant material and root exudates (Sullivan et al., 2020).

Furthermore, aboveground productivity tends to increase with rising temperatures (Figs. S6 and S7), contributing to enhanced autotrophic respiration. Warming during winter also affects the microbial communities responsible for decomposing labile and stable organic carbon in the soil, which can offset the respiration response to temperature and ultimately reduce overall soil respiration (Tian et al., 2021). The extent of belowground autotrophic respiration in response to warming, along with the supply of labile substrates through rhizodeposition and root exudates, plays a crucial role in determining net CO₂ fluxes under warming conditions (Nyberg et al., 2020). In our study, the sensitivity of ecosystem respiration (R_eco) to air temperature (Q₁₀) did not significantly change during the warm winter, remaining comparable to the Q₁₀ observed during the reference period (Fig. S10).

A decrease in snowpack and increased soil freezing can have immediate short-term impacts on plant CO₂ uptake, but these changes may also lead to lasting negative effects on tree function (Repo et al., 2021). Particularly, sites with prolonged cold winter seasons could experience significant adverse effects from winter warming. Trees in northern latitudes and at higher altitudes may be especially vulnerable to these changes, as their optimal growth temperatures are more sensitive to short-term temperature fluctuations. Conversely, in ecosystems with greater seasonal temperature variability, the optimal temperature range for growth is broader (Weng and Liao, 2010; Liu, 2020).

4.3 Winter tree physiology effect on CO₂ fluxes

Responses of coniferous species to soil warming can vary largely depending on the species' adaptive traits; the overall ecosystem context; and interactions with other environmental factors such as precipitation, temperature, and nutrient availability (Dawes et al., 2017; Oddi et al., 2022). The sites we studied here, although dominated by evergreen needleleaf species, consisted of different canopy species, and some sites were dominated by a mixture of species (Table 1). There can be significant differences in photosynthetic parameters across different species of evergreen conifers that would affect the tree and ecosystem response to warming (Fürstenau Togashi et al., 2018). The different responses of productivity to increased warming in ENFs can stem from differences in the quantity (and quality) of stored NSCs in the roots and the rate at which this C storage is mobilized within the tree during the warm winter (Bansal and Germino, 2009). Warmer temperatures and dry conditions in winter lead to stomatal closure and depletion of carbohydrate reserves for trees that are adapted to ample precipitation and conditions of low VPD (vapor pressure deficit) in winter, and this effect leads to reduced CO₂ uptake of trees during warmer winters (Earles et al., 2018).

Low temperatures are essential for signals that trigger the synthesis of soluble carbohydrates involved in osmotic and freezing protection against cold extremes (Chang et al., 2021) that otherwise impair the Calvin cycle by inhibiting the regeneration of ribulose bisphosphate (RuBP) and decrease the efficiency of ribulose-1,5-bisphosphate carboxylase/oxygenase (Rubisco) carboxylation (Ensminger et al., 2012; Crosatti et al., 2013). Non-structural carbohydrates (sugar and starch) that are accumulated during the growing season are utilized in winter to ensure the survival of trees (Zhu et al., 2012; Tixier et al., 2020), and failure to develop overwintering defenses can cause evergreen conifer needles to remain susceptible for example to photo-oxidative damage during frost events (Chang et al., 2016). Studies that combine ecosystem-scale flux measurements with tree-level observations have the potential to closely examine the adverse effects of winter warming on cold-adapted forests.

Our results represent the first analysis of the impact of winter warming on CO₂ fluxes in European evergreen needleleaf forests, highlighting the importance of understanding the various underlying mechanisms that influence forest CO₂ dynamics. Data on the responses of photosynthetic traits over ecologically relevant timescales (ranging from days to years) are limited. However, eddy covariance observations offer a valuable opportunity to construct long-term time series of canopy-level processes, enabling the investigation of the effects of extreme climatic conditions across all seasons. We also encourage further research that integrates long-term observations with plant-level experiments to explore how changes in winter functioning may influence trees' responses to early-season extremes, such as spring frost and drought, and to better understand the implications of these extremes for ecosystem carbon uptake.

5 Conclusions

Our study examined the effects of the warm winter of 2019–2020 on CO₂ fluxes in evergreen needleleaf forests across Europe. We showed that during wintertime net CO₂ emissions increased in response to warming across most sites. However, responses varied by location, with factors such as forest structure and local climatic conditions creating microclimates that either mitigated or intensified the impact of warming on CO₂ fluxes. By combining long-term eddy covariance data with plant-level experiments, we can gain valuable insights into how winter warming affects forest ecosystems. Future research should prioritize understanding the carryover effects of winter warming on tree responses to seasonal climatic extremes, as this knowledge is crucial for predicting how cold-adapted forests will respond to ongoing winter warming.