Predicted Future Changes in the Mean Seasonal Carbon Cycle Due to Climate Change

Abstract

Through photosynthesis, forests absorb annually large amounts of atmospheric CO₂. However, they also release CO₂ back through respiration. These two, opposite in sign, large fluxes determine how much of the carbon is stored or released back into the atmosphere. The mean seasonal cycle (MSC) is an interesting metric that associates phenology and carbon (C) partitioning/allocation analysis within forest stands. Here, we applied the 3D-CMCC-FEM model and analyzed its capability to represent the main C-fluxes, by validating the model against observed data, questioning if the sink/source mean seasonality is influenced under two scenarios of climate change, in five contrasting European forest sites. We found the model has, under current climate conditions, robust predictive abilities in estimating NEE. Model results also predict a consistent reduction in the forest’s capabilities to act as a C-sink under climate change and stand-aging at all sites. Such a reduction is predicted despite the number of annual days as a C-sink in evergreen forests increasing over the years, indicating a consistent downward trend. Similarly, deciduous forests, despite maintaining a relatively stable number of C-sink days throughout the year and over the century, show a reduction in their overall annual C-sink capacity. Overall, both types of forests at all sites show a consistent reduction in their future mitigating potential.

1. Introduction

Forests play a pivotal role in the biosphere–atmosphere feedback by annually absorbing large amounts of atmospheric CO₂ through photosynthesis (GPP; ~150 PgC year⁻¹) and releasing it back because of, e.g., ecosystem respiration (R_eco), a relatively close amount yet not necessarily equal, which varies year by year [1,2,3]. The net ecosystem exchange (NEE) of CO₂ between ecosystems and the atmosphere is the net balance between these two gross fluxes opposite in sign, and it governs much of the overall terrestrial annual net carbon (C) budget. Imbalances between CO₂ sources (even including carbon lost by fires and other processes) and sinks directly increase or decrease atmospheric CO₂ levels [4]. Terrestrial ecosystems—and forests in particular—are contributing substantially to climate change mitigation, provided that they are C-sinks and not C-sources [5]. Forests that might absorb more than they emit are commonly considered carbon sinks (with NEE-negative in sign), while if they emit more than they absorb are considered as carbon sources (with NEE-positive in sign). The magnitude of this exchange of CO₂, however, is subject to substantial variability and trends, in large part as a response to variations and trends in climate [6]. Indeed, forests have been shown to be extremely sensitive to changes in environmental conditions (e.g., climate, seasonality, atmospheric CO₂ concentration, nitrogen deposition), to aging [7,8], and to disturbances [9], including management practices [10,11], which can control both photosynthesis and respiration. Therefore, estimating NEE, GPP, and R_eco is a key step for better understanding the underlying mechanisms constraining ecosystem functioning [12].

Europe and the Mediterranean are expected to become in the near future a ‘Hot Spot’ of climate change [13,14,15,16]. The literature reports that under climate change scenarios forests are expected to grow faster, to mature earlier but also to die younger, curtailing their life span [17], because of, mainly, warming and increased atmospheric CO₂ concentration (the so-called ‘CO₂-fertilization effect’) [18]. Conversely, there is a general lack of evidence and knowledge on how, overall, forest ecosystems will, on the whole, react to climate change. Forest carbon balance will be impacted by climate change because various main processes are impacted, which, in turn, may react and respond differently to climate change, also because they are vulnerable and sensitive to separate environmental drivers. Indeed, while an extensive line of evidence shows that the increased availability of CO₂ may amplify the photosynthetic rate and assimilation capacity [19], such an increase is largely debated since there is no evidence that such positive changes will generally continue indefinitely [20,21]. Similarly, the effects of warming are largely discussed because, while it is documented that some species may take advantage of a longer vegetation season (e.g., deciduous species), there are also negative effects linked to warming such as heat waves and the often-associated drought events [22,23], including late frost [24] and disturbances [25], which can be detrimental to growth till tree survival. In addition, there is a general concern that the changing temperature response of respiration turns boreal forests from carbon sinks into carbon sources [26]. Indeed, warming has also been found to accelerate both autotrophic as well as heterotrophic respiration (the two components of ecosystem respiration), meaning that increased temperature may lead forests to release more carbon, potentially more than absorbed annually [27]. Conversely, drought has been shown to reduce microbial respiration and then heterotrophic respiration [28]. How these processes (i.e., photosynthesis and ecosystem respiration) will be impacted by climate change annually will determine much of the future forest annual C-budget.

The Mean Seasonal Cycle (MSC) metric, which reflects the average distribution of flux (i.e., NEE, GPP, and R_eco) throughout the days of a year, is an insightful measure linking phenology with carbon partitioning and allocation within seasonal climatic variability. By capturing the typical fluctuations in a specific region due to changing seasons, the MSC highlights the expected seasonal changes in climate data, averaged over many years, to smooth out anomalies and emphasize the regular, cyclical nature of these changes. Many studies [29,30,31,32] have indeed shown that climate change will impact both the phenology by changing the date for the beginning and the end of the growing season, as well as by changing the shape of the Leaf Area Index (LAI) distribution over the year, which, at the same time, will influence the way, among the other things, when photosynthesis can start and how recent and old photosynthates are partitioned and used to build new tissues and to replenish the reserves used for the metabolism, as well as carbon allocation [33] and the C-dynamic.

Process-based models are valuable tools to understand how and to what extent future climate change will impact these two fluxes (GPP and R_eco) in the MSC, both processes being controlled by warming and changes in precipitation regime and atmospheric CO₂ concentration [34,35]. Here, we at first applied and validated, under current observed climate conditions, the ‘Three Dimensional–Coupled Model Carbon Cycle–Forest Ecosystem Module’ (3D-CMCC-FEM), a biogeochemical, biophysical process-based forest ecosystem model designed to simulate carbon, nitrogen, and water cycle in forest ecosystems and, secondly, under climate change conditions.

