When describing the movement of water in a variably saturated plant root zone, most existing hydrological models use the assumption of quasi-steady-state flow to relate root water uptake to canopy transpiration, thereby neglecting the effect of changing plant water storage. This approach is known to be problematic, especially when considering relatively large volumes of water stored in the tissues of tall trees. We propose a simple algorithm, based on the concept of whole-plant hydraulic capacitance, to deal with the problem. The algorithm is implemented in a one-dimensional soil water flow model involving vertically distributed macroscopic root water uptake. In this study, the proposed transient storage approach was compared with the quasi-steady-state approach. Both approaches were used to simulate soil water flow and diurnal variations of transpiration at a forest site covered with Norway spruce [Picea abies (L.) H. Karst.]. The key parameter of the transient storage approach, the plant hydraulic capacitance, is estimated by comparing the variations in the potential transpiration rate, derived from micrometeorological measurements, with observed sap flow intensities. The application of the proposed algorithm leads to more realistic predictions of root water uptake rates at the site of interest. The inclusion of the plant water storage effects improved the ability of the model to capture the anticipated diurnal variations in actual transpiration rates. The algorithm can be easily implemented into existing soil water flow models and used to simulate transpiration stream responses to varying atmospheric and soil moisture conditions including isohydric and partly also anisohydric plant responses to drought stress.

RWU, root water uptake

The role of plant water storage in soil–plant–atmosphere systems has been given increasing attention in recent plant physiology studies (e.g., Borchert and Pockman, 2005; Cermak et al., 2007; Meinzer et al., 2008; Carrasco et al., 2015; Pfautsch et al., 2015). These studies have shown that water can be stored intracellularly or extracellularly in plant tissues. Stored water can be released from plant tissues due to their elasticity or by cavitation. According to Cermak et al. (2007), the amount of water released from water storage compartments inside mature trees can account for up to one third of daily water loss by transpiration. Stem water storage effects result in significant buffering of xylem water potential fluctuations, which helps to preserve the integrity of the xylem water system during water stress periods (Meinzer et al., 2003; Scholz et al., 2007).

To minimize the risk of xylem embolism, caused by cavitation, plants need to keep their xylem water potential above a certain threshold value specific for different plant tissues and plant species (Cochard, 1992; Lu et al., 1996). For plants exposed to water stress, the range of xylem water potential between the daily minimum potential and the potential corresponding to the onset of cavitation can be interpreted as a safety margin (Sperry, 2000; Meinzer et al., 2009; Carrasco et al., 2015). It is assumed that there are two possible regulation strategies in dealing with water stress—referred to as isohydric and anisohydric behaviors (Tardieu and Simonneau, 1998; McDowell et al., 2008; Pretzsch et al., 2013). Isohydric plants maintain their midday xylem water potential at a relatively constant level by closing their stomata and thus reducing transpiration. By contrast, anisohydric plants allow the midday xylem water potential to decline and maintain transpiration until the transpiration stream is disrupted by cavitation.

An increasing number of soil–plant–atmosphere continuum models include a capacitance term to account for the water storage in plant tissues. Among the more recent ones, for example, are the models of Lhomme et al. (2001), Zhuang et al. (2014), and Tardieu et al. (2015). While these models pay a great deal of attention to hydraulic resistances and capacitances of plant tissues along the transpiration stream (Lhomme et al., 2001); capacitance, inductance, and fuse effects associated with plant xylem water storages, hydraulic architecture of leaf system and imperfections in soil–root contact (Zhuang et al., 2014); and coordination of stomatal controls, transpiration forcing, leaf growth, and abscisic acid signaling (Tardieu et al., 2015), they neglect the effects of distributed soil water status and treat the soil compartment as a zero-dimensional system, characterized by effective values of soil water potential and soil water content.

