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

Abstract

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ý Kr̆íz̆ and Rájec), growing under contrasting climatic conditions, were studied. The comparative analyses of GPP were based on a five-year-long dataset of eddy covariance fluxes 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 confirmed 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-specific environmental conditions.

1. Introduction

Gross primary productivity (GPP) constitutes an essential part of net ecosystem CO₂ exchange (NEE) between the atmosphere and the forest ecosystems [1,2,3]. GPP provides important information about the photosynthetic capacity of the ecosystem and is modulated by numerous environmental factors leading to changes in forest growth [4,5,6,7]. Nonetheless, the accurate quantification of forest GPP on a global scale remains challenging due to the spatio-temporal heterogeneity of forest ecosystems [8]. At the plot scale (within the area of the flux-tower footprint), the eddy covariance (EC) technique has emerged as a direct state-of-the-art approach to study the seasonality of carbon fluxes between the atmosphere and forest ecosystems, and thus represents an indirect yet reliable measure of GPP [9,10,11,12].

Seasonal trends of GPP demonstrate the forest response to variations in environmental factors. However, these forest responses vary from site to site depending on the species, age, local climate, and additional factors. For instance, rising temperatures within the optimum range and water availability in most temperate forest ecosystems is positive to reduce the freezing stresses and increase ecosystem productivity by enhancing the rate of ecosystem photosynthesis. However, higher temperatures (beyond the optimum level) may induce varying physiological and biochemical responses to heat stress at different forest ecosystems based on the acclimation processes of specific plant species to changes in their environment [13,14]. Such contrasting GPP responses (due to the shift in plant functional traits) of specific species to different environmental stress conditions at different sites are poorly understood [15]. Therefore, through ecosystem carbon flux measurements, physiological responses of individual plant species to certain climatic anomalies such as heat and water stress can be studied [16,17]. Furthermore, as multiple climatic factors affect forest carbon uptake, there is the need to identify near-optimal environmental conditions for different species, under which GPP can reach maximal theoretically attainable values termed potential GPP (GPP_pot).

GPP is primarily driven by available photosynthetically active radiation (PAR) coupled with other physiological processes like stomatal conductance that control carbon fluxes between the ecosystem and the atmosphere [18,19]. Moreover, the physiological processes that affect carbon uptake are also affected by complex interactions between a set of environmental drivers, including the quality of incoming solar radiation, temperature, atmospheric humidity, and soil moisture [18,20]. For example, the work in [19] has shown that photosynthesis of the upper spruce canopy could be severely depressed during clear sunny days due to the abrupt closure of stomata under high values of vapour pressure deficit. However, under cloudy conditions, the rate of photosynthesis of the lower spruce canopy increases due to the effective penetration of anisotropic diffuse radiation into the canopy, thus leading to an enhanced carbon gain and light use efficiency [21].

The main goal of this study is to determine the main environmental variables influencing the ratio between observed GPP and GPP_pot of three forest ecosystems with contrasting plant species using two types of regression models. Forests dominated by Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.) were selected for this study due to their sensitivity to climate and high commercial and ecological importance in Central Europe [22,23]. Furthermore, the current distribution of spruce species at often unsuitable locations with sub-optimal climate conditions make them vulnerable to both biotic and abiotic stresses [24,25,26,27]. Thus, it is important to study the physiological responses of different plant species to the local environmental conditions across different forest sites. EC data were collected from two Norway spruce forest sites (Bílý Kr̆íz̆ and Rájec) with distinct climatic conditions, and a beech forest (Štítná) from 2012–2016 during the main growing season (May–September) to exclude the impact of phenology.

2. Materials and Methods

2.1. Sites Description

