Using Sap Flow Data to Parameterize the Feddes Water Stress Model for Norway Spruce

Abstract

Tree water use is a key variable in forest eco-hydrological studies and is often monitored by sap flow measurements. Upscaling these point measurements to the stand or catchment level, however, is still challenging. Due to the spatio-temporal heterogeneity of stand structure and soil water supply, extensive measuring campaigns are needed to determine stand water use from sap flow measurements alone. Therefore, many researchers apply water balance models to estimate stand transpiration. To account for the effects of limited soil water supply on stand transpiration, models commonly refer to plant water stress functions, which have rarely been parameterized for forest trees. The aim of this study was to parameterize the Feddes water stress model for Norway spruce (Picea abies [L.] Karst.). After successful calibration and validation of the soil hydrological model HYDRUS-1D, we combined root-zone water potential simulations with a new plant water stress factor derived from sap flow measurements at two plots of contrasting soil moisture regimes. By calibrating HYDRUS-1D against our sap flow data, we determined the critical limits of soil water supply. Drought stress reduced the transpiration activity of mature Norway spruce at root-zone pressure heads <−4100 cm, while aeration stress was not observed. Using the recalibrated Feddes parameters in HYDRUS-1D also improved our water balance simulations. We conclude that the consideration of sap flow information in soil hydrological modeling is a promising way towards more realistic water balance simulations in forest ecosystems.

1. Introduction

Tree water use is a key variable in forest eco-hydrological studies and closely linked to meteorological conditions and soil water availability [1,2,3,4]. To monitor the water use of individual trees, a number of methods are available [5,6]. Sap flow measurements are most popular, not only because they are cost-efficient and easy to apply. They also yield reasonable information on temporal transpiration patterns and deliver valuable insights into plant-physiological responses to varying hydro-climatic conditions (e.g., soil moisture limitations). Respectively, [7] identified sap flow as a “key trait in the understanding of plant hydraulic functioning”.

However, to make stem-based sap flow measurements valuable for watershed and land use managers, they need to be scaled to a ground area basis [2]. Common scalars in this regard are sapwood area, basal area, diameter or circumference at breast height and leaf area [2,8,9]. However, the scaling of sap flow measurements from the tree to the forest stand includes many uncertainties. To begin with, different measuring approaches yield differing estimates of sap flow rates per tree [9,10,11,12]. This is partly related to uncertainties in signal transformation [13,14,15], but also to radial and circumferential variations in sap flow density, which complicates the scaling from the point measurement to water use by tree [16,17,18,19,20,21]. Furthermore, sap flow activity and trans-sectional sapwood distribution vary with tree age, tree dimensions, social position of the tree and also with soil water supply [4,22,23,24,25]. Thus, extensive measurement campaigns are needed to get a solid estimate of stand water use from sap flow measurements alone.

Alternatively, water balance models can be used to estimate stand transpiration. Such models basically calculate transpiration as a function of atmospheric boundary conditions (potential evapotranspiration), soil water supply and vegetation characteristics from which leaf area index (LAI) and stomatal conductance are the most important [2,26]. Under non-limiting conditions, transpiration generally follows the dynamic of atmospheric demand. However, the response of plants to variabilities in soil water supply is highly species specific. To account for these species-specific feedbacks in the soil-vegetation-atmosphere system, water balance models commonly include plant water stress functions (e.g., [27,28]), where critical limits of soil water supply can be adapted according to the vegetation characteristics. For many agricultural crops, such thresholds have already been defined and included into models that help to optimize irrigation scheduling and crop yields [3]. However, little is known about the species-specific water-stress-response of forest trees. Since forests play an important role in regional and trans-regional water cycles, there is evidence that an improved knowledge of tree species specific water stress thresholds would greatly enhance the simulations of water fluxes from the forest stand to the watershed and landscape scale [2,9,29,30,31].

The aim of this study was: (1) to parameterize the Feddes water stress model for Norway spruce, which is a wide-spread and economically important species in Europe; (2) to investigate if the species-specific calibration of the Feddes parameters can improve water balance simulations in the soil hydrological model HYDRUS-1D [32]. In a first step, we used HYDRUS-1D to estimate the root zone water potential for two research plots with contrasting soil water regimes. We derived the Feddes parameters for Norway spruce from correlations of the simulated water potentials and a new sap flow-based plant water stress factor.

2. Materials and Methods

2.1. Study Area