Plant responses to water stress are reflected in variations of transpiration and root water uptake intensities, which are important hydrological variables. It is essential, therefore, that these responses are adequately parameterized and properly incorporated in hydrological models. In the context of vadose zone hydrology, root water uptake (RWU) is usually represented as a spatially distributed macroscopic sink of soil water. Most existing soil water flow models apply the assumption of quasi-steady-state flow in plant xylem to relate the RWU intensity to the actual transpiration rate, neglecting transient water storage effects in plant tissues (Raats, 2007; Šimůnek and Hopmans, 2009; Jarvis, 2011; de Willigen et al., 2012; Guswa, 2012; Javaux et al., 2013; de Jong van Lier et al., 2013; Couvreur et al., 2014; Manoli et al., 2014). A few models combine soil water and plant xylem flow in a coupled system in which both the soil and the xylem are treated as porous media with distributed hydraulic properties (Janott et al., 2011; Rings et al., 2013). Such an approach is, in principle, much more realistic, but the resulting models are complex and computationally demanding.

In this study, we combined the concept of lumped hydraulic capacitance, representing the hydraulic effect of variable whole-plant water storage, with the concept of spatially distributed macroscopic root water uptake to obtain an algorithmically simple transpiration stream model that could be easily implemented into existing soil water flow models. The proposed approach allows simulation of basic plant responses to varying atmospheric and soil moisture conditions, including isohydric–anisohydric adjustments to water stress. The suggested approach was tested by comparing the simulation results obtained with the transient plant water storage approach with those based on the quasi-steady-state assumption and by comparing both approaches with observed Norway spruce sap flow data.

Study Site

The Norway spruce sap flow data were acquired in a small mountainous watershed (the Liz catchment), situated in the Bohemian Forest, Czech Republic. The catchment outlet is located 830 m above sea level. The average annual precipitation is 861 mm, and the average annual temperature is 6.5°C. The prevailing soil type is oligotrophic forest Eutric Cambisol developed on biotite paragneiss bedrock. The catchment is mostly covered with Norway spruce and European beech (Fagus sylvatica L.). The majority of spruce trees are 94 yr old, with a height of about 28 m and a diameter of about 40 cm. Data used in the present study were recorded at two monitoring stations—at the spruce forest floor and at a nearby meadow (Tesar et al., 2006). Variables measured at the forest station included sap flow, air and soil temperatures, soil water pressure, soil water content, and precipitation. Complementary micrometeorological variables (Fig. 1) were monitored at the meadow station. Sap flow in the stem xylem of four Norway spruce trees was monitored at breast height using the heat field deformation method with multipoint sensors (Nadezhdina et al., 2012).

Soil Water Flow Model

To describe the vertical flow of water in a soil profile, we use the one-dimensional Richards’ equation (Richards, 1931):
where θ is the volumetric soil water content (m3 m−3), K is the unsaturated hydraulic conductivity of the soil (m s−1), S is a sink term representing the intensity of root water uptake (s−1), H is the soil water potential (m), h is the soil water pressure head (m), and z is the vertical coordinate, considered positive upward (m). The intensity of root water uptake varies with depth and time according to the transpiration demand of the atmosphere, the soil water availability, and the plant water status. The algorithm used to evaluate the RWU intensity is often referred to as the macroscopic RWU model (e.g., Raats, 2007).

Root Water Uptake Model

The uptake of soil water by plant roots is described by a macroscopic RWU model based on the traditional water-potential-gradient approach (van den Honert, 1948; Hillel et al., 1976; and many others). It is assumed that the RWU intensity is directly proportional to the difference between water potentials of the soil and the root xylem and indirectly proportional to the respective hydraulic resistances:
where σ is the specific active root surface (m−1), rroot is the root radial resistance (s), rsoil is the effective hydraulic resistance of the soil matrix (s), Hrx is the root xylem water potential and Hsoil is the soil water potential. In this study, we assumed that Hrx varies with time but is constant within the root system. This is consistent with the observation that the axial resistance to water flow in roots is one to several orders of magnitude lower than the radial resistance (e.g., Steudle and Peterson, 1998).
We assumed that the effective soil hydraulic resistance can be expressed as
where λ(z) is the characteristic length associated with the transport of water from the bulk soil to the root surface, derived from the active root system geometry (e.g., Vogel et al., 2013, 2016).
The active root system geometry is characterized by the active root length density R (m−2), the specific active root surface σ, the average active root radius r0 (m), and the rhizosphere radius r1 (m). In this study, the spatial variability of R and σ was limited to the vertical direction and r0 was spatially invariant. Furthermore, assuming the cylindrical character of the root system geometry, the vertical distributions of R, σ, and r1 can be related through simple formulas:
and, finally, the characteristic length λ can be expressed as a fraction of r1:
where a is the dimensionless fraction coefficient.