This study uses multi-year (2012–2016) EC measurements from three forest ecosystem stations located in the Czech Republic (Central Europe) that are part of the CzeCOS (Czech Carbon Observation System; http://www.czecos.cz/, accessed on 6 March 2021) and FLUXNET Network (Flux Tower Network; https://fluxnet.fluxdata.org/, accessed on 6 March 2021). FLUXNET site IDs are CZ-BK1 (wet spruce forest at Bílý Kr̆íz̆), CZ-RAJ (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. Characteristics of the investigated sites.

Site Name	CZ-BK1	CZ-RAJ	CZ-Stn
Location	Moravian-Silesian Beskids Mountains	Drahany Highland	White Carpathian Mountains
Coordinates	49${}^{\circ }$30${}^{\prime }$08${}^{\prime }$${}^{\prime }$ N, 18${}^{\circ }$32${}^{\prime }$13${}^{\prime }$${}^{\prime }$ E	49${}^{\circ }$26${}^{\prime }$37${}^{\prime }$${}^{\prime }$ N, 16${}^{\circ }$41${}^{\prime }$48${}^{\prime }$${}^{\prime }$ E	49${}^{\circ }$02${}^{\prime }$09${}^{\prime }$${}^{\prime }$ N, 17${}^{\circ }$58${}^{\prime }$12${}^{\prime }$${}^{\prime }$ E
Elevation (in m a.s.l)	875	625	540
Topography	Mountainous (13${}^{\circ }$ slope with SSW exposure, located close to a mountain ridge)	Hilly (5${}^{\circ }$ slope with NNE exposure)	Mountainous (10${}^{\circ }$ slope with WSW exposure)
Ecosystem Type	Coniferous evergreen forest	Coniferous evergreen forest	Deciduous broadleaf forest
Prevailing species	Norway spruce (Picea abies (L.) Karst.)	Norway spruce (Picea abies (L.) Karst.)	European beech (Fagus sylvatica L.)
Canopy height (m)	16 (mean, as of 2015)	33 (mean, as of 2015)	31 (mean, as of 2015)
Stand age (years)	35 (as of 2016)	113 (as of 2016)	115 (as of 2016)
Mean seasonal air temperature (May-September; ${}^{\circ }$C)	14 *	16 *	17 *
Total seasonal precipitation (May-September; mm)	2730 *	1635 *	1719 *
Seasonal sum of reference evapotranspiration (May-September; mm)	2036 *	2325 *	2166 *
Soil type	Haplic and Entic Podzol	Modal Cambisol oligotrophic	Eutric Cambisol
References	[28]	[29]	[23]

* 2012–2016 study period.

). 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.

2.2. 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. Description of the eddy covariance systems at the investigated sites.

Site Name	CZ-BK1	CZ-RAJ	CZ-Stn
Ultrasonic Anemometer
Instrument	Gill HS-50 (Gill Instruments, UK)	Gill R3-100 (Gill Instruments, UK) but later changed to Gill HS-50 on 5 June 2015	Gill R3-100 (Gill Instruments, UK)
Gas Analyser
Instrument	LI-7200	Initially LI-7000 (IRG-0226) closed-path gas analyser, but later changed to LI-7200 on 5 June 2015	LI-7000 closed-path gas analyser
Measurement Height for the Eddy covaiance set-up (m)	Initially 20.5 m, but later changed to 25 m on 7 June 2016	41 m	44 m

). 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ý Kr̆í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ý Kr̆í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. Description of ancillary climatic measurements at the investigated sites.

Site Name	CZ-BK1	CZ-RAJ	CZ-Stn
Air Temperature and Humidity Profile
Instrument	EMS33 temperature and humidity sensor	EMS33 temperature and humidity sensor	EMS33 temperature and humidity sensor
Height (m)	2.0, 7.6, 12.6, 13.5, 14.3, 14.8, 15.4, 16.5, 18.7	2.0, 11.0, 23.0, 29.0, 35.0, 42.0	2.0, 12.0, 22.0, 28.0, 33.0, 38.0, 44.0
Net Radiation
Instrument	CNR1 net radiometer (Kipp & Zonen)	CNR1 net radiometer (Kipp & Zonen)	CNR1 net radiometer (Kipp & Zonen)
Height (m)	Initially 20 m, then changed to 22 m in August, 2013	42 m	42 m
Soil Temperature
Instrument	Pt 1000	Pt 100	Pt 1000
Depth (m)	0, 0.05, 0.10, 0.20, 0.30, 0.50	0, 0.05	0, 0.05, 0.10, 0.20, 0.30, 0.50

Environmental Measuring System—EMS.

.

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 (${}^{\circ }$C) and relative air humidity (RH; %) according to [31]. The daily sum of PAR (MJ m⁻² day⁻¹) was derived following the work in [32].

2.3. 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⁶), with the first 1.8 × 10⁶ simulations based on the prior distribution and the remaining 3 × 10⁶ 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ý Kr̆í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³ m⁻³ at the wet spruce forest site in Bílý Kr̆íz̆, 0.017 m³ m⁻³ for the dry spruce forest site in Rájec and 0.032 m³ m⁻³ for the beech forest site in Štítná, suggesting realistic SVWC estimates.