The study was conducted in the 27 ha Wüstebach experimental test site [33], which is located in the Rur River catchment close to the German-Belgian border (Figure 1). The test site belongs to the TERENO Eifel/Lower Rhine Valley Observatory [34,35]. Altitudes range from 595 m a.s.l. in the north to 628 m a.s.l. in the south. While hillslopes are dominated by shallow Cambisols and Planosols, Gleysols and Histosols have developed in the groundwater-influenced riparian zone along the Wüstebach stream. The soils mainly show a silty clay loam texture with a medium to high coarse material fraction and a litter layer on top [36]. The soils in the valley bottom have slightly higher bulk densities and lower (macro-) porosities and skeleton contents than those on the hillslopes [36,37,38]. The climate is characterized by an annual mean temperature of 7 °C [39] and a mean annual precipitation of 1220 mm [33]. The annual potential evapotranspiration of the years 2010–2013 was 630 mm [40]. The catchment is forested with Norway spruce that were planted in the late 1940s [41]. The trees have now reached a height of about 28 m and the tree density is about 320 trees ha⁻¹.

2.2. Data and Experimental Design

We used sap flow data from two research plots (Figure 1) from which one (Riparian) is located in the wetter valley bottom (~600 m a.s.l.), while the other is situated on the drier hillslope (Slope see Figure 1) at ~610 m a.s.l. With a gradient of 2°, the Riparian plot is slightly inclined towards N, while the Slope plot shows a gradient of 8° in NNW direction. Prevailing soil types are Gleysol (Riparian) and Cambisol/Planosol (Slope).

At each plot, sap flow of four trees was monitored over 2 years (2009–2010) using improved Granier-type sap flow sensors (SF-L 20/33, Ecomatik, Dachau, Germany). The classic Granier system consists of two sensor probes inserted radially into the sapwood, one above the other. The upper probe is equipped with a heating element and a thermocouple, thus recording the heat dissipation due to sap flow. The lower probe measures the ambient reference temperature of the wood. Our improved sensor version includes another pair of thermocouples that is placed horizontally to the upper heated probe to account for natural inner-wood temperature variations. The mean of the inner-wood temperature variations recorded by these additional reference probes is subtracted from the values recorded by the classic Granier system before applying the Granier formula, where sap flux density is derived from the temperature gradient between the heated and the lower reference probe using the empirical equation [42,43]:

$$\mathrm{SFD}=119\times {\left(\frac{\Delta T_{max}-\Delta T}{\Delta T}\right)}^{1.231}$$

(1)

in which SFD is the sap flux density (g·m⁻²·s⁻¹), ∆T is the actual temperature gradient between the two probes and ∆T_max the maximum temperature gradient measured between the probes in a given time period. Following [15], who found a daily ∆T_max determination suitable for the study area, we determined ∆T_max on the basis of a 24 h moving window (12 h before and 12 h after the actual point of measurement). This dynamic ∆T_max determination procedure appeared the best opportunity to handle potential drifts in ∆T_max.

The sap flow sensors were installed in the outermost 3.3 cm of the sapwood on the north side of the sample trees at ~1.5 m above ground. To ensure that the sensor probes were completely surrounded by sapwood (i.e., no heartwood), we determined each tree’s sapwood depth by drill hole analyses before installing the sap flow sensors. We insulated our probes with reflective polystyrene and plastic boxes. The temperature gradients measured and recorded at a datalogger (type CR1000, Campbell Scientific Ltd., Logan, UT, USA) in 30-min time steps.

Data gaps were filled using a simple three step procedure: (1) linear regression among all sap flow series within one plot; (2) choice of the data series with the highest coefficient of determination (R²) with the data series to be gap-filled; (3) interpolation of missing values from the respective regression line. Mean daily sap flow density was calculated from the gap-filled hourly data

Soil moisture was monitored by the wireless sensor network SoilNet, which provides catchment-wide information on soil water dynamics in 5, 20 and 50 cm depths at 150 measuring locations (see Figure 1). The installed measuring devices are the capacitance soil water content sensors ECH2O EC-5 and ECH2O 5TE (Decagon Devices, Pullman, WA, USA) [44]. Sensor calibration was performed according to [37]. For this study, we used the mean daily soil moisture at 20 and 50 cm depth of two SoilNet stations located close to our sap flow plots (Figure 1). Although there are other SoilNet stations closer to the Slope plot than the chosen station SN 28, we used this station, because it is located at the same elevation as the sap flow plot. The SoilNet data from 5 cm depth was not considered in this study, because the respective sensors were located in the litter layer and not in the mineral soil we focus on (Figure 2).

Grass reference evapotranspiration ET₀ was calculated using the FAO-Penman-Monteith method [45] following [40], who first calculated hourly ET₀ values and then aggregated these to daily sums. The required input data were obtained from the TERENO weather station Schöneseiffen (3.4 km east of the Wüstebach site). For gap-filling, linear regression predictors from the same variables of the TERENO station Selhausen (40 km north) were used. To calculate potential evapotranspiration for specific crops (not grass), ET₀ is usually multiplied with the crop coefficient K_c. However, according to [45] a K_c value of 1 can be used for coniferous forests. Therefore, we assume ET₀ to represent potential forest transpiration in this study. Daily precipitation was obtained from the Kalterherberg weather station of the German Meteorological Service (DWD), which is located 9.6 km east from of the test site and representative for our study site [46]. The precipitation data were corrected for systematic measurement errors according to [47].