Plant Water Balance

Under the assumption of quasi-steady-state flow of water through the soil–plant–atmosphere continuum, the RWU intensity, integrated over the depth of the root zone, is equal to the actual transpiration rate Ta (m s−1):
where Tr is the whole-plant root water uptake (m s−1) and zR and z0 are the coordinates of the lower and upper boundary of the root zone, respectively (m).
A more realistic form of the whole-plant water balance equation can be obtained by adding a transient storage term to Eq. [6]:
where W is the whole-plant water storage (m). Assuming that W can be expressed as a single-valued function of Hrx, we obtain
where C is the whole-plant hydraulic capacitance. Because both W and Hrx are expressed in units of length, C is here defined as dimensionless.
Substitution of Eq. [8] into Eq. [7] yields the final form of the whole-plant water balance equation:

Determination of Plant Hydraulic Capacitance

A variety of methods have been used to determine the hydraulic capacitance of plants. The capacitance of sapwood can be estimated directly from measurements of sapwood water release curves (e.g., Meinzer et al., 2003; Scholz et al., 2007). Vergeynst et al. (2014) estimated whole-branch capacitance by combining measurements of acoustic emissions caused by cavitation, supplemented with branch diameter shrinkage and gravimetric water loss measurements. Whole-tree capacitance can be determined from observations of the time lag between evaporative demand and stem sap flow (e.g., Phillips et al., 1997; Kumagai et al., 2009) or the time lag between sap flow signals measured at different heights of a tree stem (e.g., Cermak et al., 2007).

Sapwood capacitance has been defined as the slope of a line fitted to the initial phase of the sapwood water release curve (Meinzer et al., 2003, 2008; Vergeynst et al., 2014). According to Meinzer et al. (2003), the measured water release curves are nearly linear in a relatively broad range of xylem water potentials, changing the slope as the water potential approaches the point of cavitation. Thus, in the linear range of water release curves, the sapwood capacitance can be approximated by a constant value. The capacitance value is expected to significantly increase at the onset of cavitation.

When simultaneous estimates of stem water storage and root xylem water potentials are available, it is possible to determine the effective whole-plant hydraulic capacitance directly from the definition of Eq. [8]:

In this study, we estimated the whole-tree hydraulic capacitance of Norway spruce by comparing the potential transpiration rate (determined from micrometeorological measurements) with the observed time series of stem sap flow data. On clear summer days, the observed course of potential transpiration has a nearly semi-sinusoidal shape, reaching zero at sunset, while the trunk xylem sap flow exhibits an after-dusk tailing shape (Fig. 2). We assumed that Hrx increases overnight from its pre-dusk value to the value that is close to equilibrium with the average soil water potential. The C estimate can be obtained by dividing the cumulative nocturnal root water uptake (interpreted as ΔW) by the corresponding difference in Hrx. Because the pre-dusk value of Hrx was not measured at our study site, we used the combined soil water flow and transpiration stream model iteratively to estimate both C and ΔHrx (setting C = 0 as an initial estimate).

Isohydric Control

In isohydric plants, the xylem water potential is limited to values above a certain threshold value, here denoted as Hcrit. The plants maintain a distinct safety margin between this value and another threshold value associated with the onset of cavitation, Hcav, by means of stomatal control (e.g., Meinzer et al., 2009).

To simulate isohydric behavior, we used a simple procedure that ensures that the root xylem water potential, Hrx, is kept greater than or equal to the threshold value Hcrit. The procedure consists of two alternative stages: (i) under favorable atmospheric and soil moisture conditions, the actual transpiration equals the atmospheric demand, i.e., Ta = Tp (where Tp is the potential transpiration rate [m s−1]), and the root xylem potential is greater than the critical value, i.e., Hrx > Hcrit; and (ii) when the atmospheric demand is too high or the soil profile too dry, the xylem water potential reaches the critical value (Hrx = Hcrit) and the actual transpiration becomes a fraction of its potential value (Ta < Tp).