2.4. Data Processing and Analysis

2.4.1. 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₂, 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₂ 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₂ 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 (($\mathsf{\mu }$mol m⁻² s⁻¹)) values. The half-hourly GPP values were aggregated to obtain daily and monthly sums of GPP values.

2.4.2. 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).

2.4.3. 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.

2.4.4. 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²) in relation to estimated and predicted GPP_norm, and root mean square error (RMSE) values were derived and compared.

3. Results

3.1. 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ý Kr̆í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ý Kr̆í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ý Kr̆í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ý Kr̆íz̆.

3.2. 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 2. Results show that the maximum GPP_pot value of approximately 15.81 gC m⁻² day⁻¹ at the wet spruce forest site in Bílý Kr̆íz̆ was the highest across all the forest sites. This maximum GPP_pot value was obtained on DOY values from 197–203 (within July) characterized by lower mean T_air and high SVWC values (Figure 1). Additionally, maximum GPP_pot values of approximately 12.97 gC m⁻² day⁻¹ and 14.48 gC m⁻² day⁻¹ were obtained during June for the dry spruce (with DOY values from 170–177) and beech (with DOY values from 165–166 and 178) forest sites in Rájec and Štítná, respectively. Furthermore, the seasonal variations in the daily derived GPP_norm values over the growing season across all forest sites were compared for each of the studied year (Figure 3). Lower GPP_norm values were found in 2015 depicting GPP depression especially at the dry spruce forest site in Rájec during the main growing season of that year. The GPP_norm values were comparatively higher at the beech forest site in Štítná than at both spruce forest sites over the entire study period (Figure 3).

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 (Figure A2 and Figure 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.

3.2.1. 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. Statistical values from the stepwise multiple linear regression model (SMLR) and random forest (RF) analyses showing the accuracy indicators in predicting the Normalized GPP (GPP_norm) values across the wet spruce forest in Bílý Kr̆íz̆ (CZ-BK1), the dry spruce forest in Rájec (CZ-RAJ) and the beech forest in Štítná (CZ-Stn) for May–September of 2012–2016. The Pearson correlation, R², and root mean square error (RMSE) values have all been shown.

Variants	CZ-BK1		CZ-RAJ		CZ-Stn
	SMLR	RF	SMLR	RF	SMLR	RF
Pearson correlation	0.63	0.76	0.71	0.81	0.67	0.75
R2	0.40	0.57	0.49	0.65	0.45	0.55
RMSE	0.18	0.15	0.14	0.12	0.15	0.13

). 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² 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 (

Table A3. Results of the stepwise multi-linear regression model on the significant environmental variables that influenced the normalized gross primary productivity (GPP_norm) at the wet spruce forest site in Bílý Kr̆íz̆ (CZ-BK1). CNI: clearness index; PAR: photosynthetic available radiation; T_air: air temperature; T_s: soil temperature; SVWC: the soil volumetric water content; and VPD: vapour pressure deficit.

	Estimate	Standard Error	t Value	Pr (>|t|)
(Intercept)	$-6.397\times {10}^{-1}$	$1.397\times {10}^{-1}$	−4.581	$5.47\times {10}^{-6}$ ***
PAR	$1.659\times {10}^{-1}$	$1.541\times {10}^{-2}$	10.772	<$2\times {10}^{-16}$ ***
PAR2	$-8.813\times {10}^{-3}$	$1.025\times {10}^{-3}$	8.601	<$2\times {10}^{-16}$ ***
SVWC2	$1.095\times {10}^1$	1.938	−5.649	$2.34\times {10}^{-8}$ ***
PAR:Tair	$-5.234\times {10}^{-3}$	$9.454\times {10}^{-4}$	−5.536	$4.36\times {10}^{-8}$ ***
SVWC	$5.048\times {10}^1$	$9.514\times {10}^{-1}$	5.305	$1.51\times {10}^{-7}$ ***
Ts:PAR	$4.884\times {10}^{-3}$	$1.062\times {10}^{-3}$	4.598	$5.04\times {10}^{-6}$ ***
Ts2	$-2.002\times {10}^{-3}$	$4.817\times {10}^{-4}$	−4.155	$3.64\times {10}^{-5}$ ***
VPD:SVWC	$9.490\times {10}^{-2}$	$2.389\times {10}^{-2}$	3.972	$7.84\times {10}^{-5}$ ***
CNI2	$-3.623\times {10}^{-1}$	$1.025\times {10}^{-1}$	−3.533	0.000437 ***
CNI	$4.589\times {10}^{-1}$	$1.326\times {10}^{-1}$	3.461	0.000571 ***
VPD2	$-2.191\times {10}^{-3}$	$7.005\times {10}^{-4}$	−3.128	0.001830 ***
VPD	$-3.745\times {10}^{-2}$	$1.525\times {10}^{-2}$	−2.456	0.014300 *
Ts	$2.550\times {10}^{-2}$	$1.077\times {10}^{-2}$	2.367	0.018214 *
VPD:PAR	$3.094\times {10}^{-3}$	$1.367\times {10}^{-3}$	2.264	0.023906 *
CNI:Tair	$1.066\times {10}^{-2}$	$5.007\times {10}^{-3}$	2.128	0.033690 *
CNI:PAR	$-1.544\times {10}^{-2}$	$7.571\times {10}^{-3}$	−2.039	0.041772 *
VPD:Tair	$1.937\times {10}^{-3}$	$9.802\times {10}^{-4}$	1.976	0.048507 *

Signif. code: p < 0.001 ‘***’; p < 0.05 ‘ *’.

and

Table A4. Results of the stepwise multi-linear regression model on the significant environmental variables that influenced the normalized gross primary productivity (GPP_norm) at the dry spruce forest site in Rájec (CZ-RAJ). CNI: clearness index; PAR: photosynthetic available radiation; T_air: air temperature; T_s: soil temperature; SVWC: the soil volumetric water content; and VPD: vapour pressure deficit.

	Estimate	Standard Error	t Value	Pr (> |t|)
(Intercept)	1.2416627	0.2834895	4.380	$1.39\times {10}^{-5}$ ***
PAR	0.1957506	0.0157165	12.455	<$2\times {10}^{-16}$ ***
PAR2	−0.0055807	0.0006615	−8.436	<$2\times {10}^{-16}$ ***
Ts:SVWC	0.2931591	0.0448089	6.542	$1.25\times {10}^{-10}$ ***
Tair2	−0.0009479	0.0001752	−5.409	$8.97\times {10}^{-8}$ ***
CNI2	−0.9090371	0.1592056	−5.710	$1.74\times {10}^{-8}$ ***
Ts2	−0.0015452	0.0002921	−5.290	$1.68\times {10}^{-7}$ ***
SVWC:CNI	4.0751744	0.8162736	4.992	$7.71\times {10}^{-7}$ ***
Ts:Tair	0.0016139	0.0003835	4.208	$2.94\times {10}^{-5}$ ***
PAR:SVWC	−0.352378	0.0896162	−3.932	$9.35\times {10}^{-5}$ ***
SVWC	−8.4537169	2.2630082	−3.736	0.000204 ***
VPD	−0.0364140	0.0102154	−3.565	0.000392 ***
Ts	0.0404096	0.0149492	−2.703	0.007053 ***
VPD:SVWC	0.1377881	0.0518878	2.656	0.008118 ***
SVWC2	10.9036492	5.0661812	2.152	0.031754 *

Signif. code: p < 0.001 ‘***’; p < 0.05 ‘ *’.

). 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ý Kr̆í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. Results of the stepwise multi-linear regression model on the significant environmental variables that influenced the normalized gross primary productivity (GPP_norm) at the beech forest site in Štítná (CZ-Stn). CNI: clearness index; PAR: photosynthetic available radiation; T_air: air temperature; T_s: soil temperature; SVWC: the soil volumetric water content; and VPD: vapour pressure deficit.

	Estimate	Standard Error	t Value	Pr (> |t|)