The mean leaf area index was determined by averaging respective measurements (SunScan SS1, Delta-T Devices Ltd., Cambridge, UK) at 60 selected SoilNet locations within the Wüstebach site.

2.3. Estimating Root-Zone Water Potential

The water potential of the root-zone is a measure for the energy needed to extract water from the soil [48]. Therefore, we consider the root-zone water potential to be a solid measure to characterize the water availability for plants. However, since water potential is not measured at the SoilNet stations, we used the 1-D Richards equation as implemented in the HYDRUS-1D software [32] to inversely determine water potential dynamics from the soil water content (SWC) observations. This is a common approach in soil hydrology, which also allows for the integration of the point scale SWC measurements into vertical soil moisture profiles. In HYDRUS-1D the soil hydraulic properties are parameterized with the Mualem-van Genuchten model [49]:

$$\theta \left(h\right)=\{\begin{array}{ll}{\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left[1+{\left|\alpha h\right|}^n\right]}^m}& h<0\\ \\ \\ {\theta }_s& h\ge 0\\ \end{array}$$

(2)

$$K\left(h\right)=K_s{S}_e^{0.5}{\left[1-{\left(1-S_e^{0.5/m}\right)}^m\right]}^2$$

(3)

where

$$m=1-1/n\hspace{1em}\hspace{1em}n>1$$

(4)

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(5)

where h is the pressure head (cm), θ_r is the residual water content (cm³ cm⁻³), θ_s the saturated water content (cm³ cm⁻³), α and n are empirical parameters, from which α is related to the air-entry pressure value and n is related to the pore size distribution, m is the retention curve shape parameter, which is dependent on n, K_s is the saturated hydraulic conductivity (cm h⁻¹), and S_e is the relative soil saturation.

Both research plots were modelled separately and in daily time steps. As top boundary conditions, we used the corrected precipitation data and the gap-filled potential evapotranspiration data (see above). The latter was divided into potential evaporation (E₀) and potential transpiration (T₀) on the basis of the leaf area index (LAI) [45]. In this study, mean LAI on site was 4.8 (with a standard deviation of 0.7) and assumed to be constant in time. Therefore, the differentiation between E₀ and T₀ was of minor importance. However, applying our approach to deciduous trees with seasonal LAI dynamics would require the consideration of the T₀ limitation by LAI. For the hillslope plot (Slope) we set the lower boundary condition to free drainage, following [50], who found this option to provide the best HYDRUS-1D simulation results for the Wüstebach catchment.

To account for the groundwater fluctuations in the riparian zone, the lower boundary condition for the wetter plot (Riparian) was set to deep drainage [51]. For the deep drainage option, the discharge rate at the bottom of the soil profile is defined as a function of the position of the groundwater table and the empirical parameters A_qh and B_qh [52]:

$$q\left(h\right)=A_{qh} exp\left(B_{qh}\times \left|h-GW{L}_{ref}\right|\right)$$

(6)

in which q(h) (cm day⁻¹) is the discharge rate, h (cm) is the pressure head at the bottom of the soil profile, A_qh (cm day⁻¹) and B_qh (cm⁻¹) are empirical parameters and GWL_ref (cm) is the reference groundwater depth. We manually optimized A_qh and B_qh with regard to the best fit (root-mean-square error) between measured and simulated soil moisture dynamics.

Depending on the local site conditions [36,38] we discretized the HYDRUS-1D soil profiles into two (Riparian) and three (Slope) materials (one/two soil horizons plus an organic litter layer on top) (Figure 2). The hydraulic properties of the litter layer were parameterized according to [53]. For the other soil materials, we inversely estimated the Mualem-van-Genuchten parameters from our measured SWC in 20 and 50 cm depth. We therefore used the inverse solution option as implemented in HYDRUS-1D, where the objective function Φ is minimized applying the Levenberg-Marquardt nonlinear least-squares fitting approach (for detailed information on the optimization procedure, see [32]).

Root water uptake in HYDRUS-1D depends on the root distribution over the soil depth and on soil water availability [32]. For our Norway spruce plots, we assumed a rooting depth of 45 cm, which is in line with [54,55]. The root zone was assumed to start below the litter layer, which had a thickness of 6 (Slope) and 9 cm (Riparian), respectively (Figure 2). The root density in spruce forests is typically highest in the upper soil layers [56,57,58]. Therefore, we assumed the potential root water uptake factor to be 1 in the uppermost 15 cm of the mineral soil (one third of the root zone). For the deeper layers, we implemented a linear decrease of the potential root water uptake factor from 1 in 15 cm depth to 0 in 45 cm depth.