Transpiration Stream Model

The proposed transpiration stream model requires a simultaneous solution of soil water flow Eq. [1], root water uptake Eq. [2], and plant water balance Eq. [9]. This is achieved by a simple sequential coupling procedure.

The time derivative on the left-hand side of Eq. [9] is discretized using a finite difference formula:
By substituting Eq. [2] and [11] into Eq. [9], we obtain
This equation is then used to derive an explicit formula for the root xylem water potential as a function of transpiration, as well as the inverse relationship, i.e., transpiration as a function of root xylem water potential (for a given vertical distribution of the soil water potential Hsoil):

The derivation of Eq. [13] and [14] requires that Hrx is constant within the root system and changes only with time. In this study, we used several additional simplifying assumptions, such as that rroot is constant, the root geometry is time invariant, and C is constant. However, these assumptions, unlike the one about the spatial invariance of Hrx, can be easily relaxed.

At each time step of the numerical solution of Richards’ Eq. [1], the root xylem water potential, Hrx, is first evaluated from Eq. [13], assuming Ta = Tp. If Hrx falls below Hcrit for a given potential transpiration rate, Hrx is set equal to Hcrit and Ta is calculated using Eq. [14], yielding Ta < Tp. The resulting Hrx value is then used to calculate the RWU intensity S(z) (from Eq. [2]) and to update the value of the sink term in Eq. [1].

Quasi-Steady-State Flow Assumption

The quasi-steady state flow condition can be invoked by setting C equal to zero in Eq. [13] and [14]. This leads to an algorithm that is similar to the one suggested by Hillel et al. (1976) and later implemented in several RWU models (e.g., de Jong van Lier et al., 2008; de Willigen et al., 2012; Vogel et al., 2013; to name the more recent ones).

Anisohydric Control

To simulate anisohydric behavior, Hcrit in the above algorithm can be set equal to Hcav, and the plant hydraulic capacitance at Hcav modified to represent a partly embolized xylem. For short periods, partial cavitation could be a reversible process. However, if unfavorable weather conditions induce irreversible changes in plant tissues, hydraulic capacitance becomes a complex function of related plant physiological processes. In such a case, a more elaborate transpiration stream model would have to be considered.

Model Application at the Study Site

The proposed transpiration stream algorithm was implemented in the numerical code S1D (Vogel et al., 2010, 2013), which is designed to simulate soil water flow in a vertical soil profile containing preferential pathways (solving a dual set of Richards’ equations for a two-domain flow system). The soil profile at the study site was assumed to consist of two hydraulically connected flow domains—the soil matrix domain and the preferential flow domain. By partly diverting storm water from the root zone, preferential flow significantly affects the soil water balance at the macroscopic scale; therefore it could not be neglected. It was assumed that roots take water up from the soil matrix, i.e., the S(z) function was evaluated in the soil matrix domain only.

The S1D code was used to simulate soil water flow and spruce tree transpiration at the Liz site from April to September 2010. The soil profile was schematized as a vertical one-dimensional 70-cm-deep soil column. The soil matrix flow domain was divided into two layers—an organic soil layer and a mineral soil layer—distinguished by different soil hydraulic properties (Table 1). The initial condition for the simulations was derived from the observed variations in soil water pressure at the beginning of the simulated period. Precipitation data, measured directly on site, were used to define the upper boundary condition. The observed soil water pressure was applied as a time-dependent Dirichlet type boundary condition at the lower boundary of the soil column. Potential transpiration of spruce trees was estimated by means of the Penman–Monteith equation (Monteith, 1981) based on the micrometeorological data observed at Liz. This was done by applying the standard FAO procedure for reference evapotranspiration (Allen et al., 1998), using hourly data from the nearby meteostation (located at the meadow site). The actual transpiration rate was assumed to drop to zero overnight as well as during rainfall events.