Specifically, we question and analyze: (1) the capability of the 3D-CMCC-FEM to represent under the current climate the main C-fluxes governing the C-cycle in terms of net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_eco), by validating the model against independent data from the Fluxnet network; and, provided that the model is close to the observed data; (2) how, and if, the sink/source mean seasonality will be influenced under two locally bias-corrected scenarios of warming (RCP 2.6 and 6.0) and atmospheric CO₂ enrichment from three CMIP5 Earth System Models, within ISIMIP-PROFOUND initiative, in five well-studied and long-monitored contrasting forest sites (three evergreens and two deciduous) on a longitudinal transect through Europe up to the end of the century.

2. Materials and Methods

2.1. Model Description (3D-CMCC-FEM ‘v.5.6’)

The ‘Three Dimensional–Couple Model Carbon Cycle–Forest Ecosystem Module’ (hereafter ‘3D-CMCC-FEM’) is a biogeochemical, biophysical, process-based forest ecosystem model (see [10,11,36,37,38,39,40,41,42,43,44] and reference therein). The model is designed to simulate carbon, nitrogen, and water cycles in forest ecosystems at commonly 1-hectare spatial resolution and the main eco-physiological processes (e.g., photosynthesis) at daily temporal resolution. The most recent code versions since Collalti et al. [17] adopt the biogeochemical photosynthesis model of Farquhar, von Caemmerer, and Berry [45] to compute gross primary productivity (GPP). The biogeochemical photosynthesis model is parameterized as in Bernacchi et al. [46,47] and temperature acclimation for leaves as in Kattge and Knorr [48]. The 3D-CMCC-FEM considers, as in De Pury and Farquhar [49], light interception, reflection, transmission, and assimilation (and leaf respiration) for both sun and shaded leaves. Autotrophic respiration (R_A) is computed mechanistically following the ‘Growth and Maintenance Respiration Paradigm’ (GMRP; [50]), which is divided into the metabolic costs for synthesizing new tissues (growth respiration, R_G) and the metabolic costs for maintaining the existing ones (maintenance respiration, R_M). In 3-D-CMCC-FEM, the maintenance respiration is based on Nitrogen amount (a fixed fraction of carbon mass varying between the six tree compartments) and is temperature-controlled by a standard Arrhenius relationship [36]. ‘Type I’ and ‘Type II’ acclimation of respiration to temperature (i.e., short- and long-term acclimation; [10,39,51]) are also accounted for. Any imbalance between carbon assimilation and carbon losses because of plants’ respiration is buffered by a seventh pool, the Non-Structural Carbon pool (NSC; starch and sugars undistinguished), which has priority in the carbon allocation all over the year. The net primary production (NPP) is the GPP minus R_A. Biomass production (BP) is the amount of NPP not used for replenishing the NSC pool. Indeed, other forms of non-structural carbon losses (e.g., biogenic volatile organic compounds, BVOCs, or root exudates to mycorrhizas) are not accounted for by the model. The phenological scheme, as well as the carbon partitioning/allocation scheme, distinguished for deciduous and evergreen tree species, is in-depth described in Collalti et al. [10,36,39,52]. Heterotrophic respiration follows a BIOME-BGC-like approach (which follows the CENTURY model; [53,54]) distinguishing decomposition for litter and soil pools with each of the four different conceptual pools characterized by different decomposability degrees (i.e., fast, medium, slow, and a recalcitrant carbon pool) [39,55]. Altogether, litter and soil decomposition emissions form heterotrophic respiration (R_H), which, summed up to R_A, constitutes ecosystem respiration (R_eco). Net ecosystem exchange (NEE) is equal to R_eco minus GPP. Therefore, negative values indicate carbon uptake from the atmosphere (i.e., the system acts as a C-sink, NEE < 0), and positive values indicate carbon release (i.e., the system acts as a C-source, NEE > 0). The 3D-CMCC-FEM’s sensitivity to its 54 species-specific parameters and how it varies along the forest development and under climate change is described in Collalti et al. [17].

In the present study, we used version 5.6 [41], which slightly differs from v.5.5, described in Collalti et al. [8] and Dalmonech et al. [11], for a new scheme (and relative parameterization) for sapwood and live wood turnover and dynamics and some additional new forest management schemes (not used here) described in Testolin et al. [43].

2.2. Case Study Areas