To account for reduced water uptake rates in case of water stress (aeration stress and drought stress), we used the Feddes plant water stress function [28] as implemented in HYDRUS-1D. In the Feddes function, a dimensionless water stress factor represents the water stress level of the plant as a function of the root zone water potential (Figure 3). This water stress factor corresponds to the ratio of actual (T_act) and potential (maximum possible) transpiration (T₀) and thus ranges from 0 (maximum stress/no root water uptake) to 1 (no stress/optimum root water uptake). Optimum root water uptake (T₀ = T_act) is given between h₂ (field capacity FC) and h₃ (onset of drought stress). The decreasing plant water stress between h₂ and h₁ considers that plants typically experience increasing aeration stress towards saturated conditions; h₄ corresponds to the permanent wilting point (WP), h₁ to the point of maximum aeration stress (anoxic conditions).

Presuming that h₃ decreases with decreasing T₀, HYDRUS-1D allows for making h₃ a function of T₀. The range of T₀, within which h₃ becomes a function of T₀ may be defined by the user, but is usually 1–5 mm day⁻¹ (default values). The h₃ value for T₀ < 1 mm day⁻¹ is called h_3,low, that for T₀ > 5 mm day⁻¹ h_3,high. If T₀ is between 1 and 5 mm day⁻¹ and the soil dries out to less than h_3,high, the critical root-zone water potential, at which potential root water uptake is reduced to the actual root water uptake, is linearly related to T₀ [32].

As start values for the model inversion, we used the Feddes parameters (in cm pressure head) as proposed by [59] for a 90-years-old spruce site in the Czech Republic, where h₁ = h₂ = 0 cm, h_3,high = −600 cm, h_3,low = −1200 cm, and h₄ = −15,000 cm. Later, we manually optimized the Feddes parameters on the basis of our sap flow data (see below).

2.4. Model Validation

We validated our model setup on the basis of the soil moisture measurements in 20 and 50 cm depth using a split-sample approach. To consider the observed drop in soil moisture between the years 2010 and 2011 in the optimization procedure, we chose the years 2009–2011 as the calibration period. For validation, we used another three year period (2012–2014). The goodness-of-fit between observed and simulated data was assessed on the basis of the root mean squared error (RMSE). Our sap flow data was restricted to the years 2009 and 2010 and was therefore not considered for model validation.

2.5. Parameterizing the Feddes Model Using Sap Flow Data

In the classical Feddes approach, plant water stress is derived from the ratio of T_act and T₀ (see above). The assumption is that this ratio is mainly dependent on soil water supply. However, [60] showed that Norway spruce reacts to high midday vapor pressure deficits (VPD) by closing stomata. Thus, the trees physiologically reduce their own transpiration activity on many days of the growing season. In order not to mix-up this atmospheric stress reaction with actual soil water stress, we applied a new approach to calculate plant water stress from our sap flow data. This approach was based on the assumption that the micro-climate is identical at both plots; thus, in case of atmospheric stress, all trees (Riparian and Slope plot) show the same physiological reaction.

For each tree we determined the 6 highest SFD values per day. We averaged these maxima by day and plot and normalized the resulting data series (SFD_n,rip for Riparian and SFD_n,sl for Slope).

For the Riparian plot, we assumed optimum soil water supply, because due to the constant influence of groundwater the soil at this plot showed relatively high soil water contents throughout the study period and critical levels of soil water supply as reported from other studies [23,24,31,61,62] were never reached. We calculated a drought stress factor SFD_drought (1: optimum water supply; 0: maximum water stress) for the Slope plot by dividing SFD_n,sl by SFD_n,rip. To ensure that potential aeration stress at the Riparian plot would not affect our drought stress function, we only used data of dry days (3rd day without precipitation). We plotted SFD_drought against the simulated Slope water potential and analyzed the correlation between SFD_drought and root zone h for different levels of soil water shortage. The data which showed the highest significant correlation (R² and Student’s t-test) between SFD_drought and root zone h was used to determine the drought-stress trendline. The intersection of this trendline with the x-axis corresponds to the simulated Slope root zone water potential at maximum drought stress (SFD_drought = 0) and was therefore interpreted as the permanent wilting point h₄.

During model calibration we found that the parameter h_3,high was independent from mean daily sap flow dynamics and hence could not be optimized using sap flow measurements. However, it is well-known that at high levels of atmospheric demand, the transpiration rate of spruce forests decreases with increasing vapor pressure deficit [23]. This was also confirmed by our own sap flow measurements (data not shown). While coupled soil–plant–atmosphere continuum models are able to account for atmospheric stress, this process is not implemented in HYDRUS-1D. Nevertheless, since the Feddes approach is differentiated into h_3,high and h_3,low, atmospheric stress can be considered by adapting the T₀ range at which h₃ is a function of T₀. The study of [60] suggests that spruce trees reduce transpiration when midday vapor pressure deficit exceeds 1250 Pa independent of soil wetness. Accordingly, we analyzed days with maximum VPD of approx. 1250 Pa and found that a T₀ value of 2.7 mm day⁻¹ can be used as lower boundary for the h₃-dependent T₀ range. Furthermore, we used T₀ at the day of maximum observed daily SFD (4.8 mm day⁻¹) as the upper boundary for the h₃-dependent T₀ range.