(Intercept)	1.6324051	0.2913214	5.603	$3.00\times {10}^{-8}$ ***
PAR	0.1035450	0.0089856	11.523	<$2\times {10}^{-16}$ ***
Ts:SVWC	0.4954805	0.0627527	7.896	$1.08\times {10}^{-14}$ ***
Ts	−0.1794898	0.0233970	−7.671	$5.56\times {10}^{-14}$ ***
CNI	0.8213389	0.13378886	5.957	$4.04\times {10}^{-9}$ ***
Tair	0.0988082	0.0178532	5.534	$4.38\times {10}^{-8}$ ***
SVWC	−6.3883130	1.3261064	−4.817	$1.78\times {10}^{-6}$ ***
Tair:SVWC	−0.2274474	0.0467787	−4.862	$1.43\times {10}^{-6}$ ***
Tair:PAR	−0.0024049	0.0005103	−4.713	$2.93\times {10}^{-6}$ ***
CNI:PAR	−0.0558883	0.0148897	−3.753	0.000189 ***
VPD2	−0.0008560	0.0002309	−3.708	0.000225 ***
VPD:Ts	0.0021074	0.0006005	3.510	0.000477 ***
VPD	−0.0561466	0.0176605	−3.179	0.001540 ***
SVWC2	4.9409576	1.5641808	3.159	0.001651 ***
VPD:SVWC	0.1048877	0.0339888	3.086	0.002107 ***
CNI2	−0.5376705	0.2101539	−2.558	0.010718 ***