Five case studies were selected as representative of the main European forest species (and climate) and at the same time because of long-monitored sites and part of the Fluxnet network [56], the ISIMIP (Inter-Sectoral Impact Model Intercomparison Project, [57]; https://www.isimip.org/, accessed on 1 January 2020) initiative and the PROFOUND database [42,58,59]: the temperate European beech (Fagus sylvatica L.) forest of Collelongo, Italy (IT-Col), and of Sorø, Denmark (DK-Sor), the maritime pine (Pinus pinaster Ait.) forest of Le Bray, France (FR-Lbr), the boreal Scots pine (Pinus sylvestris L.) forest of Hyytiälä, Finland (FI-Hyy), and the temperate Norway spruce (Picea abies (L.) H. Karst) forest of Bílý Křìž, Czech Republic (CZ-Bk1) (Figure 1). Stand characteristics are described in

Table 1. Site description and stand initialization data are used in simulations with data from the PROFOUND database [58]. Model initialization data (e.g., diameter at breast height—DBH, tree height, stand age, and stand density) correspond to the stand conditions of the first year in the historical simulations.

Fluxnet Code	Nation	Site Name	Coordinate			Climate	Species Composition	Period of Simulation and Climate Data		Mean DBH	Tree Height	Mean Stand Age	Stand Density
			Lat (°)	Lon (°)	Elevation (m)			Historical (Observed)	Historical + RCP (Modeled)	(cm)	(m)	(years)	(n. trees ha−1)
IT-Col	Italy	Collelongo	41.8	13.6	1560	Temperate	Fagus sylvatica	1997–2014	1997–2099	20.2	17.3	105	900
FR-LBr	France	Le Bray	44.7	–0.8	61	Temperate	Pinus pinaster	1996–2008	1997–2099	26.7	17.8	26	614
FI-Hyy	Finland	Hyytiälä	61.9	24.3	181	Boreal	Pinus sylvestris	1996–2014	1997–2099	13	11.3	35	870
CZ-Bk1	Czech Republic	Bílý Kříž	49.3	18.3	875	Temperate	Picea abies	2000–2008	1997–2099	9.1	7.5	19	2388
DK-Sor	Denmark	Sorø	55.3	11.4	40	Boreal	Pinus sylvestris	1996–2008	1997–2099	24.3	21	75	353

.

2.3. Input, Meteorological Data, and Climate Change Scenarios

To run, 3D-CMCC-FEM needs a set of input data from state variables representing the stand at the beginning of simulations and that account for structural characteristics (e.g., tree height, average age, diameter at breast height; see Table 1) as well as carbon and nitrogen pools (e.g., stem carbon and nitrogen); meteorological forcing data (e.g., daily maximum and minimum temperature, daily precipitation); annual atmospheric CO₂ concentrations; and species-specific parameters (e.g., maximum stomatal conductance). All the data used to initialize the model in the present study for the five stands come from the ISIMIP initiative and the PROFOUND database [58]. More specifically, daily observed meteorological data for model validation come from the Fluxnet2015 Dataset [56], while daily modeled historical (1997–2005) and future climate scenarios (2006–2099) are those from the ‘ISIMIP 2bLBC’ experiments (‘2b experiments Locally Bias Corrected’) coming from three different Earth System Models (ESMs; GFDL-ESM2M, IPSL-CM5A-LR, and MIROC5, respectively) based on the Climate Model Intercomparison Project 5 (CMIP5) driven by two Representative Concentration Pathways (RCPs) of atmospheric greenhouse gas concentration trajectories, namely, RCP 2.6 and RCP 6.0. [61,62]. The ISIMIP 2bLBC have the same structure as those in the 2b experiments but have been corrected by improving the method described in Hempel et al. [63] and subsequently by the methods described in Frieler et al. [57] and Lange [64] using the observed meteorology at the local level [58]. Therefore, the 2bLBC climate data represent the more consistent and closer modeled climate data with the observational data. The annual atmospheric CO₂ concentrations for the historical period are based on Meinshausen et al. [65] and have been extended for the period from 2006 to 2015 with data from Dlugokencky and Tans [66]. Values specific for each RCP for the period 2016 to 2099 are also based on Meinshausen et al. [65] and were used within the Farquhar, von Caemmerer, and Berry [45] photosynthesis model with values varying, at the end of this century, from 421.4 μmol mol⁻¹ (RCP 2.6) to 666.4 μmol mol⁻¹ (RCP 6.0), respectively. NEE, GPP, and R_eco data, for model validation, come from the Fluxnet2015 Dataset [56]. Other variables have been validated at these forest stands in the past (although using slightly different model versions) and described in Collalti et al. [8,10,37], Marconi et al. [40], Mahnken et al. [42], and Dalmonech et al. [11].

2.4. Model Runs, Validation, and Mean Seasonal Cycle Analyses under Climate Change

The model simulations for model validation under measured forcing climate ran for CZ-Bk1 from 2000 to 2008, for IT-Col from 1997 to 2014, for FI-Hyy from 1996 to 2014, for FR-Lbr from 1996 to 2008, and for DK-Sor from 1996 to 2008 (see Table 1). For all sites model simulations under climate change scenarios began in 1997 and finished in December 2099. Model validation was performed by comparing modeled NEE, GPP, and R_eco against measured eddy covariance estimates (for GPP and R_eco using the night-time method with constant USTAR [67], as reported in the Fluxnet2015 Dataset [56].

To analyze 3D-CMCC-FEM’s capabilities to simulate NEE, GPP, and R_eco for daily and monthly time series, a set of commonly used statistical metrics have been applied to compare measured vs. modeled data (under observed climate forcing). To avoid considering bad quality data, a filtering procedure for quality-check has been applied; days with less than 60% of valid data were not considered and excluded both in the model and in the observed datasets. Therefore, daily NEE, GPP, and R_eco eddy covariance data with low-quality values (i.e., less than 0.6; [68]) were removed. Consequently, the corresponding daily modeled data were removed too. The monthly NEE, GPP, and R_eco values (both from eddy covariance and the model) were aggregated from the daily ones. The common statistic we applied includes Pearson’s correlation coefficient (r), Relative Mean Bias (RMB), Normalized Root Mean Square Error (NRMSE), and Modeling Efficiency (ME)

In climate change projections, we considered the potential modifications in the ability of forest stands to absorb or emit carbon throughout the season and across the years under two different locally bias-corrected climate change scenarios, each coming from three Earth System Models. This involves averaging the daily values of the MSC of the three fluxes considered. It is noteworthy that NEE represents the equilibrium between carbon absorption by vegetation during photosynthesis and carbon release through vegetation and microbial respiration. Not only the length of the growing season but also the balance between the yearly amount of photosynthesis and R_eco has been shown to control much of the variability across the sites and the decades analyzed [47,49]. This calculation is derived from the variance between GPP and R_eco encompassing both autotrophic respiration (R_A, including ground components) and heterotrophic respiration (R_H) [69].

Under climate change scenarios, we account for potential changes in the sink/source capacity of the stands during the season by averaging NEE, GPP, and R_eco values every ten years up to 2100 [70] and accounting for changes in the sink/source and source/sink length during the year, computed as the number of total days of the year (DoY) where a forest stand behaves as a C-sink (GPP > R_eco with NEE < 0) or a C-source (GPP < R_eco with NEE > 0), as described by the NEE. In addition, we also account for the changes in the DoY where NEE changes its sign at least for ten consecutive days to avoid artifact effects of pulsing, e.g., the ‘Birch effect’ on R_eco [71], and to account for unstable conditions and no clear source/sink and sink/source seasonal transition during the year. Therefore, we discuss changes in the MSC under the RCP 2.6 and RCP 6.0 for NEE, being the net result of opposite fluxes (i.e., GPP and R_eco), through its changes in negative values, i.e., days in the year where NEE < 0 and describing CO₂ uptake from the atmosphere, and positive, i.e., days in the year where NEE > 0 and describing CO₂ release to the atmosphere. In this way, we account for potential changes that may lead to anticipations or delays in the switch from source/sink and sink/source capacity, which often happens during spring and autumn during the year. The analyses under climate change scenarios (2006–2099) also include the changes in the shape of the Mean Seasonal Cycle (MSC) for NEE, GPP, and R_eco values and for the changes (both in the absolute and the percentage values) in the annual value. Changes in MSC were estimated on the ten-year average values using the 1998–2008 decade as a benchmarking reference. Furthermore, we account for changes in the annual values of NEE, GPP, and R_eco due to climate change.

3. Results

3.1. Model Validation vs. Fluxnet Data

The NEE, GPP, and R_eco, as modeled by 3D-CMCC-FEM, exhibit strong correlations with the observed daily and monthly eddy covariance data and for the overall MSC across all five sites (Figure 2 and Figure A1 and

Table A1. Summary of the statistics between simulated and measured NEE from the Fluxnet2015 Dataset (Pastorello et al., 2020 [56]), calculated on the 5 cases studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-LBr, Hyytiälä—FI-Hyy, Bílý Křìž—CZ-Bk1, and Sorø—DK-Sor). The table shows the daily and monthly values for Person’s Coefficient (r—dimensionless), Relative Mean Bias (RMB—dimensionless), Normalized Root Mean Square Error (NRMSE—dimensionless), Modeling Efficiency (ME—dimensionless), and Mean Absolute Error (MAE, gC m⁻² time⁻¹).

	IT-Col	FR-LBr	FI-Hyy	CZ-Bk1	DK-Sor
Daily NEE
r	0.96	0.93	0.93	0.85	0.97
RMB	0.07	−0.86	−0.48	0.14	−0.46
NRMSE	−0.57	−1.56	−1.54	1.39	−6.59
ME	0.92	−0.10	0.58	0.69	0.46
MAE	0.67	1.20	0.67	1.15	1.42
Monthly NEE
r	0.98	0.95	0.93	0.92	0.99
RMB	0.07	−0.88	−0.49	0.21	−0.46
NRMSE	0.99	−1.55	−1.52	−0.50	−6.14
ME	0.95	−0.15	0.55	0.79	0.55
MAE	11.42	33.39	19.63	19.12	37.78

,

Table A2. Summary of the statistics between simulated and measured GPP from the Fluxnet2015 Dataset (Pastorello et al., 2020 [56]), calculated on the 5 cases studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-Lbr, Hyytiälä—FI-Hyy, Bílý Křìž—CZ-Bk1, and Sorø—DK-Sor). The table shows the daily and monthly values for Person’s Coefficient (r—dimensionless), Relative Mean Bias (RMB—dimensionless), Normalized Root Mean Square Error (NRMSE—dimensionless), Modeling Efficiency (ME—dimensionless), and Mean Absolute Error (MAE, gC m⁻² time⁻¹).

	IT-Col	FR-LBr	FI-Hyy	CZ-Bk1	DK-Sor
Daily GPP
r	0.99	0.95	0.93	0.95	0.99
RMB	0.15	0.55	0.01	−0.14	0.00
NRMSE	0.34	0.42	0.39	1.21	0.19
ME	0.91	0.28	0.86	0.85	0.97
MAE	0.91	1.51	0.91	0.97	0.75
Monthly GPP
r	0.99	0.96	0.93	0.99	1.00
RMB	0.15	0.52	0.01	−0.06	0.00
NRMSE	0.31	0.40	0.38	0.15	0.16
ME	0.92	0.34	0.86	0.97	0.97
MAE	20.68	41.69	26.05	14.25	20.26

and

Table A3. Summary of the statistics between simulated and measured R_eco from the Fluxnet2015 Dataset (Pastorello et al., 2020 [56]), calculated on the 5 cases studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-LBr, Hyytiälä—FI-Hyy, Bílý Křìž—CZ-Bk1, and Sorø—DK-Sor). The table shows the daily and monthly values for Person’s Coefficient (r—dimensionless), Relative Mean Bias (RMB—dimensionless), Normalized Root Mean Square Error (NRMSE—dimensionless), Modeling Efficiency (ME—dimensionless), and Mean Absolute Error (MAE, gC m⁻² time⁻¹).

	IT-Col	FR-LBr	FI-Hyy	CZ-Bk1	DK-Sor
