Environmental Effects on Normalized Gross Primary Productivity in Beech and Norway Spruce Forests

: The strong effects of climate change are expected to negatively impact the long-term resilience and function of forest ecosystems, which could lead to changes in forest carbon balance and productivity. However, these forest responses may vary with local conditions and forest types. Accordingly, this study was carried out to determine gross primary productivity (GPP) sensitivity to changes in environmental parameters. Central European beech (at Štítná) and spruce species (at Bílý K˘rí˘z and Rájec), growing under contrasting climatic conditions, were studied. The comparative analyses of GPP were based on a ﬁve-year-long dataset of eddy covariance ﬂuxes during the main growing season (2012–2016). Results of forest GPP responses with changes in environmental factors from a traditional Stepwise multiple linear regression model (SMLR) were used and compared with Random forest (RF) analyses. To demonstrate how actual GPP trends compare to potential GPP (GPP pot ) courses expected under near-optimal environmental conditions, we computed normalized GPP (GPP norm ) with values between 0 and 1 as the ratio of the estimated daily sum of GPP to GPP pot . The study conﬁrmed the well-known effect of total intensity of the photosynthetically active radiation and its diffuse fraction on GPP norm across all the forest types. However, the study also showed the secondary effects of other environmental variables on forest productivity depending on the species and local climatic conditions. The reduction in forest productivity at the beech forest in Štítná was presumed to be mainly induced by edaphic drought (anisohydric behaviour). In contrast, reduced forest productivity at the spruce forest sites was presumably induced by both meteorological and hydrological drought events, especially at the moderately dry climate in Rájec. Overall, our analyses call for more studies on forest productivity across different forest types and contrasting climatic conditions, as this productivity is strongly dependent on species type and site-speciﬁc environmental conditions.