Signif. code: p < 0.001 ‘***’.

).

Moreover, similar quadratic terms such as PAR², T_s², and CNI² 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 (Table A3–Table A5). Furthermore, at the beech forest site in Štítná, VPD² and CNI² proved to be inversely related to GPP_norm and statistically significant (Table A5). However, the positive coefficient value for SVWC² across all the forest sites showed an increase in GPP_norm values with increasing SVWC values (Table A3, Table A4 and Table 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ý Kr̆í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.

3.2.2. GPP_norm Prediction through Random Forest Analyses (RF)

The validation statistics of the RF modeling are shown in Table 2. Note that higher R² (>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).

3.2.3. 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ý Kr̆í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ý Kr̆íz̆ (with the highest variable importance score of 0.18 for VPD).

A decision tree representation showing regions with high and low GPP_norm daily values from the RF analyses were plotted for each forest ecosystem (Figure 7).

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⁻² day⁻¹), while high GPP_norm values were attainable at the beech forest in Štítná above a PAR threshold of approximately 4 MJ m⁻² day⁻¹. 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ý Kr̆í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⁻² day⁻¹) values at the wet spruce forest in Bílý Kr̆íz̆ (6 cases) and the beech forest site in Štítná (7 cases). Unlike in the wet spruce and beech forest sites, lower GPP_norm values were obtained at the dry spruce forest site in Rájec on periods with low SVWC (<0.16 m³ m⁻³; 4 cases) and CNI (<0.2; 13 cases).

Moreover, apart from PAR, other environmental variables were also observed to be important for group separation. SVWC (≤0.16 m³ m⁻³; 250 cases) and T_air (≤19 ${}^{\circ }$C; 28 cases) values had secondary effects on GPP_norm values at the wet spruce forest Bílý Kr̆íz̆. At the dry spruce forest site in Rájec, other limiting factors for GPP_norm values were also found to be T_s (≤20 ${}^{\circ }$C; 255 cases) and SVWC (≤0.16 m³ m⁻³; 24 cases). There were periods (193 cases) when higher GPP_norm values were realized with T_s< 20 ${}^{\circ }$C and PAR threshold value of approximately 6 MJ m⁻² day⁻¹ at the dry spruce forest site. As was observed in both Figure 6 and Figure 7, SVWC (≤0.26 m³ m⁻³; 255 cases) and VPD (≤9 hPa; 51 cases) values were found to have secondary effects on GPP_norm values at the beech forest site. However, at similar PAR threshold values of approximately 5 MJ m⁻² day⁻¹ and different SVWC threshold values, higher GPP_norm values were realized at both the wet spruce (SVWC values of approximately 0. 16 m³ m⁻³; 197 cases) and beech forest sites (SVWC values of approximately 0.26 m³ m⁻³; 192 cases).

4. 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ý Kr̆í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ý Kr̆íz̆, T_air split the subspace with lower values of GPP_norm at high T_air threshold values above 19 ${}^{\circ }$C during days with lower SVWC values (<0.16 m³ m⁻³). This suggests that high T_air values limit the forest productivity at the mountainous wet spruce forest site in Bílý Kr̆í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ý Kr̆í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 (Figure A1 and Figure 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ý Kr̆í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 (Figure A3 and Figure A4) [19,76].

5. 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ý Kr̆í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.