Daily Reco
r	0.90	0.89	0.94	0.86	0.89
RMB	0.82	0.16	−0.30	0.30	−0.42
NRMSE	0.72	0.24	0.41	1.16	0.35
ME	−1.17	0.50	0.75	0.63	0.62
MAE	1.14	0.71	0.74	0.88	1.29
Monthly Reco
r	0.99	0.94	0.95	0.99	0.95
RMB	0.77	0.14	−0.29	0.19	−0.42
NRMSE	0.83	0.20	0.39	0.24	0.29
ME	−0.64	0.68	0.77	0.93	0.71
MAE	25.42	18.01	21.02	11.61	31.82

in Appendix A). Some slight differences are observed for the daily values for temperate European beech forests of IT-Col and DK-Sor when representing the NEE and R_eco between the 100 and 200 DoY (Figure 2(1.a–1.c) and Figure 2(5.a–5.c)). Modeled overestimations for GPP of about 5 gC m⁻² day⁻¹ at the peak of production (~180 DoY) for both the European beech forest at IT-Col (Figure 2(1.b)) and the maritime pine forest at FR-Lbr (Figure 2(2.b)) was observed. The highest correlation coefficients between modeled and observed data were observed for DK-Sor and IT-Col (r = 0.97 and 0.96, respectively) for the daily NEE and GPP (Table A1, Appendix A), while FI-Hyy shows the best correlation for daily R_eco (r = 0.94) and monthly NEE (r = 0.99) (Table A2 and Table A3).