To detect potential aeration stress at the Riparian plot, we selected days, for which no drought stress was observed on the Slope plot. For these days, SFD_n,sl was used as a reference for optimum water supply. We calculated an aeration stress factor SFD_aeration (1: optimum water supply; 0: maximum water stress) by dividing SFD_n,rip by SFD_n,sl. To determine the Feddes parameters h₁ and h₂, we plotted SFD_aeration against the simulated Riparian root-zone water potential and analyzed the observed data trend as described for SFD_drought.

2.6. Relative Extractable Soil Water (REW)

It is well known that absolute soil water contents do not fully reflect the actual soil water status [29,63]. Therefore, many researchers refer to the relative extractable soil water (REW) as a measure of soil water status (e.g., [4,23,63,64,65,66]). To compare our Feddes parameters with critical REW values from the literature, we calculated the mean root-zone water contents corresponding to the reported critical REW according to the relations described by [64,66]:

$$\mathrm{REW}=\frac{\theta -{\theta }_{PWP}}{{\theta }_{FC}-{\theta }_{PWP}}$$

(7)

in which θ is the actual water content, θ_FC is the water content at field capacity, and θ_PWP is the water content at permanent wilting point. To be consistent with the literature [64,66], we chose θ at −300 cm pressure head (−0.03 MPa) for θ_FC and θ at −16,300 cm pressure head (−1.6 MPa) for θ_PWP.

2.7. Uncertainty Assessment

To assess the uncertainty of the simulated root zone water potential related to the assumed rooting depth and the shape of the potential root water uptake function, we performed a sensitivity analysis.

First, we ran the model for each plot assuming constant rooting depth (45 cm), but differing root distributions (factor 1 in the uppermost 15/30/45 cm of the mineral soil and linear decrease from 1 to 0 towards the bottom of the root zone). Then, we varied the rooting depth (30/40/50/60 cm) and simulated the respective root zone water potentials assuming a root distribution factor of 1 for the uppermost 15 cm of the mineral soil linearly decreasing from 1 to 0 between 15 cm depth and the final rooting depth.

Test statistics (Mann-Whitney U-test) on the similarity of the simulated root zone water potentials were conducted on the log scale.

To evaluate the small-scale variability of soil water supply around the sap flow stations and to ensure that the selected SoilNet measuring locations are representative for our sap flow plots, we simulated the root zone water potentials of 4 (Riparian) and 7 (Slope) surrounding SoilNet locations (Figure 1). These simulations were again performed on the basis of calibrated and validated Mualem-van Genuchten parameters.

From the average root zone water potentials of the surrounding SoilNet locations, we determined the Feddes parameters as described before and compared the results with the water stress thresholds as resulting from the initial research setup.

3. Results

3.1. Sap Flow Data

The observed SFD dynamic of both research plots corresponds well to the dynamic of the calculated ET₀ (Figure 4). R² between mean daily SFD and ET₀ was 0.78 for the Slope plot and 0.84 for Riparian plot. The differing absolute SFD among trees and plots result from variations in tree diameter and crown size (

Table 1. Variation of diameter at breast height (DBH), sap wood depth (SWD), projected crown area (CA) and mean daily sap flow density between May and September (SFD) by plot and sample tree.

Sample Tree	Riparian				Slope
	DBH	SWD	CA	SFD	DBH	SWD	CA	SFD
	cm	cm	m²	g m−2 s−1	cm	cm	m²	g m−2 s−1
tree 1	57	6	49.2	4.48	44.7	4.7	56.3	4.52
tree 2	40.3	4.2	23.3	2.25	40.9	4.3	60.4	5.09
tree 3	46.3	4.9	22.2	1.56	44.1	4.6	33.8	4.69
tree 4	40.4	4.2	- *	3.13	48.2	5.1	- *	5.24
mean	46	4.8	31.6	2.9	44.5	4.7	50.2	4.89

* These trees were cut before our survey on CA determination; therefore we estimated CA from the allometric relationships of [67], where CA = π × (0.6122 + 0.0536 × DBH)².

), which can be explained by small-scale variations in planting density and uneven forest thinning. This becomes particularly apparent for tree 1 of the Riparian plot (Figure 4b).

At the Riparian plot, we observed some data inconsistencies in the winter that were probably caused by frost events and rodent activity. Therefore, the sap flow sensors were checked and readjusted at the beginning of the following vegetation period (May 2010). Since we did not use the noisy data until sensor readjustment, the winter gap at the Riparian plot is slightly longer than at the Slope plot (Figure 4).

