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

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.


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 Figure 1. Overview of the Wüstebach test site and its location in Germany. The research plots (Riparian, Slope) cover differing soil types and topographic situations, which results in contrasting soil moisture regimes (wet and dry). The grey dots indicate the SoilNet sensor network of which two stations have been selected for this study (yellow dots). Black dots show the SoilNet stations we selected to analyze soil moisture variability around the sap flow stations. The difference between the displayed contour lines is 2.5 m.

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]: (1) Figure 1. Overview of the Wüstebach test site and its location in Germany. The research plots (Riparian, Slope) cover differing soil types and topographic situations, which results in contrasting soil moisture regimes (wet and dry). The grey dots indicate the SoilNet sensor network of which two stations have been selected for this study (yellow dots). Black dots show the SoilNet stations we selected to analyze soil moisture variability around the sap flow stations. The difference between the displayed contour lines is 2.5 m.

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]: (1) in which SFD is the sap flux density (g·m −2 ·s −1 ), ∆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 2 ) 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 catchmentwide 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 0 was calculated using the FAO-Penman-Monteith method [45] following [40], who first calculated hourly ET 0 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 0 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 0 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.

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]: where where h is the pressure head (cm), θ r is the residual water content (cm 3 cm −3 ), θ s the saturated water content (cm 3 cm −3 ), α 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 −1 ), 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 0 ) and potential transpiration (T 0 ) 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 0 and T 0 was of minor importance. However, applying our approach to deciduous trees with seasonal LAI dynamics would require the consideration of the T 0 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]: in which q(h) (cm day −1 ) is the discharge rate, h (cm) is the pressure head at the bottom of the soil profile, A qh (cm day −1 ) and B qh (cm −1 ) 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 (Tact) and potential (maximum possible) transpiration (T0) and thus ranges from 0 (maximum stress/no root water uptake) to 1 (no stress/optimum root water uptake). Optimum root water uptake (T0 = Tact) is given between h2 (field capacity FC) and h3 (onset of drought stress). The decreasing plant water stress between h2 and h1 considers that plants typically experience increasing aeration stress towards saturated conditions; h4 corresponds to the permanent wilting point (WP), h1 to the point of maximum aeration stress (anoxic conditions).
Presuming that h3 decreases with decreasing T0, HYDRUS-1D allows for making h3 a function of T0. The range of T0, within which h3 becomes a function of T0 may be defined by the user, but is usually 1-5 mm day −1 (default values). The h3 value for T0 < 1 mm day −1 is called h3,low, that for T0 > 5 mm day −1 h3,high. If T0 is between 1 and 5 mm day −1 and the soil dries out to less than h3,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 T0 [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 h1 = h2 = 0 cm, h3,high = −600 cm, h3,low = −1200 cm, and h4 = −15,000 cm. Later, we manually optimized the Feddes parameters on the basis of our sap flow data (see below). 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 0 ) and thus ranges from 0 (maximum stress/no root water uptake) to 1 (no stress/optimum root water uptake). Optimum root water uptake (T 0 = T act ) is given between h 2 (field capacity FC) and h 3 (onset of drought stress). The decreasing plant water stress between h 2 and h 1 considers that plants typically experience increasing aeration stress towards saturated conditions; h 4 corresponds to the permanent wilting point (WP), h 1 to the point of maximum aeration stress (anoxic conditions).
Presuming that h 3 decreases with decreasing T 0 , HYDRUS-1D allows for making h 3 a function of T 0 . The range of T 0 , within which h 3 becomes a function of T 0 may be defined by the user, but is usually 1-5 mm day −1 (default values). The h 3 value for T 0 < 1 mm day −1 is called h 3,low , that for T 0 > 5 mm day −1 h 3,high . If T 0 is between 1 and 5 mm day −1 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 0 [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 1 = h 2 = 0 cm, h 3,high = −600 cm, h 3,low = −1200 cm, and h 4 = −15,000 cm. Later, we manually optimized the Feddes parameters on the basis of our sap flow data (see below).

Figure 3.
Feddes plant water stress function modified after [28]. Root water uptake commences with the onset of desaturation at h1 (anoxic moisture conditions), optimum root water uptake is given between h2 (field capacity, FC) and h3 (onset of drought-stress), where indiceshigh and low refer to the magnitude of potential transpiration rates. From optimum root water uptake towards h4 (permanent wilting point, WP) the plants experience increasing drought stress.

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.

Parameterizing the Feddes Model Using Sap Flow Data
In the classical Feddes approach, plant water stress is derived from the ratio of Tact and T0 (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 (SFDn,rip for Riparian and SFDn,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 SFDdrought (1: optimum water supply; 0: maximum water stress) for the Slope plot by dividing SFDn,sl by SFDn,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 SFDdrought against the simulated Slope water potential and analyzed the correlation between SFDdrought and root zone h for different levels of soil water shortage. The data which showed the highest significant correlation (R²  [28]. Root water uptake commences with the onset of desaturation at h 1 (anoxic moisture conditions), optimum root water uptake is given between h 2 (field capacity, FC) and h 3 (onset of drought-stress), where indices high and low refer to the magnitude of potential transpiration rates. From optimum root water uptake towards h 4 (permanent wilting point, WP) the plants experience increasing drought stress.

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.

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 0 (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 2 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 4 .
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 0 range at which h 3 is a function of T 0 . 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 0 value of 2.7 mm day −1 can be used as lower boundary for the h 3 -dependent T 0 range. Furthermore, we used T 0 at the day of maximum observed daily SFD (4.8 mm day −1 ) as the upper boundary for the h 3 -dependent T 0 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 1 and h 2 , we plotted SFD aeration against the simulated Riparian root-zone water potential and analyzed the observed data trend as described for SFD drought .

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]: 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 .

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.

Sap Flow Data
The observed SFD dynamic of both research plots corresponds well to the dynamic of the calculated ET 0 (Figure 4). R 2 between mean daily SFD and ET 0 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), 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).
Water 2017, 9, x FOR PEER REVIEW 9 of 19 ( 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.

Sap Flow Data
The observed SFD dynamic of both research plots corresponds well to the dynamic of the calculated ET0 (Figure 4). R² between mean daily SFD and ET0 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), 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).

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

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 Ks decrease with depth ( Table 2).
After inverse modelling, and manual calibration of Aqh and Bqh 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).
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).  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, Ks saturated hydraulic conductivity (cm day −1 ).  Norway Spruce (h1, h2, h3,low, h4) 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 SFDdrought and simulated root-zone h was observed below −4100 cm (R² = 0.77, p < 0.001). Thus, we set h3,low = −4100 cm. Extrapolating the trend-line between SFDdrought and simulated root-zone h towards the x-axis, where SFDdrought = 0 (maximum water stress), resulted in an h4 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 h1 and h2 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).  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 −1 ).