At the IT-Col and CZ-Bk1, the model shows the best performances in simulating NEE (RMB = 0.07 and 0.14, respectively—Table A1), while the lowest values were reached by DK-Sor and IT-Col for the GPP (RMB ranging between −0.06 and 0.55, Table A2). Overall, RMB values for GPP are relatively low across all case studies and time scales. Regarding R_eco at FI-Hyy and DK-Sor, the model tended to underestimate both daily (RMB of −0.3 and −0.42, respectively) and monthly (−0.29 and −0.42, respectively). Conversely, at IT-Col, the model exhibited a slight overestimation (RMB daily 0.82 and monthly 0.77, Table A3).

The model reports negative NRMSE values for NEE across all time scales, indicating a slight overestimation (values ranging from −0.57 to −6.59) with the exception of the CZ-Bk1 site for the daily NEE value (1.39), and at the IT-Col site for the monthly NEE value (0.99). Regarding the GPP, (Table A2) DK-Sor showed the lowest NRMSE for daily and monthly values (0.19 and 0.16, respectively). Conversely, at CZ-Bk1 and FR-Lbr forests, the model displays the highest NRMSE for daily (1.21) and monthly (0.40) values, respectively. Last, Table A3 shows the validation results for R_eco. At FR-LBR, the highest accuracy and precision were reported with the lowest NRMSE for daily and monthly values (0.24 and 0.20 respectively). For the other sites, we found almost the same model capability described for the above fluxes (i.e., NEE and GPP) with a lower correlation and a slightly higher error in terms of daily R_eco at CZ-Bk1 (NRMSE = 1.16).

ME exhibits values close to one across all time scales and sites for NEE, with values ranging from −0.10 (FR-LBr) to 0.92 (IT-Col) for daily NEE values and −0.15 (FR-LBr) to 0.95 (IT-Col) for monthly values. The modeled GPP has a similar trend of the observed NEE with values ranging from 0.28 for FR-LBr to 0.97 for DK-Sor for daily values and FR-LBr = 0.34 to DK-Sor = 0.92 for monthly values (Table A2). The R_eco achieved lower performances than GPP and NEE in terms of ME; Table A3 displays values from 1.17 for IT-Col to 0.75 for FI-Hyy and IT-Col = −0.64 to FI-Hyy = 0.77, for daily and monthly time scale, respectively.

Finally, the lowest MAE was found at IT-Col, for NEE (0.67 gC m⁻² day⁻¹ and 11.42 gC m⁻² month⁻¹) and GPP (0.91 gC m⁻² day⁻¹ and 20.68 gC m⁻² month⁻¹), while the lowest MAE for R_eco was for FR-LBr (0.71 gC m⁻² day⁻¹ and 18.01 gC m⁻² month⁻¹) (Table A1, Table A2 and Table A3 in Appendix A).

3.2. Mean Seasonal NEE Cycle under Climate Change Scenario

Figure 3 and Figure 4 display the 10-year average NEE seasonal cycle under climate scenarios (i.e., RCP 2.6 and RCP 6.0) for the five case studies. Overall, across all the sites and scenarios considered, there is a consistent reduction in the absolute NEE over time (i.e., NEE is ‘less negative’ showing a reduction in the sink capacity) with changes in the source/sink (NEE becomes negative and the site turns into a C-sink) and sink/source (NEE becomes positive and the site turns into a C-source) switch over the year(s). However, RCP 2.6 generally exhibits lower reductions in annual and mean seasonal NEE when compared to RCP 6.0 across most study sites and time intervals. The loss in the modeled sink capacity is because R_eco increases more than GPP, and the differences between the scenarios are related to an increase in R_eco higher than that for GPP under RCP 6.0.

The rate of decrease in the annual NEE varies among locations and forest species. For example, the beech forests at IT-Col and DK-Sor show a reduction in the NEE (i.e., the site is less C-sink) of about 68%, and 43% (Figure 3 and Figure 4—panel 1 and panel 5) at the end of the century, generally exhibiting a more moderate decrease compared to evergreen sites. The Scots pine forest shows, at the end of the century, a reduction in the sink capacity of 95%, standing out as the most significant decrease in NEE (Figure 3 and Figure 4—panel 3). Over time, the rate of decline shows a tendency to speed up, hinting at a reduction in the capacity of the carbon sink because of, despite an overall increase in the GPP (Figure A2 and Figure A3), the increased ecosystem respiration due to increased temperatures (Figure A4 and Figure A5). By the end of the century, substantial reductions in NEE, i.e., the sites become less C-sink, across all locations and scenarios were simulated, with some locations experiencing over 70–90% reduction compared to the 1999–2009 decade (

Table A4. Summary of the NEE variation (%) from the first decade of simulation (1999–2009), considering both climate change scenarios (RCP 2.6 and 6.0) for 5 case studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-LBr, Hyytiälä—FI-Hyy, Bílý Kříž—CZ-Bk1, and Sorø—DK-Sor). In bold values where changes were the highest between the decades while underlined the lowest ones. Note that negative values indicate that NEE becomes less negative, e.g., −100% indicates a reduction in the negative values of NEE.

NEE
RCP 2.6
year	CZ-Bk1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	−37%	−19%	−4%	−24%	−9%
2019–2029	−64%	−35%	−16%	−44%	−14%
2029–2039	−70%	−31%	−27%	−53%	−22%
2039–2049	−76%	−39%	−39%	−59%	−29%
2049–2059	−71%	−48%	−47%	−58%	−33%
2059–2069	−79%	−59%	−58%	−65%	−32%
2069–2079	−75%	−61%	−66%	−70%	−37%
2079–2089	−76%	−69%	−76%	−75%	−44%
2089–2099	−74%	−68%	−95%	−72%	−43%
RCP 6.0
year	CZ-Bk1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	−35%	−11%	0%	−25%	−11%
2019–2029	−60%	−27%	−11%	−43%	−14%
2029–2039	−74%	−36%	−23%	−52%	−26%
2039–2049	−76%	−43%	−36%	−57%	−34%
2049–2059	−75%	−44%	−46%	−61%	−31%
2059–2069	−75%	−43%	−57%	−64%	−23%
2069–2079	−72%	−54%	−69%	−69%	−30%
2079–2089	−76%	−59%	−79%	−75%	−34%
2089–2099	−71%	−60%	−87%	−68%	−36%

).