3.2. Simulated Soil Water Dynamics

The Mualem-van Genuchten parameters as resulting from the parameter estimation by HYDRUS-1D show lower residual, but higher saturated water contents for the Riparian plot than for the Slope plot. At the Slope plot, α and K_s decrease with depth (

Table 2. Mualem-van Genuchten parameters as obtained from the HYDRUS-1D inverse modeling procedure by plot and soil layer. θ_r: residual water content, θ_s: saturated water content, α: empirical parameter related to the air-entry pressure value, n: empirical parameter related to the pore size distribution width, K_s saturated hydraulic conductivity (cm day⁻¹).

Soil Layer	θr	θs	α	n	Ks
	Riparian
B *	0	0.63	0.008	1.58	4
	Slope
A	0.1	0.52	0.008	1.22	165
B	0.12	0.5	0.001	1.29	78

* The A horizon was very thin and could therefore not be parameterized for this plot.

).

After inverse modelling, and manual calibration of A_qh and B_qh for the deep drainage option (−0.85 and −0.012, respectively), the root mean squared errors (RMSE) between observed and simulated data (Figure 5) were 3.6 vol % (Slope), and 4.6 vol % (Riparian). In the validation period, the RMSE of the Riparian plot slightly increased (5.3 vol %), while that of the Slope remained the same. The manual calibration of the Feddes parameters slightly improved the RMSE at the hillslope plot (3.4 vol % for both calibration and validation period), but had no effect on the RMSE of the Riparian plot simulations.

As expected, SWC in the groundwater-influenced plot (Riparian) were generally higher than those on the hillslope. The simulated mean annual SWC was 54% for the Riparian plot and 36% for the Slope). The simulated root zone water potentials (Figure 6) highlight the contrasting soil moisture regimes of the two research plots. The mean simulated annual root zone pressure heads were −88 cm (Riparian) and −3400 cm (Slope). All simulated daily root-zone h values ranged from +27 cm (fully saturated soil column; Riparian) to −12,577 cm (very dry conditions near wilting point; Slope).

3.3. Feddes Parameters for Norway Spruce (h₁, h₂, h_3,low, h₄)

For the Slope plot, we found the measured water stress to be dependent on the simulated root-zone water potential (Figure 7a). The strongest significant correlation between SFD_drought and simulated root-zone h was observed below −4100 cm (R² = 0.77, p < 0.001). Thus, we set h_3,low = −4100 cm. Extrapolating the trend-line between SFD_drought and simulated root-zone h towards the x-axis, where SFD_drought = 0 (maximum water stress), resulted in an h₄ value of −15,000 cm which is often used as—a very general—permanent wilting point [68]. Our data did not indicate the occurrence of aeration stress at our research site (Figure 7b). Parameters h₁ and h₂ were therefore set to 0. The Feddes water stress function as resulting from our parameterization shows a higher drought resistance of Norway spruce than previously assumed (Figure 7c).

3.4. Simulated Water Balance before and after Calibration

The mean annual simulated water balance of the two research plots indicates generally higher transpiration at the Riparian plot compared to the Slope plot (

Table 3. Mean annual simulated water balance (October 2009–September 2014) before and after calibration of the Feddes parameters in HYDRUS-1D. Percentage change of calibrated water balance component in comparison to un-calibrated value is given in brackets. ET₀: potential evapotranspiration, ET_act: actual evapotranspiration; all details in mm year⁻¹.

Water Balance Component	Riparian		Slope
Precipitation	1321
ET0	597
	un-calibrated	calibrated	un-calibrated	calibrated
ETact	597	597	403	445 (+10%)
Surface runoff	54	54	0	0
Drainage	634	634	841	797 (−5%)

). While the water balance of the Riparian plot was not affected by the calibration of the Feddes parameters, the calibration raised ET_act by 10% for the Slope plot. At the same time, drainage decreased by 5%. For all simulation runs, we achieved a very good fit between mean daily sap flow and simulated actual transpiration (R² = 0.82 at Slope and 0.84 at Riparian plot). At the Slope plot, the correlation between sap flow data and simulated actual transpiration was slightly improved by the calibration of the Feddes parameters (R² = 0.83). For the Riparian plot, R² was not affected from the calibration.

3.5. Sensitivity of Feddes Parameters to Root Zone Variation and SoilNet Measuring Location

Variations of the potential root water uptake factor’s depth distribution did not lead to significant variations of the simulated root zone h (Mann-Whitney U-test). The same applies for variation of the rooting depths.

