**2. Materials and Methods**

#### *2.1. Station Description*

The study used multiyear (2014–2016) measurements from Bílý Krí˘ z and Rájec ecosys- ˘ tem stations that are part of the CzeCOS (Czech Carbon Observation System; http:// www.czecos.cz/ (accessed on 6 March 2021)) and FLUXNET (Flux Tower Network; https: //fluxnet.fluxdata.org/ (accessed on 6 March 2021)) station network. The FLUXNET station IDs for the wet and dry spruce forest are CZ-BK1 and CZ-RAJ, respectively. CZ-BK1 is also a candidate for ICOS (Integrated Carbon Observation System; https://www.icos-cp.eu/ (accessed on 6 March 2021)) network. The main vegetation cover in both stations is the even-aged stand of Norway spruce. Their main characteristics are presented in Table 1. CZ-BK1 is situated in the Moravian-Silesian Beskids Mountains in the Czech Republic and has a cold and humid climate. CZ-RAJ is situated in the Drahany Highland at a lower altitude and has a moderately warmer and drier climate than CZ-BK1. 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 2014 to 2016 (Table 1).

**Table 1.** Characteristics of the study sites.


#### *2.2. Eddy Covariance and Ancillary Measurements*

At each station, the EC system consisted of an infrared gas analyser (LI-COR, Lincoln, NE, U.S.A.) and an ultrasonic anemometer (Gill, Hampshire, U.K.) measuring at a 20 Hz frequency (models are specified in Table 2). Each system was mounted on a tower several meters above the forest canopy. The EC measurements were complemented by an extensive set of sensors to collect the required auxiliary meteorological data. Measurements included the air temperature (Tair) and humidity profile with EMS33 temperature and humidity sensors (EMS, Brno, Czech Republic) and precipitation using Precipitation Gauge 386C (Met One Instruments, Grants Pass, OR, U.S.A.). Measurements for the incoming photosynthetic active radiation was made using LI-190R Quantum Sensor (LI-COR, Lincoln, NE, U.S.A.) at Bílý Krí˘ z and EMS12 sensor (EMS, Brno, Czech Republic) at Rájec. Additionally, profiles ˘ of soil moisture were measured by using the CS616 (Campbell Scientific, Inc., Logan, UT, U.S.A.) sensors. An overall description of the instrumentation at the study sites is given in Table 2.

The *VPD* for each site was computed from Tair and relative humidity (*RH*) as follows [39]:

$$VPD = \frac{100 - RH}{100}(SVP) \tag{1}$$

where *SVP* (Pa) is the saturated vapour pressure given as follows:

$$SVP = 610.7 \times 107.5^{\frac{lur}{277.5 + lur}} \tag{2}$$

**Table 2.** Description of the eddy covariance systems at the investigated sites.


#### *2.3. Data Processing and Analysis*

2.3.1. Turbulent Flux Measurements

At both study sites, the EC data from the main growing season (May–September) of 2014–2016 were used for the analysis. Gas analyser measurements of the density of the water vapour in the air and ultrasonic anemometer measurements of three-dimensional wind components and sonic temperature were made [35,36]. The eddy covariance processing was performed, using an open-source software, EddyPro 7.0.6 (Li-COR, Lincoln, NE, U.S.A.). The most recent methods for flux corrections, conversions, and a thorough quality control scheme as proposed by [40] were also applied. This process involves despiking and raw data statistical screening, basic quality checking (QC) of turbulent fluxes (flux stationarity and integral turbulence characteristics tests), coordinate rotation using the planar fit method [41], spectral correction [42–44], detecting and compensating for time lags of signals from the ultrasonic anemometer and the gas analyser, footprint estimations and calculating half-hourly fluxes. In EC data post-processing, determining periods with low turbulent mixing is a critical step [35]. Flux measurements over periods with insufficiently developed turbulence, i.e., low friction velocity (u∗) need to be detected and filtered out. This filtering procedure assured the exclusion of CO2 measurements not representative of the biotic flux, i.e., net ecosystem exchange (NEE) [31,45,46].

The R software [47] package 'REddyProc' [48] was used to gap-fill the EC data at both sites, using marginal distribution sampling (MDS) [46]. This way, all missing values were replaced by the average value estimated under similar meteorological conditions within a specific time window. If no similar conditions were available within the starting time window of about 7 days, the window was increased. NEE was partitioned into GPP and

ecosystem respiration (Reco). The flux partitioning approach by [49] using daytime data was applied to estimate half-hourly GPP (μmol m−<sup>2</sup> s−1) values. The half-hourly GPP values were aggregated to obtain daily and monthly GPP sums.

#### 2.3.2. Soil Water Content Simulations

Since the soil volumetric water content (SVWC) was not measured throughout the entire period of our analyses (the measurements started at the beginning of 2016), we used a simulated daily SVWC instead. We applied the soil water balance model R-4ET, an R package for Empirical Estimate of Ecosystem EvapoTranspiration as used in [50]. The calibration of the soil water balance model was carried out using Bayesian statistics implemented in the R package BayesianTools [51] with a Differential-Evolution Markov chain Monte Carlo sampler. In Bayesian calibration, the model input parameters are iteratively updated to provide a probability distribution of the calibrated parameters that represents the uncertainty in the measured data and in the model structure. The simulations were repeated 4.8 × 106 times, where the first 1.8 × 106 runs were based on the prior distribution, while the remaining 3 × 106 were constrained by the posterior distribution resulting from the first set of runs. A number of iterations treated as burn-in were set to 0, and the thinning parameter determining the interval in which values are recorded was set to 1. The high number of iterations was necessary to attain a good input parameter convergence with a narrow distribution. For inspecting the convergence, we used Gelman diagnostics with a criterion for the potential scale reduction factor being less than 1.2 for all parameters [52,53]. The final selection of parameters from their probability distributions was conducted by using a maximum a posteriori probability estimate, i.e., the value representing the mode of the posterior distribution.

The main model input variables included meteorological data and the leaf area index (see [50], for more details). The soil was stratified into 12 layers, which increased in thickness with increasing depth down to 3 m. Note that the soil layers were selected in the way that the layers used for optimization matched the depths of all sensors with an extra ±2.5 cm, considering the volume measured by the sensors. At both sites, the available SVWC measurements were available throughout the main part of the root zone. At CZ-BK1, the sensors CS616 (Campbell Scientific, Inc., Logan, UT, U.S.A.) were placed at 0.05, 0.1, 0.22, 0.34, and 0.42 m soil depths and measured since 2016. At CZ-RAJ, the same type of sensors were also placed at 0.05, 0.1, 0.2, 0.5, and 0.8 m depths and measured since 2016. The soil parameters, including SVWC at saturation, field capacity, wilting point, and saturated hydraulic conductivity, were optimized at all 12 depths [50]. Additional single parameters relevant for the SVWC simulations that were optimized included the following: rooting depth with the Beta parameter describing the root profile shape [54], surface resistance and the degree of isohydricity [55], the water interception capacity of leaf and bark area, and curve number representing a runoff parameter [50]. The uniform distribution of priors was used, where the lower and upper limits were set to be within ±50% of the values based on field measurements of the wilting point, field capacity, and saturated water content and ±100% for the remaining parameters that were estimated from the literature or from previous anecdotal analyses. To provide the model with sufficient spin up time to stabilize and ensure reliable and robust parametrization, the simulation was initiated at the beginning 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 at both sites, according to a computed standardized precipitation-evapotranspiration index over the region). The overall simulation was done for the period 2010–2019 where the observed data for Bayesian calibration spanned the period 2016–2019.

The simulated SVWC averaged over all depths yielded a root mean square error of 0.037 m<sup>3</sup> m−<sup>3</sup> and 0.017 m3 m−<sup>3</sup> for CZ-BK1 and CZ-RAJ, respectively, suggesting realistic SVWC estimates.