Feddes Parameters for Norway Spruce (h 1 , h 2 , h 3,low , h 4 )
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 2 = 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 4 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 1 and h 2 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).

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). While the water balance of the Riparian plot was not affected by the calibration of the Feddes parameters, the calibration raised ETact 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. 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. ET0: potential evapotranspiration, ETact: actual evapotranspiration; all details in mm year −1 .

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). 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 2 = 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 2 = 0.83). For the Riparian plot, R 2 was not affected from the calibration. 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 0 : potential evapotranspiration, ET act : actual evapotranspiration; all details in mm year −1 .

Water Balance Component
Riparian Slope

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 4 decreased from −15,000 cm (SN28) to −13,600 cm pressure head (mean of the surrounding SoilNet stations), while h 2 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.

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 h3,low as determined from the mean root zone water potentials of the surrounding SoilNet locations was −3800 cm, which is 300 cm below the h3,low value we determined from SN28. h4 decreased from −15,000 cm (SN28) to −13,600 cm pressure head (mean of the surrounding SoilNet stations), while h2 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.  Figure 1). Simulated root zone water potentials of SN17 (Riparian) and SN28 (Slope) are shown in black. Note that the pressure head (-cm) is plotted on the log-scale (single positive values cannot be shown).

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 Ks 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.  Figure 1). Simulated root zone water potentials of SN17 (Riparian) and SN28 (Slope) are shown in black. Note that the pressure head (-cm) is plotted on the log-scale (single positive values cannot be shown).

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.

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 3 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 2 > 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 0 . 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 1 = h 2 = 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.

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.