The GPP, similar to NEE, generally increases across forest types and scenarios over time, albeit with varying degrees (Figure A2 and Figure A3). Specifically, RCP 6.0 shows higher growth rates, particularly in later years (2059–2099), with increases ranging from 28 to 51% across all studied forests, compared to 10–26% for RCP 2.6 (

Table A5. Summary of the GPP variation (%) from the first decade of simulation (1999–2009), considering both climate change scenarios (RCP 2.6 and 6.0) for 5 case studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-LBr, Hyytiälä—FI-Hyy, Bílý Kříž—CZ-Bk1, and Sorø—DK-Sor). In bold values where changes were the highest between the decades while underlined the lowest ones.

GPP
RCP 2.6
year	CZ-Bk1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	7%	6%	22%	12%	8%
2019–2029	10%	5%	34%	15%	15%
2029–2039	16%	11%	39%	18%	18%
2039–2049	17%	13%	42%	19%	19%
2049–2059	18%	12%	43%	20%	20%
2059–2069	18%	11%	41%	18%	23%
2069–2079	17%	12%	36%	18%	18%
2079–2089	16%	9%	32%	18%	15%
2089–2099	17%	10%	26%	16%	15%
RCP 6.0
year	CZ-BK1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	7%	7%	23%	11%	7%
2019–2029	10%	9%	36%	14%	13%
2029–2039	13%	10%	43%	18%	14%
2039–2049	15%	12%	44%	20%	15%
2049–2059	20%	17%	49%	20%	20%
2059–2069	23%	23%	50%	23%	28%
2069–2079	27%	26%	51%	26%	31%
2079–2089	28%	29%	50%	26%	31%
2089–2099	31%	30%	51%	28%	34%

). Analysis of long-term trends suggests saturation and subsequent slight decreases in GPP growth rates under RCP 2.6, with this phenomenon being most noticeable in the forests of FI-Hyy. RCP 6.0 consistently shows a higher percentage of change in R_eco compared to RCP 2.6. By the end of the century (2089–2099), R_eco changes from 37 to 106% for RCP 2.6 and 60 to 121% for RCP 6.0 across all forests examined (

Table A6. Summary of the R_eco changes (%) from the first decade of simulation (1999–2009), considering both climate change scenarios (RCP 2.6 and 6.0) for 5 case studies selected (i.e., Collelongo—IT-Col, Le Bray—FR-LBr, Hyytiälä—FI-Hyy, Bílý Kříž—CZ-Bk1, and Sorø—DK-Sor). In bold values where changes were the highest between the decades while underlined the lowest ones.

Reco
RCP 2.6
year	CZ-Bk1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	29%	16%	39%	47%	15%
2019–2029	47%	22%	66%	71%	26%
2029–2039	59%	29%	83%	87%	33%
2039–2049	64%	35%	96%	94%	36%
2049–2059	63%	38%	102%	95%	40%
2059–2069	67%	40%	107%	99%	43%
2069–2079	63%	44%	104%	103%	39%
2079–2089	62%	43%	104%	106%	37%
2089–2099	63%	44%	106%	101%	37%
RCP 6.0
year	CZ-Bk1	IT-Col	FI-Hyy	FR-LBr	DK-Sor
2009–2019	27%	14%	39%	45%	14%
2019–2029	45%	24%	66%	69%	23%
2029–2039	56%	29%	87%	85%	28%
2039–2049	60%	36%	97%	93%	33%
2049–2059	67%	43%	112%	99%	39%
2059–2069	71%	51%	121%	108%	47%
2069–2079	76%	60%	130%	117%	54%
2079–2089	80%	66%	136%	123%	55%
2089–2099	82%	67%	142%	121%	60%

). Boreal and maritime pine forests experienced higher changes compared to deciduous forests. For example, by 2089–2099, under RCP 6.0, changes reach 142% and 121% at FI-Hyy and FR-Lbr forests, while deciduous forests like European beech in IT-Col and DK-Sor experience lower changes at 67% and 60% respectively.

3.3. Changes in NEE Dynamics under Different Climate Scenarios