The values of plant-related parameters (with the exception of the whole-plant hydraulic capacitance) were adopted from our previous study (Vogel et al., 2016). The parameter values are summarized in Table 2. We assumed that water is taken up mostly by fine roots (commonly defined as roots with a diameter <0.2 cm). Based on values reported in the literature for spruce trees (e.g., Ruess et al., 2006; Kucbel et al., 2011), the mean radius of active roots, r0, was set equal to 0.01 cm and the maximum root length density to 0.2 cm−2. The shape of the root length density function, R(z), was estimated based on information found in the literature and a limited survey performed at the site of interest. We assumed a constant value of R(z) equal to the estimated maximum root length density for depths 1 to 20 cm, followed by a linear decrease, with R(z) equal to zero at the lower boundary of the root zone at the depth of 70 cm. No roots were assumed in the uppermost 1 cm of organic soil. Furthermore, we estimated that the coefficient a in Eq. [5] was equal to 0.05. This value was determined using the mesoscopic analysis of radial flow toward roots (Vogel et al., 2016).

To facilitate the comparison with simulated actual transpiration rates (m s−1), the time series of the sap flow densities (kg m−2 s−1) obtained for the individual Norway spruce tree specimen were upscaled to the forest stand using a procedure explained in our previous study (Vogel et al., 2016). The threshold value at the onset of cavitation, Hcav, was reported around −250 m for Norway spruce (estimated from the leaf water potential by Cochard [1992] and Lu et al. [1996]). Assuming a safety margin of about 50 m (Bréda et al., 2006) and the difference between leaf and root water potentials, the Hcrit value can be expected to be in the range of about −200 to −150 m. In this study, we assumed Hcrit = −150 m.

The whole-tree hydraulic capacitance was set equal to 9 × 10−6. The value was chosen to produce a reasonable agreement between the simulated variation of Tr and the observed sap flow. Expressed in units of water mass per unit pressure change, it corresponds to 12 kg MPa−1 per tree (assuming 13 m2 of ground area per tree, as determined for the study site), or 7 kg MPa−1 m−3 of sapwood (assuming 1.6 m3 of sapwood per tree).

Effects of Plant Water Storage

Precipitation data, together with the key micrometeorological variables affecting the magnitude of spruce tree transpiration (Fig. 1a), indicate that the soil was well supplied with water throughout the season studied, with two drier periods in the first half of July and the second half of August. That is also reflected in the soil water pressure (Fig. 1b) and soil water content data (not shown).

The effects of two different approaches to plant water storage in coupled soil-water-flow and root-water-uptake modeling are presented in Fig. 2. A 12-d period with intensive transpiration demand, interrupted by a rainy day, was selected to illustrate the functioning of the model.

In Fig. 2a, the actual transpiration rate, Ta, calculated assuming quasi-steady-state flow through the transpiration stream, is presented together with the applied potential transpiration, Tp, and the observed sap flow in a selected spruce tree specimen. Under the quasi-steady-state assumption, the whole-plant RWU rate equals the actual transpiration rate (the plant water storage is constant) and both are nonzero only when nonzero potential transpiration is applied (i.e., during the daytime in our simulation). In contrast, the sap flow measurements indicate nonzero nocturnal fluxes in the stem xylem.

Simulation results obtained assuming transient plant water storage are presented in Fig. 2b. Allowing the plant water storage to change causes separation of the actual transpiration rate from the actual RWU rate. The predicted sunny-day transpiration stream variations can be divided into three stages. Starting in the morning, the plant actual transpiration fully meets the potential transpiration demand while the root water uptake rate lags behind and the plant water storage decreases. When the root xylem water potential reaches the critical value, Hcrit, the plant water storage becomes fixed and the plant transpiration becomes limited by the plant RWU. After the potential transpiration rate reduces below the limiting RWU rate, the actual transpiration follows the potential transpiration demand, while the RWU rate reduces gradually, allowing the nocturnal replenishment of the plant water storage.

In Fig. 2c, the simulated root xylem water potential variations, assuming either constant or transient plant water storage, are presented. The value of Hrx fluctuates between the soil water potential (depicted by dashed lines) and the critical xylem potential maintained by the plant (set to −150 m). Considering nonzero hydraulic capacitance helps to smooth the xylem potential variations and to delay or even eliminate the critical state at which the transpiration becomes reduced.