The difference between the mean soil water supply of the surrounding SoilNet measuring locations and that of the stations we selected for our study (SN17 and SN28; cf. Figure 1) was statistically significant (p < 0.001). However, the absolute difference among the simulated SoilNet locations was small and the contrasting soil moisture regimes between the two plots are still clearly visible (Figure 8). Respectively, the resulting Feddes parameters were only slightly affected by the selection of the SoilNet measuring location. Parameter h_3,low as determined from the mean root zone water potentials of the surrounding SoilNet locations was −3800 cm, which is 300 cm below the h_3,low value we determined from SN28. h₄ decreased from −15,000 cm (SN28) to −13,600 cm pressure head (mean of the surrounding SoilNet stations), while h₂ was not affected from the choice of the SoilNet measuring location. This indicates that soil type and elevation were reasonable criteria for the choice of our measuring locations.

4. Discussion

4.1. Soil Moisture Simulations

The water balance simulation captured the prevailing soil moisture dynamics reasonably well. Considering the measurement uncertainty of the installed SWC sensors (4–5 vol %; [37]) the RMSE of 3.4 to 5.3 vol % appear very low and are within the range of other HYDRUS-1D studies on the Wüstebach test site [50,53]. The varying soil hydraulic properties among plots and layers (Table 2) fit well to the observations of [36,38]. The low K_s value in the riparian zone is plausible, because the soils in the valley bottom have a higher bulk density and a lower macro-porosity than the soils on the hillslopes [37,38].

The simulated mean annual root-zone water potentials (Figure 6) highlight the contrasting soil moisture regimes of our research plots (wet and dry). These characteristics become even more pronounced when referring only to the vegetation periods (May-September), where mean root-zone pressure heads are −124 (Riparian) and −4308 cm (Slope). Covering particularly the stages of critical soil water supply (undersupply/oversupply) our data seems very suitable for investigating the on-set of drought and aeration stress in plants.

4.2. Water Stress Response of Norway Spruce

There is evidence that the shallow rooting Norway spruce is more vulnerable to drought than other species with a deeper rooting system [4,23,69]. Also, [4,25] found the sap flow activity of Norway spruce to decouple from climate variables with increasing soil water shortage. However, critical levels of soil water supply that induce a decline in sap flow rates and transpiration considerably vary among studies.

While [4] observed gradually decreasing sap flow rates with increasing soil water shortage for REW < 0.4 (which corresponds to a pressure head of −3240 cm at our Slope plot), [23] found Norway spruce to be more drought resistant. Depending on the soil depth considered, they report mean critical REW values of between 0.21 (0–20 cm) and 0.24 (0–50 cm). However, the observed critical REW level strongly varied among trees (REW: 0.18–0.37), which [23] attribute to the heterogeneous soil conditions on site. This is in line with [4], who found the stress response of Norway spruce to become increasingly variable between individuals with continuing drought stress. Comparing the water stress response of Norway spruce among more than 60 pure and mixed stands, [24] report difficulties in determining a consistent critical level of REW_. For most of the plots, drought limitation of sap flow began around a REW of 0.3. However, the data scatter around this critical limit was large and a clear decrease in transpiration for all stands except peat stands was only observed for REW below 0.25.

In simulation studies, a REW level of 0.3–0.4 is commonly used to indicate the on-set of drought stress in Norway spruce [31,61,62]. However, this critical REW range is widely accepted for many species (not only for spruce) and our data show that REW does not represent the actual energy state of the soil. For example, at the Slope plot a REW value of 0.3 corresponds to a root-zone pressure head of −4680 cm, which is in the same range (pF 3.7) as the critical h value of −4100 cm (pF 3.6) we observed in our study (Figure 7a). However, a REW threshold of 0.3 at the Riparian corresponds to a root-zone pressure head of −2140 cm suggesting that drought stress would occur much easier at this site (i.e., at lower soil matric potentials). On the other hand, a critical root-zone water potential of −4100 cm pressure head would here correspond to a REW of 0.17.

One could argue that the poor fit between our critical h and the critical reported REW values for the Riparian plot is of minor importance, because critical water contents are never reached [24,70]. Moreover, the concept of REW seems to apply fairly well to the conditions at the Slope plot and to other soils in central European high and low mountain ranges [4,31,61,62]. However, Norway spruce is also widely distributed in boreal forest landscapes, where typical critical REW values seem to be lower than 0.3 [23,24] and soils are often peaty [71], hence showing similar characteristics (low residual water content, high saturated water content; [72,73,74]) as our plot in the riparian zone. Our data indicates that for such soils the concept of REW is likely to fail. This emphasizes the value of soil water potential in comparison to other measures of soil water supply, because, respective water stress thresholds can be applied to any type of soil as long as soil water retention and plant characteristics are known.

However, while critical soil water potentials have been determined for many agricultural crops, little effort has been put into the quantification of Feddes parameters for forest trees. To our knowledge, the concept of transpiration limiting soil water potentials has hardly been applied in spruce site simulation studies. In the few studies we found, the drought stress threshold was generally assumed to be less negative than the value we determined from our sap flow data (h_3,low = −4100 cm).