We considered fluctuations in the length of sink/source and source/sink forest stand behaviors throughout the year. This calculation involved determining the total number of days in a year (‘N. days year⁻¹’) where a forest stand exhibited either C-sink (NEE < 0) or C-source behavior (NEE > 0). The 10-year average number of days as a C-sink and as a C-source under the RCP 2.6 and 6.0 for the five case studies selected is presented in Figure 5a. The number of days identified as a C-sink in the evergreen forests for the scenario RCP 2.6 (i.e., FR-LBr, CZ-Bk1, and FI-Hyy) starts relatively low in the first decade (1999–2009) but increases significantly over time, showing a consistent upward trend (Figure 5(2.a,3.a,4.a). At the CZ-Bk1 site, the number of days considered as a C-sink starts at 47 N. days year⁻¹ in the 1999–2009 decade and increases steadily over the decades, peaking at 199 N. days year⁻¹ in 2059–2069 and remaining relatively high thereafter. For FI-Hyy forests, the number of days as a C-sink starts at 107 N. days year⁻¹ in 1999–2009 and rises to 217 days in 2089–2099, with a peak of 225 in 2079–2089. Lastly, at the forests of FR-LBr, there were no days as a C-sink in 1999–2009, but they increased steadily over time, reaching 114 days by 2079–2089. The number of days as a C-source varies inversely to those as a C-sink across decades, demonstrating a general reduction over time for evergreen forests. The number of days functioning as a C-source decreased, from ~300 in 1999–2009 to ~200 N. days year⁻¹ in 2089–2099 (Table A1), for evergreen forests. Deciduous forests on the contrary almost maintain a relatively stable number of days as a C-sink through the century (Figure 5(1.a,5.a). The number of days as a C-sink at IT-Col forests ranges from 212 to 224 N. days year⁻¹, while at DK-Sor from 214 to 222 N. days year⁻¹, with slight fluctuations observed across decades and no clear overall trend. The same trend for the capacity of the stand to act as a C-source, but with slightly different days in the range from 137 to 159 N. days year⁻¹ recorded in the forest of Sorø, and from 140 to 149 N. days year⁻¹ for the IT-Col site. The discrepancies between the already described RCP 2.6 and RCP 6.0 scenarios are minimal as they exhibit a similar trend for evergreen and deciduous, respectively, with only a slight change in the number of days, mainly in the last decade, with a magnitude of ~10 N. days year⁻¹.

To assess the changes in MSC due to climate change scenarios, we focused on the shifting patterns of the DoY over decades, particularly examining when forest stands transitioned from a C-source to a C-sink and vice versa (see Figure 5b). Regardless of the scenario considered, the European beech forests revealed constant transition periods across decades. In IT-Col, the shift from being a C-source to a C-sink occurred between DoY 136 and 150, with the opposite transition from sink to source around DoY 285. Meanwhile, in DK-Sor, the same transitions happened between DoY 115 and 127, and reversed between DoY 272 and 280 (Figure 5(1.b,5.b)). Even for evergreen forests, there are no noticeable differences in the shift corresponding to RCP 2.6 and RCP 6.0 scenarios. At CK-BZ1, the shift to a C-sink occurs relatively early in the year, spanning from DoY 2 to 31 across the decades. At FI-Hyy, the transition timings vary widely, ranging from DoY 30 to 50 over the decades. Similarly, at FR-LBr, the C-sink transitions occur from DoY 1 to 16 across the decades, with some decades exhibiting earlier shifts. On the other hand, transitions to C-sources at CZ-BK1 occur from late November to early December (DoY 201–329), displaying a decreasing trend over the years. In contrast, at FI-Hyy, transitions take place from mid-September to late September (DoY 197–285), while at FR-LBr, they occur from late September to early October (DoY 255–285) (Figure 5(2.a,3.a,4.a)).

4. Discussion

4.1. Model Validation

To sum up statistical metrics, the 3D-CMCC-FEM performs best in replicating the mean seasonal patterns of the three fluxes in European beech forests DK-Sor and IT-Col. R_eco also reaches a satisfactory performance for the boreal Scots pine forest of FI-Hyy, and the Scots pine forest of CZ-Bk1. The robust predictive ability of 3D-CMCC-FEM in estimating NEE across different timeframes, forest species, and climates, as proved by its alignment with the eddy covariance Fluxnet2015 Dataset [56], underscores its effectiveness in capturing the complex dynamics of carbon fluxes within forest ecosystems, as documented by previous works (see, e.g., [10,11,36,40,42]). Nevertheless, there are slight inconsistencies, particularly during peak photosynthesis periods (around the 200th day of the year) in both deciduous and evergreen forests, which have been acknowledged in the existing literature. Studies indicate that estimates of ecosystem respiration derived from eddy covariance often underestimate actual values for forest ecosystems [72,73,74,75].

4.2. Mean Seasonal NEE Cycle under Climate Change Scenario

The sensitivity of forest ecosystems to changes in environmental factors such as climate change, seasonal variations, and atmospheric CO₂ levels has been thoroughly evidenced in the literature. Indeed, the interaction among these variables plays a crucial role in sha** the carbon exchange dynamics [6,12,17,31,76].

The analysis of the mean NEE seasonal cycle under different climate change scenarios presents several fascinating findings. The first key observation is the consistent reduction in forest stand capabilities to act as carbon sinks from the atmosphere in the coming years, across all study sites and climate change scenarios. Despite the overall lack of understanding regarding how forest ecosystems will respond to climate change, several new studies concur regarding the diminishing forest carbon uptake capabilities, thus confirming our statement [26,31,77,78,79].

In broader terms, climate change affects forest carbon balance by influencing key processes, which can respond differently due to their sensitivity to various environmental drivers [80]. Indeed, the reduction in NEE, within the forests observed in this study, has a disparity in the intensity of climate impacts across these scenarios: RCP 2.6 exhibits less intense NEE reduction compared to RCP 6.0 scenarios. The warming accelerates both autotrophic and heterotrophic respiration meaning that increased temperature may lead forests to an increase in R_eco [27,58,81,82]. On the other hand, the increasing atmospheric CO₂ concentration intensifies the GPP through the ‘carbon fertilization effect’ (i.e., reported to be the cause of 44% of the GPP increase since the 2000s) [18,19]. Under the RCP 2.6 scenario, R_eco exhibits a steady linear increase until the century’s end, while GPP follows a bell-shaped curve, reaching saturation around the mid-century. Consequently, as GPP saturates, its compensatory capacity reduces, while R_eco continues to rise due to further temperature increases in the latter half of the century resulting in a further decline in NEE. Conversely, for the RCP 6.0 scenario, GPP saturation occurs only towards the end of the century, and while one might anticipate a reversal in NEE trends considering this factor, it does not materialize. This is because, under the RCP 6.0 scenario, a higher temperature increase is predicted compared to RCP 2.6 (i.e., [58]), leading to a more pronounced rise in ecosystem respiration relative to photosynthesis. As a consequence, there is a greater decrease in the stand forest’s carbon sink capacity. The result refers to all forests studied in this work, but in the forest of FI-Hyy, the phenomena look more pronounced (Figure A2 and Figure A3 in Appendix A).

The fate of ecosystem carbon flux depends not only on atmospheric and climate conditions but also on the type of forests analyzed. The reduction in C-sink capabilities is particularly notable in evergreen forests, which exhibit a higher decrease in NEE compared to evergreen sites. The boreal Scots pine forest of FI-Hyy stands out with the most significant reduction in NEE, indicating a heightened vulnerability to climate change effects. This finding is supported by the study of Hadden and Grelle [26] who found that, over a 17-year period, the forest ecosystem in a boreal forest stand in northern Sweden transitioned from being a carbon sink to a carbon source. This could mean that past efforts to validate the neutrality hypothesis [83] with climate change impact show limitations, and we need new research directions and new perspectives to better capture changes in the carbon fluxes between the ecosystem and atmosphere [84]. Indeed, a long-standing debate around the C-neutrality of old-growth forests (and some of the sites become old-growth at the end of simulations) raises concerns, and debates will increase about the reduction of the sink capacity of aging forests, as assumed by Odum’s theory. As found in [7] (but see [85]), we also found that even a >200-year-old stand (as IT-Col in 2099) still has sink capacity. The annual NEE decrease (and much less the changes in the MSC), as shown by the model results, is certainly a function not only of climate but also of the inherent effects related to, e.g., biomass, both live and dead, accumulation (which led to increases in respiratory costs), and changes in, e.g., the forest structure (which led to decreases in the carbon assimilation), as stands age. Nevertheless, such effects on annual NEE are expected to be greatly exacerbated by climate change.

4.3. Changes in NEE Dynamics under Different Climate Scenarios

Climate change impacts plant phenology by altering the start and end dates of the growing season, which influences when photosynthesis can begin and consequently affects C-fluxes [20,29,30,32,86].

A primary finding is that, regardless of the scenario analyzed, the number of days identified as a C-sink in evergreen forests increases significantly over time, indicating a consistent upward trend. Similarly, the number of days classified as a C-source decreases over the decades, showing a general reduction. The second finding is that for evergreen forests, the DoY to C-sink tends to increase (indicating a forward shift in the year when the system becomes a sink), and the DoY to C-source decreases (indicating a backward shift in the year when the system becomes a source), aligning with the overall trend of fewer C-sink days and more C-source days over time. In contrast, deciduous forests maintain a relatively stable number of C-sink (and C-source) days throughout the century, reflecting a steady DoY when the system becomes a C-sink (or a C-source), despite an anticipated beginning of the growing season but compensated by higher respiration rates. Indeed, for the deciduous, the 3D-CMCC-FEM simulates the bud-brake through a thermic sum function and leaf and fine root development (and the relative growth respiration peak) in a well-defined and short time during spring [36]. Conversely, for the evergreens, leaf and fine root growth development is spread all over the spring. At the same time, leaf fall in the deciduous starts under certain hours of solar radiation, and thus, this is not under the control of climate, while in the evergreen, it happens all over the year, and under the control of climate, balanced by incoming photosynthates for new leaves and fine roots. Ultimately, deciduous spring C-sink capacity is counterbalanced by high C-emissions mainly because of growth respiration in spring. Such behavior is different for evergreens, which lengthen their C-sink capacity during spring. However, the lengthening of the growing season does not automatically mirror an increase in the net sink capacity because R_eco shows an increase much more than GPP. It is generally acknowledged that the changing temperature response of respiration transforms forests from C-sinks to C-sources [26,76,87,88,89], while the stability of the carbon sink/source dynamics over the decades for deciduous ecosystems is a relatively recent finding. The DoY for deciduous forests remains unchanged because the earlier start of the growing season, triggered by rising temperatures, is balanced by an earlier increase in respiration. This compensates for the earlier rise in GPP at the level of NEE. Overall, the lack of change in the number of C-sink (and C-source) days across decades and the reduction in the NEE suggest that, over the long term, deciduous forests are more efficient in using photosynthates compared to evergreen forests [90,91].

5. Limitations

The modeling framework presented has certain limitations that must be acknowledged. First, we deliberately decided to not simulate the effects of anthropogenic disturbances, e.g., forest management, nor the ones from natural disturbances caused by climate change, such as windstorms, forest fires, and insect outbreaks, so as to concentrate on the effects of climate change alone and avoid these potentially confounding effects. Climate extreme events are presumed to be incorporated into the climate scenarios used to drive the model and are, therefore, somewhat already accounted for. Moreover, indirect changes due to climate change in key factors like nitrogen deposition, phosphorus, or ozone—which could potentially amplify or mitigate our findings—were not evaluated. Nonetheless, some research (e.g., [92]) indicates that this issue might not be universally relevant. These studies highlight the strong response of various tree species to CO₂ fertilization across different levels of nutrient availability. Lastly, the potential for species migration to and from the study areas was not considered. However, such dynamics might require longer timescales than those covered in this study and it is unlikely (although still possible in theory) that species composition may completely change throughout the simulations.

6. Conclusions

The MSC metric is an interesting and descriptive metric that associates phenology and carbon partitioning allocation within forest stands. Climate change impacts both the phenology, by changing the date for the beginning and the end of the growing season, and the ecosystem carbon allocation.

We applied the process-based forest model 3D-CMCC-FEM to evaluate the potential modifications on the ability of different forest stands to absorb or emit carbon throughout the season and across the years up to 2100. Before that, we validated the model under current climate conditions and found a robust predictive ability of 3D-CMCC-FEM in estimating NEE, GPP, and R_eco across different timeframes, forest species, and climates.

The analysis of the mean NEE seasonal cycle under different climate change scenarios presents a consistent reduction in the forest stand capabilities to act as a carbon sink in the coming years, across all study sites, and climate change scenarios. The reduction in NEE ability has different intensities of climate impacts across these scenarios. The RCP 2.6 scenario demonstrates a less pronounced reduction in NEE compared to the RCP 6.0 scenario. This disparity primarily stems from variations in key variables, such as the differing rates of temperature increase between the two scenarios, as well as the CO₂ fertilization effect, while, in all sites, age effects depend on the age at the beginning of the simulations. The reduction in C-sink capabilities is mainly notable in evergreen forests, which exhibit a higher decrease in NEE compared to deciduous forest sites.

Finally, we found that the number of days as a C-sink in evergreen forests increases over the years, indicating a consistent upward trend. Conversely, the number of days as a C-source decreases over the decades, showing a general reduction. This statement aligns with the forward shift of DoY to C-sink, and the backward shift of DoY to C-source. In contrast, deciduous forests maintain a relatively stable number of C-sink (and C-source) days throughout the century, reflecting a fixed DoY when the system becomes a C-sink (or a C-source). The DoY for deciduous forests remains constant as the earlier onset of the growing season, driven by warming temperatures, is offset by an earlier uptick in respiration. Decades pass with little change in the number of days as a C-sink (and a C-source), alongside a decrease in NEE. This indicates that deciduous forests, over the long haul, demonstrate greater efficiency in utilizing photosynthates when compared to evergreen forests.