Introduction
Gross primary productivity (GPP) constitutes an essential part of net ecosystem CO 2 exchange (NEE) between the atmosphere and the forest ecosystems [1][2][3]. GPP provides (dry spruce forest at Rájec), and CZ-Stn (beech forest at Štítná). CZ-BK1 is also a candidate for ICOS (Integrated Carbon Observation System; https://www.icos-cp.eu/, accessed on 6 March 2021) network. All forest stands are even-aged monocultures. Their main characteristics are presented in (Table 1). The reference evapotranspiration amount at each forest stand was computed from in situ data and was additionally used for quantifying the atmospheric evaporative demand from 2012 to 2016. References [28] [29] [23] * 2012-2016 study period.

Eddy Covariance and Ancillary Measurements
At each station, the EC system consisted of an infrared gas analyzer (LI-COR, Lincoln, NE, USA) and an ultrasonic anemometer (Gill Instruments, Lymington, UK) measuring at 20 Hz frequency (detailed description in Table A1). Each system was mounted on a tower at height above the forest canopy (as specified in Table A1). EC measurements were complemented by an extensive set of sensors to collect the required auxiliary meteorological data.
Measurements included air temperature (T air ) and relative humidity at the top of the forest canopy with the EMS33 temperature and humidity sensors (Embedded Moisture Sensor, Vancouver, BC, Canada). Hourly precipitation (P) was determined using a Precipitation Gauge 386C (Met One Instruments, Grants Pass, OR, USA). Measurements for the incoming photosynthetic active radiation (PAR) were made with a LI-190R Quantum Sensor (LI-COR, NE, USA) at both Bílý Kȓíz and Štítná and with an EMS12 sensor (EMS, CZ) at Rájec. Additionally, profiles of soil moisture were measured using the CS616 (Campbell Scientific, North Logan, UT, USA) sensors at both spruce forest sites (in Bílý Kȓíz and Rájec) and with the ThetaProbe (ML2x 355, Delta-T, Burwell, UK) sensors at Štítná. An overall description of the instrumentation at the study sites including soil temperature measurements (T s ) is given in Table A2.
To analyze the effects of diffuse and direct radiation on daily mean GPP values, the dataset was divided into sunny and cloudy days based on the clearness index (CNI). CNI is defined as the ratio of solar irradiation transmitted through the atmosphere onto the Earth's surface relative to extraterrestrial irradiation. Sunny days were classified as days with CNI > 0.7 during daytime hours, cloudy days were characterized with CNI < 0.4, and CNI values of 0.4-0.7 were grouped as partly cloudy days [30].
The vapour pressure deficit (VPD; hPa) for each site was computed from T air ( • C) and relative air humidity (RH; %) according to [31]. The daily sum of PAR (MJ m −2 day −1 ) was derived following the work in [32].

Soil Water Content Simulations
As the soil volumetric water content (SVWC) was not measured over the entire period of the study (the measurements started at the beginning of 2016), the simulated daily SVWC was used. The simulations were performed by soil water balance model R-4ET (R package for Empirical Estimate of Ecosystem EvapoTranspiration, [33]). This soil water balance model was calibrated using the Bayesian statistics implemented in the R package BayesianTools [34] with Differential-Evolution Markov chain Monte Carlo sampler. In calibration using the Bayesian statistics, the model input parameters were updated iteratively to provide a probability distribution of the calibrated parameters representing the uncertainty in the measured data and structure of the model. These simulations were repeated several times (4.8 × 10 6 ), with the first 1.8 × 10 6 simulations based on the prior distribution and the remaining 3 × 10 6 runs restricted by the posterior distribution that resulted from the first set of processes. This high number of iterations was significant in attaining a good input parameter convergence with narrow distribution. The Gelman diagnostics based on a criterion for potential scale reduction factor (less than 1.2 for all parameters) as used in [35,36] was used to inspect the convergence. Additionally, a final selection of parameters from their probability distributions was made using a maximum posteriori probability estimate.
Moreover, the input variables for the soil water balance model comprised meteorological data and leaf area index (see in [33] for more details). The soil at each of the sites was then stratified into 12 layers, increasing in thickness with depth to~3 m. These soil layers were selected in a manner that matched the depths of all sensors with an extra ±2.5 cm due to the volume measured by these sensors. At all sites, the available SVWC measurements were available throughout the main root zone region. At the wet spruce forest site in Bílý Kȓíz, the CS616 (Campbell Scientific, Inc., Logan, UT, USA) sensors were placed and measured at depths of 0.05, 0.1, 0.22, 0.34, and 0.42 m. Furthermore, at the dry spruce forest site in Rájec, the same type of sensors was placed and measured at depths of 0.05, 0.1, 0.2, 0.5, and 0.8 m. At the beech forest site in Štítná, the ThetaProbe (ML2x 355, Delta-T, UK) sensors were placed and measured at depths of 0.05, 0.1, 0.3, 0.6, and 0.9 m. The optimized model parameters included soil parameters such as the SVWC at saturation, field capacity, wilting point, and saturated hydraulic conductivity-all optimized at all 12 depths [33]. There were also additional single parameters of relevance to the SVWC simulations, which were optimized. These included the rooting depth with the Beta parameter, which described the root profile shape [37], surface resistance and the degree of isohydricity [38], the water interception capacity of leaf and bark area, and curve number representing a runoff parameter [33]. A uniform distribution of priors was applied if the lower and upper limits were set to be within ±50% of the values based on field measurements of wilting point, field capacity, and saturated water content and ±100% for the remaining parameters that were estimated from literature or previous anecdotal analysis. In providing the model with sufficient spin-up time to stabilize and provide reliable and robust parameterization, the model simulations were initiated at the start of 2010 with initial conditions of SVWC set to the field capacity estimated from soil texture (note that 2010 was one of the wettest years across all forest sites according to a computed standardized precipitation-evapotranspiration index over the region). Overall, the entire simulation was conducted from 2010-2019, with the observed data for the Bayesian calibration spanning from 2016 to 2019. The simulated SVWC averaged over all the depths produced a root mean square error of 0.037 m 3 m −3 at the wet spruce forest site in Bílý Kȓíz, 0.017 m 3 m −3 for the dry spruce forest site in Rájec and 0.032 m 3 m −3 for the beech forest site in Štítná, suggesting realistic SVWC estimates.

Turbulent Flux Measurements
Across all the three study forests, the EC data from the main growing season (May-September) of 2012-2016 were used for the analyses. Processing of EC data (wind components, sonic temperature, CO 2 , and water vapour mixing ratios) was performed using an open-source software EddyPro (Li-COR, Lincoln, NE, USA). The most recent methods for flux corrections, conversions, and thorough quality control scheme [39,40] were applied. This process involves raw data despiking and statistical screening, basic quality checking of turbulent fluxes (i.e., flux stationarity and integral turbulence characteristics tests), coordinate rotation using the planar fit method [41], spectral correction [42][43][44], detecting and compensating for time lags of signals from the ultrasonic anemometer and the gas analyzer [43], footprint estimation and calculating half-hourly final fluxes. Flux measurements over periods with insufficiently developed turbulence, i.e., low friction velocity (u*) were detected and filtered out. This filtering procedure, according to the most current methodology [45], assured the exclusion of CO 2 fluxes not necessarily representative of the ecosystem-scale biotic flux due to insufficient mixing across the canopy [46][47][48].
The R software [49] package 'REddyProc' [45] was used to gap-fill the EC data at all sites, using marginal distribution sampling [47]. CO 2 flux (subsequently considered to be equal to NEE) was partitioned into GPP and ecosystem respiration (R eco ). The flux partitioning approach by [47] using daytime data was applied to estimate half-hourly GPP ((µmol m −2 s −1 )) values. The half-hourly GPP values were aggregated to obtain daily and monthly sums of GPP values.

Estimation of Potential and Normalized GPP
GPP pot is defined here as the estimate of daily GPP sum that would be attainable under near-optimal environmental conditions (PAR, CNI, SVWC, VPD, T air , T s , and P) at given site and day of the year (DOY [50]). The GPP pot thus forms the boundary line of a scatter plot of GPP against DOY (especially for the growing season periods) pooled over multiple years of data. Specifically, to estimate the GPP pot , the following procedures were applied until there were no outliers as described in [50]: • Compute the 95th percentile from the daily sum of GPP pooled over years 2012-2016 for each DOY over a 7-day window (applied iteratively); • Outliers were detected and removed based on the percentiles method (all observations that positioned outside the interval formed by the 1 and 99 percentiles were considered as outliers).
To better compare the variation in GPP response across all the investigated forest stands, daily normalized GPP (GPP norm ) was derived for each forest station as the ratio of the estimated GPP to the GPP pot . Therefore, GPP norm values~1 represent days with maximum assimilation rates, whereas GPP norm values~0 represent days with extreme adverse effects on forest productivity. A smoothing spline curve was applied to depict the main trends in seasonal courses of GPP norm during the growing seasons of 2012-2016.
A Pearson's correlation coefficient matrix was calculated to determine the statistical relationship between GPP norm and the environmental variables based on a covariance method at a significance level of 0.05 ( Figure A2).

Multi-Linear and Tree-Based Regression Model Analyses
In order to assess the response of GPP norm to the environmental factors (PAR, CNI, SVWC, VPD, T air , and T s ) at each site, two methods were tested: (i) a stepwise multi-linear regression (SMLR) and (ii) random forest algorithm (RF). P was excluded from the analyses as it is not a direct measure of soil moisture; thus, SVWC was used for the regression analysis. The SMLR selection with interaction terms was designed using the stepwise regression method (in both forward and backward direction), to determine the significant terms in the model and eliminate the nuisance (i.e., non-significant) variables with 95% probability. In addition, quadratic terms were included to test for nonlinearity of the environmental variables with GPP norm . All statistical analyses were performed using the R software package 'stepAIC' (Version 7.3-54, [51]) for computing stepwise regression in both forward and backward direction [49].
We trained and applied RF to predict GPP norm at each of the forest ecosystems and for benchmarking the performance of the SMLR model. We used the R software package 'randomForest' (Version 4.6-14, [52]), with the following parameters established: number of trees of the model (ntrees) = 300, number of variables in each mode (nodesize) = 5, and the number of variables used in each tree (mtry) = one third of the total number of samples, as in [53] for regression RF. As part of the byproduct of the RF, the built-in variable importance measure that ranks the environmental variables (i.e., the features) according to their relevance in predicting GPP norm at each of the forest ecosystems has been provided [54]. These feature scores are computed from permuting data points that were not included in the RF analyses (Out-Of-Bag samples) based on the mean decrease in accuracy. The higher the value of the mean decrease in accuracy, the higher the importance of the variable in the model.

Accuracy Test of Regression Models
To assess the performance of both the SMLR model and RF techniques in predicting GPP norm at each of the forest ecosystems, model prediction accuracy indicators such as the Pearson correlation, the percentage of the variance explained (R 2 ) in relation to estimated and predicted GPP norm , and root mean square error (RMSE) values were derived and compared.

Variation in Meteorological Conditions at the Experimental Stations
Seasonal sum of precipitation during the main growing season period of the studied years at the wet spruce forest site in Bílý Kȓíz was approximately 40% and 37% higher than the seasonal sum of precipitation values at the dry spruce and the beech forest sites in Rájec and Štítná, respectively (Table 1 and Figure A1). Furthermore, the months of July-August were characterized mainly by hot and dry conditions across all forest stations, especially at the dry spruce forest in Rájec with high mean monthly T air (10% higher than at the wet spruce forest site and 6% lower than at the beech forest) and low mean monthly SVWC values (14% and 52% lower than at the wet spruce and beech forest sites respectively) ( Figure 1). However, the highest T air values were consistently observed at the beech forest in Štítná and the lowest at the spruce forest in Bílý Kȓíz. There were also statistically significant differences (p < 0.01) in the mean monthly VPD and SVWC values between the months of June-September in Bílý Kȓíz (with comparatively low monthly VPD values; Figure 1c) as compared to that in both Rájec and Štítná. These results indicate a drier climate at the dry spruce forest site in Rájec as compared to a humid climate in Bílý Kȓíz.

Effect of Meteorological Conditions on GPP
The annual course of mean daily sums of estimated GPP values in all the forests studied are shown in Figure 2. The area under the red curve represents the seasonal course of GPP pot , indicating how much the forest ecosystem could potentially assimilate carbon over the year under near-optimal environmental conditions. The systematic changes in estimated GPP during the growing season representing phenological responses to warming and other microclimatic conditions were well captured in Figure  Our results from the Pearson correlation matrix revealed a moderate positive linear relationship between GPP norm and PAR across all forest stands ( Figure A2). However, CNI only showed a moderate positive linear relationship with GPP norm in only the dry spruce and beech forest with a nonlinear relationship with GPP norm at the wet spruce forest ( Figures A2 and A3). Additionally, we found statistically significant correlation coefficients (p < 0.05) between GPP norm and other environmental variables such as PAR, P, VPD, SVWC, and T air at the wet spruce forest. Furthermore, at the dry spruce forest, there were also statistically signifcant correlation coefficients between GPP norm and PAR, CNI, SVWC, P, and T s . Statistically significant correlation coefficients (p < 0.05) between GPP norm and PAR, CNI, VPD, SVWC, T air , and P were also observed at the beech forest site.

GPP norm Prediction through Stepwise Multi-Linear Regression (SMLR)
To identify the main environmental variables influencing GPP norm across the beech and both spruce forest sites, an SMLR model was built and evaluated ( Figure 4). The accuracy indicators of the SMLR showed a good model quality in predicting significant environmental variables that influenced GPP norm values ( Table 2). The model over predicted for low GPP norm values and under predicted for high GPP norm values for all the three forest sites.
The Pearson correlation for all forest stations shows strong linear relationship between the predicted and estimated GPP norm values (Table 2). Furthermore, the R 2 of the SMLR model used across each forest ecosystem revealed that approximately 40-49% of the variance in GPP norm was well predicted by the environmental variables used (PAR, CNI, SVWC, VPD, T air , and T s ). The prediction error RMSE of the interaction model was found to be 0.18, 0.14 and 0.15 for the wet spruce, dry spruce and beech forests, respectively. This shows that model performance was comparable for all sites.
At both spruce forest sites, similar statistically significant (p < 0.01) environmental variables affecting GPP norm were PAR, SVWC, VPD, and T s (Tables A3 and A4). In addition to these similar environmental variables, CNI was also found to have secondary effects on GPP norm values at the wet spruce forest site in Bílý Kȓíz (Table A3). At the beech forest site in Štítná, all the environmental variables used in the SMLR analyses were found to be statistically significant (p < 0.01) in influencing GPP norm values, especially for PAR, T s , and CNI (Table A5).
Moreover, similar quadratic terms such as PAR 2 , T s 2 , and CNI 2 with negative coefficients were found to be statistically significant (p < 0.01) across both spruce forest sites, indicating a decrease in GPP norm values (downward sloping) with an increase in PAR, T s and CNI values (Tables A3-A5). Furthermore, at the beech forest site in Štítná, VPD 2 and CNI 2 proved to be inversely related to GPP norm and statistically significant (Table A5). However, the positive coefficient value for SVWC 2 across all the forest sites showed an increase in GPP norm values with increasing SVWC values (Tables A3-A5).   For the statistically significant (p < 0.01) interaction terms of SMLR, the VPD:SVWC term was common across all the forest sites showing the interactive effects of VPD and soil water availability on GPP norm values. At the wet spruce forest site in Bílý Kȓíz, interactive effects between VPD and other environmental variables such as SVWC, PAR, and T air were statistically significant. At the dry spruce forest site in Rájec, statistically significant (p < 0.01) interactive terms between SVWC and other environmental variables such as CNI, PAR, and VPD were observed. Furthermore, at the beech forest site in Štítná, statistically significant (p < 0.01) interactive terms such as T s :SVWC, T air :SVWC, T air :PAR, CNI:PAR, VPD:T s , and VPD:SVWC were all found to have effects on GPP norm values.

GPP norm Prediction through Random Forest Analyses (RF)
The validation statistics of the RF modeling are shown in Table 2. Note that higher R 2 (>0.54) and Pearson correlation values (>0.70) with low RMSE (<0.15) were obtained in GPP norm predictions with RF than with SMLR. This suggests that the RF model presents a greater potential for predicting changes in GPP norm with the environmental factors across the forest ecosystems studied ( Figure 5).

Importance of Environmental Variables in Random Forest Analyses
All the forest ecosystems predicted PAR as the most important environmental variable ( Figure 6). The variable importance score of PAR was higher at the beech forest site in Štítná than at both spruce forest sites. Additionally, at the wet spruce forest site in Bílý Kȓíz, the second and third most important environmental variables in GPP norm prediction were VPD and SVWC. On the other hand, SVWC and VPD were also the second and third most significant environmental variables, respectively, in GPP norm prediction at the beech forest site in Štítná. This shows the significant effect of VPD at the wet spruce forest site and soil water availability at the beech forest site on GPP norm prediction. Furthermore, at the dry spruce forest site in Rájec, the second and third most important variables in GPP norm prediction were VPD and T s , respectively. Overall, VPD had more significant effect on GPP norm prediction at both spruce forest sites, especially in Bílý Kȓíz (with the highest variable importance score of 0.18 for VPD).  As already observed from both the SMLR and RF analyses, PAR is the most important environmental variable in separating GPP norm values into two distinct groups (high and low values) across all forest ecosystems (Figure 7). For both spruce forest sites, the group with mostly higher GPP norm values (right branches in Figure 7) were separated above similar PAR threshold values (approximately 4 MJ m −2 day −1 ), while high GPP norm values were attainable at the beech forest in Štítná above a PAR threshold of approximately 4 MJ m −2 day −1 . The left branch in Figure 7 likely represents exceptional cases with low PAR and GPP norm values. Such conditions are the most common in the wet spruce forest in Bílý Kȓíz with 72 cases as compared to 38 cases in the dry spruce forest and 62 cases in the beech forest. Lower GPP norm values were realized during periods with extremely low PAR (~1 MJ m −2 day −1 ) values at the wet spruce forest in Bílý Kȓíz (6 cases) and the beech forest site in Štítná (7 cases

Discussion
This study used two sets of regression models (SMLR and RF) to identify the main environmental variables influencing the ratio between actual and optimal gross primary productivity (termed as GPP norm ) across central European beech (at Štítná) and Norway spruce (at Bílý Kȓíz and Rájec) species growing under contrasting climatic conditions. The presented results illustrate that the RF regression model outperformed SMLR and was highly effective for GPP norm prediction at all the studied forest ecosystems (Table 2). Though these statistical analyses may not show the mechanistic cause of the observed patterns, the RF model also provided useful information about the variable importance ( Figure 6) and the effect of their significant interactions (that are highly correlated; Figure 7) on GPP norm . Similarly, the authors of [17] used a path analysis approach to evaluate the sensitivity of GPP to environmental drivers during the summer drought of 2018, but our study evaluated in broader terms the general limiting environmental conditions of GPP across each of the studied forest sites.
As PAR is already a well-known primary driver of forest GPP, the results from our study revealed its overall significance in influencing GPP norm across all the forest sites. However, aside from PAR, there were other site-specific environmental variables that also had significant impact on GPP norm prediction. It is mostly hard to distinguish the direct impact and importance of T air , VPD, and SVWC as they are mutually dependent [55][56][57]. However, our results show that the second most important environmental variable in GPP norm prediction was VPD at the wet spruce forest site and SVWC at the beech forest site. Thus, indicating the contrasting water use strategies of both spruce and beech species. This is in line with earlier findings which show that the isohydric Norway spruce species has tighter stomatal control in response to increasing VPD than the anisohydric European beech species, which can tolerate drought stress better [29,[58][59][60]. Therefore, as more frequent heatwaves and drought spells characterize the proceeding climate change, an unavoidable shift of the Central Europe silviculture strategy is expected, providing likely higher opportunity for growing European beech while a higher risk for the coniferous Norway spruce [61,62]. Additionally, GPP norm reduction within the beech forest site was presumably mainly induced by edaphic drought stress (low soil moisture) conditions, further stressing the significant impact of SVWC on forest GPP, especially at the beech forest stand due to its different soil properties (more clayey soil) from the other forest sites [23,[63][64][65][66][67].
Moreover, at both spruce forest sites, a more detailed analysis of the variables splitting branches within RF (Figure 7) showed that T air or T s typically grouped subspaces with high and low (drought-reduced) GPP norm values. Interestingly, at the wet spruce forest site in Bílý Kȓíz, T air split the subspace with lower values of GPP norm at high T air threshold values above 19 • C during days with lower SVWC values (<0.16 m 3 m −3 ). This suggests that high T air values limit the forest productivity at the mountainous wet spruce forest site in Bílý Kȓíz. Furthermore, at the dry spruce forest in Rájec, the three most important variable (VPD, T s , and SVWC) had the highest scores among all the forest sites (scores > 0.1; Figure 6). This shows that the spruce forest productivity within a drier climate is often significantly impacted by both meteorological and hydrological drought events (decrease in atmospheric water content) [68][69][70][71][72]. Thus, under warmer climatic conditions, the impact of drought stress conditions on forest productivity will be more severe at spruce forests sites situated within drier climates [73].
Comparatively across all forest sites, the variable importance score of PAR was higher at the beech forest site in Štítná than at both spruce forest sites in Bílý Kȓíz and Rájec ( Figure 6). Although incoming PAR with a larger fraction of diffuse light was previously shown to be more effectively utilized by the forest canopy [19,74], generally lower scores of CNI suggest limited modulation of GPP norm by sky conditions. Relative to other assessed variables, CNI played the most important role at the beech forest site in Štítná. At both the beech and dry spruce forest sites, GPP norm was observed to increase with CNI until it became saturated at some point ( Figures A1 and A2). Thus, high GPP norm values are realized under partly cloudy conditions, especially within the beech forest ecosystems, possibly due to the deep penetration of light into its canopy due to its west-southwest slope orientation [75]. On the contrary, within the wet spruce forest site in Bílý Kȓíz with a mountainous ridge, high GPP norm values were mainly realized under cloudy conditions due to the effective penetration of anisotropic diffuse radiation at all layers within the forest canopy, especially in shaded leaves ( Figures A3 and A4) [19,76].

Conclusions
In this study, we identified the main environmental variables influencing the ratio between actual and optimal gross primary productivity across Central European beech (at Štítná) and Norway spruce (at Bílý Kȓíz and Rájec) species by comparing statistical analyses from the traditional Stepwise multiple linear regression model (SMLR) and Random forest (RF) regressional model. Across all the studied forest ecosystems, anisotropic diffuse radiation entering the forest canopy under cloudy conditions within the wet spruce forest and partly conditions within the dry spruce and beech forest was the main limiting environmental factor of photosynthesis. However, other environmental factors relating to water availability, such as the vapour pressure deficit, soil moisture, and temperature (air and soil), also had significant effects on the ecosystem photosynthesis depending on the local conditions and forest type. The sensitivity of both spruce forest sites to the vapour pressure deficit and high temperatures, especially in the drier climate at Rájec, indicate that the Norway spruce is currently one of the most threatened commercial tree species in Central Europe under the recent changes in climatic conditions. Moreover, reduction in forest productivity at the beech forest was presumed to be mainly induced by edaphic drought than meteorological and hydrological drought events, as in the spruce forest sites. Our findings provide useful insights potentially applicable for land surface models and further assessment of the impacts of climate change on forests in the Central Europe and their sequestration capacity. At the same time, our findings call for more studies on forest productivity across different forest types and contrasting climatic conditions, as productivity is strongly dependent on species type and site-specific environmental conditions.