Our model assumes a linear relationship between the root xylem water potential, Hrx, and the whole-plant water storage, W. Thus the Hrx variations directly reflect the development of plant water storage depletion (note the secondary vertical axis on the right-hand side of Fig. 2c). Assuming that W is full at Hrx equal to zero, the maximum possible depletion of W is determined by the critical value of the root xylem water potential, Hcrit, and by the value of plant hydraulic capacitance. In our case, the maximum reduction is about 1.4 mm. Considering 13 m2 of ground area per tree, the maximum possible reduction in the xylem water storage is 18 L per tree. Cermak et al. (2007) detected a maximum daily variation of 53 L for a 56-m-high Douglas-fir tree [Pseudotsuga menziesii (Mirb.) Franco]; Zweifel and Häsler (2001) reported a maximum daily water depletion of 5 L for a mature Norway spruce tree; Waring and Running (1978) evaluated the daily depletion of the xylem water storage of a mature Douglas-fir stand to be 1.7 mm. On a sunny day with high transpiration demand, the maximum amount of water that can be taken from the plant water storage is determined by the difference between the pre-dawn value of Hrx and Hcrit.

Considering the predicted actual transpiration, the nonzero plant hydraulic capacitance approach allows the system to satisfy a higher transpiration demand. Estimated actual transpiration amounts during the simulated season (Fig. 3) increased from 301 mm, assuming zero capacitance, to 330 mm assuming a capacitance of 9 × 10−6. The change represents 8% of the cumulative potential transpiration. The average daily contribution to the transpiration due to the elastic plant water storage is 0.2 mm d−1.

Assuming a nonzero hydraulic capacitance significantly reduces the simulated soil water redistribution via roots (Fig. 4). The zero-capacitance simulations predict regular nightly redistribution of the soil water demonstrated by negative root water uptake intensities (Fig. 4c). If a nonzero plant capacitance is assumed, the water redistribution via roots is predicted only occasionally during major rainfall events (Fig. 4d).

As expected, the transpiration stream model is quite sensitive to the magnitude of the whole-plant hydraulic capacitance. Figure 5 shows differences in predictions of root water uptake and actual transpiration rates obtained with different values of C, ranging from 0 to 2 × 10−5.

Isohydric vs. Anisohydric Control

To a first approximation, an anisohydric plant strategy can be simulated by setting the critical value of the xylem water potential, Hcrit, equal to Hcav, which marks conditions under which the plant becomes susceptible to xylem embolism. Figure 6 compares the results of simulations performed assuming isohydric and anisohydric strategies of two different transpiring plants—Norway spruce vs. a hypothetical anisohydric tree under the same environmental conditions. The value of Hcrit = Hcav was set equal to −250 m for the anisohydric species, while Hcrit was fixed at −150 m for the isohydric species. In our example, the Hrx value, resulting from the simulation of an anisohydric strategy, did not reach the value of Hcav, thus the effect of increased hydraulic capacitance due to cavitation need not be considered. That, however, would not be true under the conditions of more pronounced drought stress.

To simulate the effects of xylem embolism under severe long-lasting drought stress, the hydraulic capacitance would need to be varied based on not only the Hrx value but also the history of Hrx development. Such a description would require a quantitative evaluation of the related plant physiological processes (e.g., the cavitation-induced loss of hydraulically active tissues, decrease of xylem axial conductance, and the subsequent process of regeneration).

A previously developed one-dimensional numerical model of soil water flow, involving macroscopic vertically distributed plant root water uptake, was extended to allow simulations of transient plant water storage effects. The resulting model predicted reasonably well the evening recharge of plant water storage by root water uptake as well as morning discharge to the atmosphere. The proposed approach can be used to simulate plant responses to water stress caused by high midday potential transpiration demand and/or by soil drought.

The comparison of transpiration stream fluxes with sap flow observations revealed improved model performance when the conventional assumption of quasi-steady-state flow along the transpiration stream was replaced by the more realistic transient plant water storage assumption.

The suggested transpiration stream model is based on an explicit evaluation of the root xylem water potential as a function of the actual transpiration rate and therefore can be easily incorporated into existing models of soil water flow to facilitate more realistic simulations of plant root water uptake and transpiration.

The research was funded by the Czech Science Foundation, Project no. 16-05665S. Sap flow and meteorological data from the Liz experimental site were available through the courtesy of the Institute of Hydrodynamics of the Czech Academy of Sciences.