Huang et al. [75] set h₃ to the root-zone h at a SWC of 27.5% of the SWC at field capacity. This corresponds to the mean critical SWC they reviewed from different studies on drought-stress in evergreen tree species. Applied to our data and considering field capacity to be represented by SWC at −300 cm pressure head (see above), this threshold corresponds to a pressure head of −3890 cm at the Riparian plot, but to h < −10 million cm at the Slope site, which is far below the common used value for the permanent wilting point. Dependent on the soil type, Rosenqvist et al. [76] assumed the onset of drought stress for Norway spruce at water potentials of −500 (coarse sandy soils) to −1000 cm pressure head (clayey soils). However, these texture dependent water stress thresholds contradict the idea that the soil water potential is already a function of SWC and soil properties. From a plant physiological perspective, species-specific water stress thresholds should be used instead. Vogel et al. [59] set their drought stress thresholds to −1200 (h_3,low) and −600 cm pressure head (h_3,high). These parameters were derived from other hydrological studies that were conducted on similar soils, but had no plant physiological basis. Nevertheless, applying the parameter set of [59] to our simulations, we achieved a very good fit between measured and simulated SWC (RMSE < 5.3 vol %). Also, mean daily sap flow and simulated actual transpiration showed high correlations (R² > 0.82). However, although the adaptation of the Feddes function had little impact on the general dynamic of the simulated fluxes, the adjustment of h_{3, low} significantly changed transpiration (p < 0.001) and drainage fluxes (p < 0.01) at the Slope plot (two-sided Mann-Whitney U-test; Table 3) and made it more realistic: Based on eddy-covariance measurements, [40] found the annual ET_act in the Wüstebach catchment to cover approximately 90% of ET₀. In our Slope plot simulations this ratio was 80% for the uncalibrated Feddes model, but 87% for the calibrated simulation run in the samestudy period (2010–2013). The runoff ratio (runoff/precipitation) observed by [40] was 56%, while that in our simulations ranged between 62% (Slope plot after calibration of the Feddes model) and 65% (Slope plot before calibration of the Feddes model). This result confirms: (1) the finding of [77] that a good fit between simulated and measured SWC does not necessarily imply that the simulated water fluxes are correct; and (2) the assumption of [78] that sap flow measurements can help improving root water uptake simulations.

It has to be noted that the impacts of the Feddes model on simulated water fluxes depend on local site conditions. The well-watered Riparian plot was not affected by the calibration, because critical pressure heads were never reached. Drought stress simply did not occur and although there is evidence on the vulnerability of Norway spruce to waterlogging [4,79], our data did not indicate arising aeration stress in the riparian zone (Figure 7b). One reason could be that spruce indeed only suffers from aeration stress at fully saturated soils (h₁ = h₂ = 0), which does not apply to the vegetation periods in this study (Figure 5). Another possible explanation is that the transpiration limiting effects of aeration stress occur on longer time-scales than the daily scale we were operating at.

Our data show that soil water supply strongly influences the transpiration of Norway spruce. Due to the observed drought stress on the hillslope (Figure 7a), the simulated actual transpiration at the Slope plot was 13% lower than that of the riparian zone (Table 3). This result highlights the importance of considering soil moisture patterns in the modeling of transpiration fluxes from the forest plot to the catchment scale.

5. Conclusions

In this study, we confirm the potential of sap flow measurements for the determination of water stress in forest trees. By combining soil hydrological simulations with a water stress factor derived from sap flow data, we were able to parameterize the Feddes water stress function for Norway spruce and improve our water balance simulations on site. Our results show that small-scale variations in soil water supply significantly influence the water balance of spruce sites. Considering soil moisture patterns in the model setup can thus improve the simulation of transpiration fluxes on the catchment scale.

Although additional sap flow stations and a bigger sample size per plot would have provided a more detailed view on the conditions on site, the sampling concept (wet plot against dry plot) delivered valuable insights into the water stress response of Norway spruce. The advantage of our sampling approach is that an upscaling of sap flow density to the tree and plot scale is unnecessary. Hence, uncertainties related to the upscaling procedure (cf. e.g., [7,66,67,68,69,70]) can be avoided. However, to transfer our finding to other sites and settings, more research on species-specific feedbacks in the soil-vegetation-atmosphere system is needed. There is still a lack of studies determining the critical limits of soil water supply for trees. Since forests play a vital role in regional and trans-regional water cycles and, against the background of climate change, the assessment of respective fluxes becomes increasingly important, future research should focus on: (1) more species-specific investigations on the water stress response of trees; and (2) an improved assessment of small-scale variabilities in soil water supply. With our sampling strategy, we present a simple and cost-efficient approach to achieve these goals, even with a limited number of sample trees.