Introduction

An approach to model plant water relations has developed during the last two decades based on two main concepts: the cohesion-tension theory (Tyree and Zimmermann 2002) and the electrical circuit analogy applied for simulating water transport in plants using resistances, capacitances, water potentials and flow (Cruiziat et al. 2002). By applying these two concepts, the description of tree hydraulic architecture has made a strong improvement towards a more realistic and comprehensive vision of tree water relationships (Cruiziat et al. 2002). Only recently, this approach was further improved by substituting the electrical circuit analogy by a porous media description for sap flow in wood that is based on Darcy’s law and includes a water capacity term to better account for the dynamic behaviour of the hydraulic storage in trees in a way that mass conservation is considered, e.g. by Bohrer et al. (2005), see also Arbogast et al. (1993), Früh and Kurth (1999), Kumagai (2001) and Aumann and Ford (2002). This avoids calculations of negative water contents and unlimited water withdrawal from the tree as occurring with applications of the electrical circuit model (Chuang et al. 2006). By using a three-dimensional (3D) representation of the tree structure and the physically based representation of tree hydro-dynamics, questions can now be adressed of how the variability in crown architecture (i.e. inter-specific, age or ecosystem dependent structure) would lead to different transpiration responses (Bohrer et al. 2005). In particular by considering the stomatal response to the corresponding branch water potential, we can then simulate the stomatal closure, the corresponding change in leaf water conductance and the resulting actual transpiration rate. Moreover, considering a loss of conductivity with decreasing water potential allows the representation of vulnerability to cavitation at critical highly negative water potential values (Hölttä et al. 2005; Cochard et al. 1999).

Since water transport from the soil through the plant into the atmosphere takes place in a soil-plant-air continuum that is interconnected by a continuous film of water, modelling of plant water transport has to consider both the water exchange at the leaf-air interface and the water flux at the soil-root interface. Therefore, also root water uptake including root hydraulic architecture has to be described to integrate the tree hydro-dynamic model into ecosystem models (Doussan et al. 2003), a step still to be taken to replace current representations of plant water transport, e.g. big leaf or resistor-capacitor approaches (Bohrer et al. 2005) that are often coupled to one-dimensional (1D) effective root water uptake models (Cowan 1965; Gardner 1960; Nimah and Hanks 1973; Feddes et al. 1978; Campbell 1985). New root water uptake models are available that explicitly describe root architecture and related soil-plant processes in three dimensions by explicitly considering the 3D distribution of the uptake (Clausnitzer and Hopmans 1994; Somma et al. 1998; Vrugt et al. 2001) and also considering water flow in the root system (Doussan et al. 2006; Javaux et al. 2008). They can be based on root growth models that allow the integration of a great diversity of environmental conditions and their impact on root system development (Dunbabin et al. 2002).

Despite the more realistic assumptions for predicting soil-root interactions, a disadvantage of the 3D models is their high level of complexity and consequently their high computational demand. Therefore, to describe and simulate transpiration at the stand level a less complex but still realistic approach is needed. By extending the hydro-dynamics model for aboveground parts of trees of Bohrer et al. (2005), we developed a simple root water uptake and transpiration model that couples a 1D soil water flow model with a 1D plant xylem water flow model. Notwithstanding its 1D character the xylem water flow model takes 3D plant architecture into account and can consider different root properties for each root node and different root types. Moreover, by specifying different properties of the soil directly surrounding each root, to a certain extent also horizontal variability might be represented.

Model development

In our model both the water flow within the plant and the water flow in the soil are described by applying the 1D porous medium or Richards equation (Richards 1931). This is based on the assumption that similar to the soil matrix, also the plant xylem can be conceived as a porous medium. The xylem is seen as a bundle of small parallel capillary pores, which are filled with water and air. A capillary can conduct water only if it is completely saturated. With increasing suction head more and more pores cavitate and the ability of the xylem to conduct water decreases. To simulate the xylem water flow and the water uptake from the soil by the roots we need a representation of the plant architecture as model input to define the flow domain and the interfaces between the plant and its environment consisting of the soil and the atmosphere.

Description of plant architecture

The plant architecture is described as a combination of stem, branches, leaves, gross and fine roots including branching roots of several orders (Fig. 1a). The area of leaves on an outer branch is prescribed by a relative leaf area distribution. Stem, branches and roots are divided into slices that are approximated by finite cylinders of different lengths and diameters (Fig. 1a). Each of these cylinders consists of an inner cylinder representing the heartwood, an outer hollow cylinder for the sapwood or xylem and further outer hollow cylinders that include the cambium and the phloem (Fig. 1c). It is assumed that the diameters of these cylinders are maximal diameters in the sense that they do not increase anymore due to an increase in water content of the plant and that a loss of plant water can lead to smaller diameters (Fig. 1d). The flow domain of the xylem is further abstracted and represented by a 3D graph consisting of 1D elements each standing for an inner or hollow cylinder of xylem tissue (Fig. 1b).

Fig. 1
figure 1

Schematic illustration of the tree model: a Tree architecture and considered water flows: soil water flow, rootwater uptake, xylem water flow and transpiration. The volume V e of the xylem element e is given by its length l e and its maximal xylem area s x, max . b Tree graph representing the solution domain of the Richards equation for xylem water flow. c Descriptions of above- and below-ground cylindrical tree elements. d Water flow through a xylem cylinder of length l e and changing diameter between maximal diameter d e, max and minimal diameter d e, min

Water flow within the plant

The basic state variable to describe water flow in plants is the hydraulic or water potential of the xylem which quantifies the energy status of the liquid phase. The hydraulic potential can be defined as the sum of the xylem potential or xylem hydrostatic potential and the gravitational potential assuming that changes of the osmotic potential are not important for long distance water flow in the xylem (Früh and Kurth 1999). Hence it is assumed that the only component of the xylem potential comes from the energy to move water isothermally and reversibly into or out of the porous matrix represented by the xylem tissue. Because this potential is determined by the porous matrix of the xylem in the following it is called xylem matric potential in analogy to the soil matric potential (Passioura 1980). If the xylem as the water conducting tissue of the plant is conceived as a porous medium (Siau 1984), water flow can be described by Darcy’s law, the relation between the water potential gradient of the xylem and the mass flow of water (Früh and Kurth 1999; Chuang et al. 2006). Following further Früh and Kurth (1999) we assume homogeneity of the xylem hydraulic characteristics, the capacitance or water retention and the axial hydraulic conductivity for each subunit of the plant xylem as represented by the xylem cylinders in Fig. 1a. By further assuming that the xylem volumetric water contents are defined in relation to fixed volumes of the xylem cylinders, we can apply the principle of mass conservation as expressed mathematically by the 1D continuity equation (Früh and Kurth 1999). For the sake of simplicity, in a first step, we neglect water exchange between xylem, phloem, cambium and heartwood and confine the water flow model mainly to the longitudinal direction of the stem and branch axes. By adding a sink-source term to represent transpiration and in our case also root water uptake, we get the following 1D Richards equation as outlined in detail by Früh and Kurth (1999) and Chuang et al. (2006). This equation can then be solved on the domain given by the graph, which represents the plant xylem porous medium by straight lines not only for stem and branches, but also below-ground for roots (Fig. 1b) extending in this way the approach of Bohrer et al. (2005):

$$ \frac{\partial \theta_x(\psi_x)}{\partial t} = \frac{\partial}{\partial z} \left[k_x(\psi_x) \left(\frac{\partial \psi_x}{\partial z} + \cos\alpha_x\right)\right] - S_x $$
(1)

where for all of the tree elements represented by straight lines θ x denotes the volumetric water content (m3 m − 3) of the xylem as function of the xylem matric potential ψ x (mm). The xylem matric potential is given on a weight basis, i.e. it is expressed as hydraulic head (mm). t (s) denotes the time, z (mm) the axial length of the element and k x (ψ x ) represents the xylem hydraulic conductivity (mm s − 1) defined as function of the xylem matric potential using the maximal xylem cross sectional area as reference surface of the water flux. The vertical position is given by the height above (positive upward) or the depth below the soil surface (negative downward). For a branch or root element α x is the zenith angle (–) and S x (s − 1) the sink or source term, which refers to roots for root water uptake and outer branches for transpiration. In a further model development also the radial exchange of water between the xylem and the phloem, the cambium and the heartwood of the tree might be added to this term.

Hydraulic characteristics of the xylem

Due to the form of the mass conservation principle involved in the formulation of the Richards’ equation, xylem water contents are defined in reference to fixed xylem volumes, i.e. in our case to the maximal xylem volumes. Therefore, if one wants to apply Richards equation and also consider xylem desaturation which keeps the xylem vessels still saturated but leads to significantly smaller xylem volumes, then this desaturation has yet to be described in relation to the maximal xylem volumes. In this case the usual form of soil water retention curves, that are nearly constant between zero matric potential and air entry value, cannot be applied to describe the xylem water–xylem matric potential relationship on a maximal sapwood volume basis. Hence to account for changes in xylem volume, in view of measurements for European beech (Fagus silvatica L.) branches by Oertli (1993), see Fig. 3a, we assume a linearly decreasing xylem water retention curve between zero matric potential at maximal saturation and air entry value. Thus we extend the xylem water retention curve of Chuang et al. (2006) describing for each tree element e of stem, branch or root the dependency of the volumetric xylem water content on the xylem matric potential at a certain height or depth z (mm) by

$$ \theta_x(\psi_x) = \begin{cases} \theta_{x,a} \left(\displaystyle\frac{\psi_x}{a}\right)^{-\lambda},\ \ \ \mbox{if}\ \psi_x <a &\cr \cr (\epsilon_x-\theta_{x,a})\left(\displaystyle\frac{a-\psi_x}{a}\right)+\ \theta_{x,a},& \text{else} \end{cases} $$
(2)

where the maximal volumetric water content of the xylem is assumed to be equal to the xylem sapwood porosity ϵ x (m3 m − 3) given for each tree element by the fraction of the maximal volume of xylem water at saturation to the maximal total volume of the xylem sapwood in the tree element or due to the assumed cylindrical geometry of the xylem elements by the fraction of the maximal water-filled xylem sapwood area s xw, max ,e (m2) and the maximal xylem sapwood area s x, max ,e (m2) of the xylem element e. In a similar way we define the xylem water content at air entry θ x,a (m3 m − 3) by use of the corresponding water-filled cross sectional xylem sapwood area s xw,a (m2) at the xylem matric potential when cavitation begins to occur:

$$ \theta_{x,a}(z) = \frac{s_{xw,a}(z)}{s_{x,\max,e}(z)} $$
(3)

Assuming that tension changes the volume of water conducting sapwood only (Perämäki et al. 2001) and neglecting corresponding changes in phloem water content, which are considered to be comparatively small, then, in combination with the xylem water retention curve, the xylem water content θ x,a at air entry can be derived using the relation between xylem matric potential and relative volume change of the xylem element (Steppe and Lemeur 2007) expressed by

$$ \frac{1}{V_e}\frac{dV_e}{d\psi_{x}} = \frac{d\theta_{x}}{d\psi_{x}} = \frac{1}{E} $$
(4)

where E (mm) denotes the elastic modulus of the xylem and V e the volume of the xylem element. By evaluating this relation, which is valid in the range of xylem water contents between ϵ x and θ x,a when the vessels are saturated, we get

$$ \theta_{x,a}=\epsilon_x+\frac{a}{E} $$
(5)

and hence a value for the xylem water content at air entry can be determined for a given value of the elastic modulus. Moreover, inserting Eq. 5 into Eq. 2 we get

$$ \theta_{x}=\epsilon_x+\frac{\psi_x}{E}\ \ \ \ \ \ \mbox{for $a\le \psi_x \le 0$} $$
(6)

and from this the actual xylem cross-sectional area s x,e (m2):

$$ \begin{array}{rll} s_{x,e}&=&s_{xw,e} + (s_{x,\max,e} - s_{xw,\max,e}) \\ &=& s_{x,\max,e} \left(1+\frac{\psi_x}{E}\right)\ \ \ \ \ \ \mbox{for $a\le \psi_x \le 0$.} \end{array} $$
(7)

We then can calculate the actual xylem diameter exploiting the assumed cylindrical volume geometry of the xylem element e (Fig. 1d).

The elastic modulus E (mm), the porosity ϵ x (m3 m − 3), the air entry value a (mm) and the exponent λ (–) that determine the xylem water retention curve, are model input parameters and assumed to be fixed values for the whole tree.

From the xylem water retention curve we can derive the xylem hydraulic conductivity k x based on the law of Hagen and Poiseuille for the mass flow rate of water in a cylindrical pipe. According to Burdine (1953) by considering capillary bundles of such pipes, the relative conductivity K(θ x ) (–) of the volumetric water flux in the capillary bundle as defined by Eq. 11 can be calculated by the following integral equation (Hoffmann-Riem et al. 1999):

$$ K(\theta_x)= (\theta_x/\epsilon_x)^p\ \ \frac{\displaystyle \int_0^{\theta_x/\epsilon_x}\ \psi_x(\theta)^{-2}\ d\theta}{\displaystyle \int_0^1\ \psi_x(\theta)^{-2}\ d\theta}, $$
(8)

Here also other functional forms such as Mualem’s equation (Mualem 1976) might be applicable. To account for the tortuosity of soil pores usually p = 2 is assumed, but for the sake of simplicity and lack of appropriate data to estimate an adequate value of p we assume p = 0, although tortuosity of water films in xylem vessels might occur during embolism. The evaluation of Burdine’s integral equation for θ x  ≤ θ x,a and the application of the Hagen-Poiseuille law to consider the effect of an increase ΔR of the effective mean hydraulic radius of the capillary bundle for θ x  > θ x,a leads to the following function of the xylem water content to represent the xylem hydraulic conductivity:

$$ K(\theta_x) = \begin{cases} \left(\displaystyle\frac{\theta_x}{\epsilon_x}\right)^{\frac{2}{\lambda}+1}, \quad \mbox{if}\ \theta_x \le \theta_{x,a} \\ \left(\displaystyle\frac{\theta_{x,a}}{\epsilon_x}\right)^{\frac{2}{\lambda}+1} +\left[1-\left(\displaystyle\frac{\theta_{x,a}}{\epsilon_x}\right)^{\frac{2}{\lambda}+1}\right]\\ \qquad\qquad \times\left(\displaystyle\frac{\theta_x-\theta_{x,a}}{\epsilon_x-\theta_{x,a}}\right)^2,\quad\text{else} \end{cases} $$
(9)

For the saturated case, i.e. for θ x  > θ x,a we consider the increase of the mean effective hydraulic radius ΔR of the capillary bundle with increasing radii of the water-filled tubes resulting in xylem volume increase, and get for the xylem radii r = r(θ x ), r max  = r(ϵ x ), r a  = r(θ x,a) (mm) at the respective xylem water contents θ x , ϵ x and θ x,a and the mean value θ m  = (ϵ x  − θ x,a)/2:

$$\begin{array}{rll} \Delta R^2&=&\frac{1}{\theta_{m}} \int_{\theta_{x,a}}^{\theta_x} \frac{r(\theta)^2-r_a^2}{r_{\max}^2-r_a^2}\ d\theta \\ &=&\frac{2}{\epsilon_x-\theta_{x,a}}\int_{\theta_{x,a}}^{\theta_x} \frac{\theta-\theta_a}{\epsilon_x-\theta_a}\ d\theta\\ &=&\left(\frac{\theta_x-\theta_{x,a}}{\epsilon_x-\theta_{x,a}}\right)^2 \end{array}$$
(10)

To account for the change of cross-sections between tree elements, for each tree element at depth or height z the hydraulic conductivity is calculated in relation to the water-filled cross sectional area of the xylem sapwood (Bohrer et al. 2005), normalised by the maximal water-filled cross sectional xylem sapwood area of the stem element at stem-base, i.e. at height zero s xw, max (0) (m2):

$$ k_x(\theta_x(z)) = \frac{s_{xw,\max,e}(z)}{s_{xw,\max}(0)}\ k_{\max}\ K(\theta_x)\ \ \ , $$
(11)

where k max (mm s − 1) denotes the maximal hydraulic conductivity of the xylem. k max is an additional input parameter, which for the sake of simplicity is assumed to have fixed values for either shoot and root. For young European beech trees we apply in the following different values for k max (see Table 1).

Table 1 Input parameters

Water flow in the soil

To simulate vertical soil water flow we apply again the 1D Richards equation, an approach frequently used to model soil-plant-atmosphere systems (Priesack and Gayler 2009; Simunek et al. 2008; Kroes and van Dam 2003):

$$ \frac{\partial \theta_s(\psi_s)}{\partial t} = \frac{\partial }{\partial z} \left[ k_s(\psi_s) \left(\frac{\partial \psi_s}{\partial z} - 1 \right)\right] - S_w $$
(12)

where θ s is the volumetric soil water content (m3 m − 3) as a function of the soil matric potential ψ s (mm), t (s) denotes time and z (mm) soil depth (here positive downward). k s (ψ s ) (mm s − 1) is the soil hydraulic conductivity which is given as a function of the soil matric potential. Both soil hydraulic property functions θ s and k s are described by parameterisations according to van Genuchten (1980). The last term S w (s − 1) represents the sink term of root water uptake (per unit soil depth). Typical boundary conditions of the flow equation are a mixed boundary condition at the soil surface to represent atmospheric conditions leading either to infiltration, evaporation or water ponding and a free drainage condition at the bottom boundary (Simunek et al. 2008; Kroes and van Dam 2003).

Soil water uptake and transpiration

Water flow between the soil-root interface and the root xylem is defined by the 1D soil-root volumetric flux j r,e (mm3 s − 1) which describes the volume of water exchanged between a root element and a soil layer per unit of time

$$ j_{r,e} = k_r\ s_{r.e}(z)\ [\psi_s(z) - \psi_x(z)] $$
(13)

where k r (s − 1) is the radial conductivity of the root through the root cell structure between root xylem and soil and s r,e(z) (mm2) the outer surface of the root element e within the soil layer at depth z (mm). ψ s (z) and ψ x (z) denote the matric potentials in soil at the root surface and in the xylem (mm) at soil depth z. Instead of applying only the radial conductivity k r , we alternatively consider the rhizosphere radial conductivity k rs (s − 1), defined by

$$ k_{rs}(z) = \sqrt{k_r\ k_s(\hat{\psi}_s(z))/l_{rs}} $$
(14)

the geometric mean of the radial conductivity k r of the root and the soil hydraulic conductivity k s at the soil matric potential near the roots divided by the radial thickness l rs (mm) of the rhizosphere soil around the roots assumed to be 2.0 mm (Hinsinger et al. 2005). To account for the impact of the strong root-soil interaction in the rhizosphere on the water status of the soil directly surrounding the roots, we take instead of the soil water matric potential ψ s a weighted mean \(\hat{\psi}_s\) between xylem and soil matric potentials

$$ \hat{\psi}_s(z)=\frac{\gamma_1\ \psi_x(z) + \gamma_2\ \psi_s(z)}{\gamma_1 + \gamma_2} $$
(15)

where γ 1 (–) and γ 2 (–) are weighting constants.

The calculated flux j r,e into or out of the cylindrical root element e is then used to calculate the below ground water uptake sink or source terms S x and S w of the Richards equations of xylem and soil water flow. The below ground xylem sink or source term S x (mm3 mm − 3 s − 1) is obtained by relating the fluxes to the corresponding xylem cylinder volumes V e  = s x, max ,e(z) l e of the root elements of maximal xylem sapwood area s x, max ,e (mm2) and element length l e (mm):

$$ S_x(z) = \sum\limits_e \frac{j_{r,e}(z)}{s_{x,\max,e}(z) l_e} $$
(16)

The water fluxes between soil and root cylinders located in a certain numerical soil layer i are related to the corresponding soil volume ΔV i given by the thickness of the soil layer Δz i and the unit area of a square meter and are summed up to give the sink or source term S w,i (mm3 mm − 3 s − 1) of the soil layer i. If the axis of a root cylinder element e intersects with more than one of the numerical soil layers, the water flux is divided between these layers proportionally to the fractions f e,i (–) of the cylinder axis intersections with the particular layers:

$$ S_{w,i} = \sum\limits_e \frac{j_{r,e}(z) f_{e,i}}{\Delta V_i} $$
(17)

Similar to water uptake, also transpiration, the loss of water to the atmosphere, is described by a sink term. At the outer branches, where the leaves are located, the corresponding water flux is prescribed by the atmospheric evaporative demand and by the stomatal hydraulic conductivity, which similar to the xylem hydraulic conductivity is assumed to depend on the xylem matric potential of the branch element (Bohrer et al. 2005). For simplicity the atmospheric water demand is prescribed by a constant flux condition or calculated from the potential evapotranspiration rate (mm s − 1) estimated by the FAO Penman-Monteith method (Allen 2000) and is distributed according to the leaf area. Therefore, in contrast to the water exchange with the soil, the exchange with the atmosphere is prescribed and no direct feedback of transpiration to the micro-climate at the leaf level is considered.

The prescribed flux condition or the potential transpiration T pot (mm s − 1) of the total tree is then distributed to single branch elements according to the relative leaf area distribution LAD e (z) (m2 m − 2) at each element e, given by the leaf area of all leaves that are connected to the branch element divided by the total leaf area of the tree. Hence the potential transpiration T pot,lf,e (mm s − 1) at the branch element e can be given by:

$$ T_{{\rm pot},lf,e} = T_{\rm pot}\ LAD_e(z) $$
(18)

The stomatal hydraulic conductivity k x,stom (mm s − 1) is modelled in dependence of the xylem matric potential in the leaf ψ x,lf (mm) by

$$ k_{x,\rm stom}(\psi_{x,lf}) = k_{\max,\rm stom}\ \exp\Bigl(-\Bigl[\frac{-\psi_{x,lf}}{b}\Bigr]^{c}\Bigr) $$
(19)

following Bohrer et al. (2005) with a stomatal maximal hydraulic conductivity k max ,stom = 1.0 mm s − 1 and stomatal hydraulic conductivity parameters b (mm) and c (–) and by additionally assuming ψ x,lf = ψ x (z ep ) at the end point z ep of the branch element where the leaf is anchored. The actual prescribed transpiration rate T act,lf,e (mm s − 1) is then calculated from the corresponding potential transpiration rate T pot,lf,e (mm s − 1) at the branch element by

$$ T_{{\rm act},lf,e} \!=\! \min\left\{T_{{\rm pot},lf,e}\ \displaystyle\frac{k_{x,\rm stom}(\psi_{x,lf})}{k_{\max,\rm stom}};\ q_{\max,lf}\right\} $$
(20)

for which the maximal evaporation rate from the leaf q max ,lf (mm s − 1) at the actual leaf xylem potential ψ x,lf is given by

$$ q_{\max,lf} = k_x(\psi_{x,lf})\ \frac{\psi_{x,lf}-\psi_{0,lf}}{\Delta z_{ep}} $$
(21)

where \(\psi_{0,lf}= - 3.0\ \cdot 10^5\) mm denotes a potential that corresponds to a minimal potential in the leaf, which can be reached under dry air conditions. Δz ep is the length of the element belonging to the end point where the leaf is anchored.

Finally, the above ground xylem sink term S x,e (mm3 mm − 3 s − 1) which represents transpiration of the outer branch element e is calculated by relating the volumetric water flux out of the xylem given by T act,lf,e s x, max ,e (mm3 s − 1) to the xylem volume V e (mm3) of the branch element:

$$ S_{x,e} = \left.T_{{\rm act},lf,e} s_{x,\max,e}\right/V_e = \left.T_{{\rm act},lf,e}\right/l_e $$
(22)

Numerical implementation

The soil and tree water flow equations are coupled via the below ground water uptake sink terms and each solved by iterations using a standard Newton–Picard fixed point iteration scheme (Celia et al. 1990; Priesack 2006) until tree and soil water contents and potentials converge below a threshold value (Huang et al. 1996). Although both flow equations may be coupled by a further fix-point iteration via the water uptake term to account for the non-linearity of the systems, similar to Tseng et al. (1995) we use a single-pass solution scheme without iteration between the two subsystems, since for the considered simulation scenarios the additional iteration does not lead to significantly different results. The solution procedures of both flow equations 1 and 12 are based on a finite element discretisation following the approach of the HYDRUS-1D model (Priesack 2006; Simunek et al. 2008). Since the solution method of the Richards equation on a graph as represented in Fig. 1b has not been documented in more detail, we describe the applied finite element method in Appendix A.

Statistical criteria

Three statistical criteria are used to assess the adequateness of simulations. The root mean square error RMSE as defined by Mayer and Butler (1993) is used to assess over the whole period of the lysimeter experiment the deviation between predicted, P i , and observed values, O i , proportioned against the mean observed value \(\bar{O}\):

$$ RMSE = \frac{1}{\bar{O}} \sqrt{\frac{\sum_{i=1}^n \left( O_i-P_i\right)^2}{n}} $$
(23)

where n is the number of observations during the experiment. The lower the value of RMSE, the better is the agreement between simulation and observation. The Modelling efficiency ME (Willmott 1982) is used to assess the model adequacy in describing the observed variability.

$$ ME = 1 - \frac{\sum_{j=1}^n \left(O_j-P_j\right)^2}{\sum_{j=1}^n \left(O_j-\bar{O}\right)^2+\sum_{j=1}^n\left(P_j-\bar{P}\right)^2} $$
(24)

where P j and O j are the predicted and observed values and \(\bar{O}\) and \(\bar{P}\) are the corresponding mean values. A value of ME < 0 indicates that \(P_j = \bar{O}\) would be a better estimation for observed variability than the simulation results. If simulated values are completely in accordance with observed values, ME equals 1.

Finally, a necessary criterion to evaluate the accuracy of a numerical scheme is given by the global mass balance error \(MB = \sum_{m=1}^{\rm end} MB(t^m)\) at the end of the simulation, where at each time step m we have for the water balance of the soil-plant system:

$$ MB(t^m) = 100 \cdot \left|1.0 - \frac{\displaystyle \sum\limits_{e=1}^n V_e\ (\theta_{x,e}^m-\theta_{x,e}^0) + \sum\limits_{i=1}^k \Delta V_i\ (\theta_{s,i}^m-\theta_{s,i}^0)}{\displaystyle \sum\limits_{j=1}^m \sum\limits_{e=1}^n V_e\ S_{x,e}^j \Delta t^j + \sum\limits_{j=1}^m (Q_0^j-Q_k^j+ \sum\limits_{i=1}^k \Delta V_i\ S_{w,i}^j) \Delta t^j}\right| $$
(25)

and where \(\theta_l^m\) denotes the volumetric water content at time step m of the tree at element l = e or of the soil at node l = i, ΔV i is the volume of the soil layer i and \(Q_0^j\) and \(Q_k^j\) denote the volumetric water fluxes (mm3 s − 1) across the boundaries of the soil profile.

Simulation scenarios and parameterisation

The new model of water flow in soil-plant systems is applied to three different scenarios mainly to illustrate that important features of water flow in plants and soils can be simulated. In a first scenario, the drying soil scenario, we simulate soil water uptake and transpiration under conditions of a prescribed constant transpiration rate without refill of soil water. By this scenario the effects of stomatal conductivity, xylem conductivity, radial root conductivity and soil conductivity on water uptake and transpiration are considered. In particular the onset of cavitation is examined. In a second scenario, the hydraulic lift scenario, we demonstrate the basic ability of the model to simulate transport of water through the roots from wetter to dryer soil regions. In a third scenario, the lysimeter scenario, we test if the model can simulate daily and seasonal water balances based on measured input parameters. In this case for the xylem water flow model only parameters of the tree architecture are changed: the stem diameter and the distributions of leaf and root area index based on measurements from the lysimeter experiment.

In the drying soil scenario, we simulate soil water uptake and transpiration also for comparison with a 3D root water uptake modelling approach. In the scenario, the water is taken up by the roots from a homogeneous soil column of 0.4 m depth and of 0.1 m by 0.1 m soil surface, which is assumed to be initially in hydrostatic equilibrium with an aquifer located 3.0 m below the soil surface (Javaux et al. 2008). The soil boundary conditions are a no flux condition at the bottom of the soil profile and no rainfall or evaporation at the soil surface. For reasons of comparison with the numerical experiment of Javaux et al. (2008) we use the same parameters to study the effect of different soil types and to analyse the sensitivity of different soil and xylem hydraulic conductivity values (Table 1). For the reference case input parameter values are taken from the second column of Table 1, and for the sensitivity analysis only the conductivity parameters, k sat for the soil, k max for the xylem and k r for root water uptake, were changed as indicated in the third column. This perturbation of parameter values corresponds to a realistic variation in observed conductivity values for soil or xylem (Javaux et al. 2008; Doussan et al. 1998). In both cases for the above-ground tree a constant potential transpiration rate of 1.0 mm d − 1 was assumed and prescribed along the tree at the ends of stem or branch elements according to the assumed leaf area distribution. Related to the 0.01 m2 soil column surface this corresponds to a volumetric water flux rate of 1.0 ·104 mm3 d − 1 from the soil to the atmosphere.

In the hydraulic lift scenario, instead of the no flux of the first scenario we consider a Dirichlet boundary condition which prescribes saturated conditions at the lower boundary of the soil column now assumed to be 1.5 m deep. Furthermore, we now prescribe a daily sinusoidal cycle of the potential transpiration according to (Childs and Hanks 1975; Priesack 2006), but keep the daily rate at 1.0 mm d − 1. All other conditions and parameters are kept the same as in the loam reference case of the first scenario.

In the lysimeter scenario, the soil-plant systems are represented by 2–3 years old European beech trees growing in cylindrical lysimeters of 2.0 m depth and 1.0 m2 surface area. Here the soil boundary conditions are a seepage flow condition at the bottom of the soil profile and an atmospheric flow condition at the top soil prescribing either infiltration, evaporation or water ponding. The third scenario is based on a lysimeter study about effects of elevated ozone and below-ground pathogen infection on juvenile European beech (Winkler et al. 2009). From this study we take hourly and daily water balance data of one of the control lysimeters, i.e. lysimeter S3, that has seen neither an ozone nor a pathogen treatment. The input data that characterize the hydraulic properties of the lysimeter soil profile are derived from measurements given in Gayler et al. (2009) and summarized by Table 2. The stem diameter d stem is 15.8 mm, the xylem conductivity parameter values are the reference values given in Table 1.

Table 2 Soil input parameters of the third scenario

In all three scenarios we use a simplified plant architecture of a young European beech tree (Fig. 2) obtained using digital photographs and photogrammetric evaluations. In the second scenario we extend a central root to the bottom of the soil column to facilitate water uptake from deeper soil layers and in the third scenario the stem diameter is increased keeping the diameter fractions between stem and branch elements the same, whereas the root system is now adapted to represent the measured root biomass distribution.

Fig. 2
figure 2

Simplified plant architecture of a young European beech tree (two different perspectives). Tree xylem hydraulic potential distribution (between −3,600 mm and −36,000 mm suction head) under still non-limiting soil conditions at day 30 of the reference drying soil scenario simulation

We estimate the maximum water-filled cross sectional xylem sapwood area s x, max of each tree element e from the diameter d e (mm) of the tree element according to Aranda et al. (2005). For young European beech trees we fit the new parameterisation of the xylem water retention curve given by Eq. 2 to the measured retention data of Oertli (1993) for shoots of European beech. This fit results in an elastic modulus value of E = 1.5 ·106 mm and an air entry value of a = − 3.0 ·105 mm with exponent λ = 0.854 for each tree element (Fig. 3). Alternatively to the xylem water retention curve given by Eq. 2 we also apply and fit a van Genuchten parameterisation (van Genuchten 1980) obtaining the parameters α = 2.8 ·10 − 6 mm − 1, n = 2.77, θ sat = ε x and θ res = 0 (Fig. 3).

Fig. 3
figure 3

a Relative xylem water retention curve: measured values taken from shoots of European beech by Oertli (1993), new parameterisation (solid line) according to Eq. 2 and van Genuchten parameterisation (dashed line). b Relative xylem hydraulic conductivity curve: measured values taken from shoots of European beech by Magnani and Borghetti (1995): (closed triangles), Cochard et al. (1999): (open triangles), Lemoine et al. (2002): for a sun branch (empty circles) and a shade branch (filled circles), conductivity curves derived from the new parameterisation (solid line) according to Eq. 9 and derived from the van Genuchten parameterisation (dashed line) according to Mualem (van Genuchten 1980)

Results

Numerical accuracy and convergence

All numerical simulations were conducted on an ordinary personal computer with a single core 1.5 GHz AMD processor. CPU time is given in seconds (Table 3). During each simulation the water mass balances are recorded to compare for every individual time step the input-output balances of the whole tree respectively of the soil profile with the corresponding internal mass change in the tree or soil. For the simulation of the soil drying scenario using the reference parameters the global mass balance error of the soil-plant system is MB = 1.3 ·10 − 3%, for the hydraulic lift scenario MB = 3.1 ·10 − 3% and for the lysimeter scenario MB = 3.2 ·10 − 4%. Table 3 summarizes results obtained for the xylem water flow simulation in the drying soil scenario using the reference parameters. By relaxation of the iteration tolerance δ tol it is seen, that already at the criterion of δ tol = 0.01 a rather low global mass balance error is reached without strongly increasing CPU time, iteration and time step numbers. In contrast, the CPU time is considerably increased by increasing the number of tree elements up to 482 or 903 by prescribing a maximal axial element length of only 20 or 10 mm instead of 100 mm as in the reference case. This increase indicates the first order accuracy in space of the applied numerical solver.

Table 3 Summary of numerical results for the xylem water flow model in case of the soil drying scenario with reference parameters

Drying soil scenario

Figures 4 and 5 show the simulated course of transpiration and root water uptake with related plant water storage for different parameterisations of the stomatal hydraulic conductivity curve. Because initially, at simulation start, tree water content is not in equilibrium with soil water and atmospheric conditions, the tree first looses water to the atmosphere before the hydraulic gradients of xylem water potentials get strong enough and force the roots to take up the available water from the soil. Since a constant maximum water flux is prescribed at the tree leaves and no water is added to the soil, the water flow in the soil-plant system quickly reaches an almost steady state as long as the soil provides enough water to fulfil the prescribed flux. After some time, when the soil gets dry, the root water uptake rate decreases due to decreased soil hydraulic conductivities and low soil water availability. This low water availability induces low xylem water potentials and consequently low stomatal hydraulic conductivities such that the transpiration rate is reduced. The transpiration rate is strongly determined by the assumed parameterisation of the stomatal hydraulic conductivity curve which reflects the importance of stomatal control on xylem water flow. For the lower value of the stomatal hydraulic conductivity parameter b the stomatal hydraulic conductivity decreases already at higher xylem potentials and hence protects the tree from water losses earlier during the drying experiment (Fig. 4). If we then assume that the closure of the stomata would lead to zero stomatal hydraulic conductivity the plant could almost keep its water content (Fig. 4b, solid line) without reaching the air entry value. But if we assume a maximal possible reduction of stomatal hydraulic conductivity to only 10% of its maximal value, then the plant internal water is almost completely transpired (Fig. 4b, dotted and dashed lines) and cavitation occurs. Generally, a higher b value leads to a higher stomatal hydraulic conductivity also at lower xylem potentials. Therefore higher water losses of tree water storage take place, which results in critically low tree water contents. In our example in three cases xylem water content falls below the air entry value θ x,a = 0.8, so that cavitation is induced (Fig. 4b). In the cases where the stomatal hydraulic conductivity remains above a non-zero limit, the transpiration gets finally limited when the plant internal water is almost exhausted (Fig. 4a).

Fig. 4
figure 4

Drying soil scenario (loam): Time course of a transpiration rates and b relative xylem water contents θ x /ϵ x of the total tree using stomatal conductivity parameter b = 1.0 · 105 mm with (solid line) and without (dotted line with open circles) reduction to zero conductivity, and b = 3.0 · 105 mm with (dashed dotted line) and without (dashed line with filled circles) reduction to zero conductivity

Fig. 5
figure 5

Drying soil scenario: Conductivity sensitivity simulation (loam): Time course of the root water uptake: reference (solid line), 10-fold xylem conductivity k max (dashed line), 0.1-fold xylem conductivity k max (dotted line), 10-fold radial root conductivity k r (dot-dot-dashed line), 10-fold saturated soil conductivity k sat (dot-dashed line), a using the radial root conductivity k r , b using the radial rhizosphere conductivity k rs with weight coefficients γ 1 = 1 = γ 2. Soil effect simulation: Loam (solid line) compared to clay loam (dashed line) and clay (dotted line), c using the radial root conductivity k r , d using the radial rhizosphere conductivity k rs with weight coefficients γ 1 = 1 = γ 2

Due to the evaporative water flux from the leaves to the atmosphere, the simulated xylem matric potentials show a strong vertical gradient along the tree height (Fig. 2) having lower potentials at the tree crown, where the leaf area is highest and higher potentials towards the roots and within the roots which take up the water from the soil.

The perturbation of conductivity values gives only very small differences between the simulations of the reference case and the cases with one-tenth- or ten-fold maximal xylem conductivity (Fig. 5). This small deviation is similar to results obtained by Javaux et al. (2008) when using the same reference conductivity values for root water uptake simulations (Table 1). Additionally, no difference can be detected if we neglect the second term of the relative conductivity curve in the case θ x  > θ x,a, i.e. if we neglect the impact of an increased mean radius of the water-filled xylem vessel elements. Moreover, for the small trees the exact form of the conductivity curve, which is difficult to evaluate considering the high variability of the measured conductivity data (Fig. 3), may play a less important role, if only the decrease of the curve with decreasing xylem matric potential is appropriately described. Because of the low sensitivity of maximal xylem conductivity variation on the simulated actual transpiration in case of the considered small trees, the relative xylem conductivity curve given by Eq. 9 could be further simplified in the saturated range above the air entry value by omitting the second term of the right hand side without significantly changing our simulation results of xylem water flow dynamics (results not shown). This low sensitivity also means, that the impact of the elastic change on the diameter and hence the conductivity of water-filled tubes of the capillary bundle representing the xylem can be neglected in our model description of xylem water flow in small European beech trees. In this case the low sensitivity of the axial xylem conductivities seems to be caused by their relatively high values compared to the low stomatal and radial root conductivity values.

When compared to the reference case in the cases of increased soil or radial root conductivity the begin of limited water supply for root water uptake is 2–5 days delayed (Fig. 5a). Applying the radial rhizosphere conductivity the delay of the onset of limitation of transpiration is 9–10 days (Fig. 5b). In the sequel, the higher conductivities lead to earlier soil water exhaustion and reduced transpiration rates below 0.2 mm d − 1.

To represent effects of horizontal soil variability we introduce a possibility to account for different soil properties in the direct surroundings of different roots or root elements. This representation is achieved by considering the radial rhizosphere conductivity defined by Eq. 14. The root water uptake rates for the three soil types applying either the constant radial root conductivity k r or the rhizosphere conductivity k rs are shown in Fig. 5c and d. Generally, we observe, that if we apply the radial rhizosphere conductivity k rs instead of the radial root conductivity k r in Eq. 13, the simulated root water extraction decreases earlier and remains slightly higher when the soil gets drier (Fig. 5b and d). For k r the root water uptake limitation appears first for the clay, then for the loam and eventually for the clay loam (Fig. 5c). For k rs the picture changes in case of the loam soil, since the time-span until transpiration gets limited is now shortest for this soil type (Fig. 5d). In both cases the clay loam soil supplies the transpiration water demand by far the best, up to 15 d longer than the clay and loam soils.

The simulated depth distribution of root water uptake is shown by Fig. 6. Under non limiting soil water conditions the root demand is met by available soil water and the depth distribution of the uptake rate almost coincides with the root area index distribution. Only in the case of the increased radial root conductivity available soil water gets already scarce in the region of high root area index, i.e. m2 root surface area m − 2 soil surface area (Fig. 6a). In this region water uptake by roots begins to be limited after 15 d, and is already very low after 30 d (Fig. 6b and c). Then water is increasingly taken up from regions of lower root area index, where soil water is not completely exhausted and still can satisfy almost half of the prescribed demand (Fig. 6d).

Fig. 6
figure 6

Drying soil scenario (loam): Depth distribution of root water uptake and soil water content using different root or soil parameterizations: reference (solid line), 10-fold xylem conductivity (dotted line), 10-fold radial root conductivity (dashed line), and 10-fold soil hydraulic conductivity (dashed-dotted line). The grey line represents the depth distribution of the root area index RAI (m2 root surface area m  − 2 soil surface). Figures ac represent different times t after the start of the drying experiment: a t = 5 d, b t = 15 d, c t = 30 d, and d corresponding volumetric soil water contents

In the cases with increased soil and radial root conductivity, similar to the 3D root water uptake simulations of Javaux et al. (2008), the limitation of root water uptake is significantly delayed by 2–5 days (Fig. 5a). Furthermore, the finding that the clay loam best sustains the evaporative demand of the plant, about 10–15 d longer than the other soil types could be confirmed by our simulations (Fig. 5c and d). But, in contrast to the 3D water uptake simulations under the second collar boundary condition (CBC 2) by Javaux et al. (2008), in our 1D simulations water stress occurred earlier for clay than for loam soil, if the radial root conductivity was used (Fig. 5c). Only if the radial rhizosphere conductivity is applied, the same order of soil types concerning the appearance of water stress is found, i.e. transpiration is reduced first for loam, then for clay, and eventually for clay loam soil (Fig. 5d).

Hydraulic lift scenario

Plant roots do not only take up water but can also release water to dry soil (van Bavel and Baker 1985; Richards and Caldwell 1987). In this way the plant root system may transport water from moist into drier, often upper soil horizons and hydraulic lift occurs. Assuming non polar flow across root membranes we are interested to see if our proposed soil-plant water flow model is in principle capable to simulate such a hydraulic lift. Under the conditions of the second scenario with periodic daily transpiration cycle the simulated water uptake rates show a diurnal variation (Fig. 7). During the day when potential transpiration rates are above zero, water is taken up mainly from the upper soil at 0.0–0.4 m depth (Fig. 7a), where most of the root biomass is located. During the night, when no transpiration occurs, deep roots take up water from the bottom of the soil column, where the soil is nearly water saturated, replenish xylem water and additionally release water to the drier soil near the top of the soil column also at 0.0–0.4 m depth (Fig. 7a). Simulations considering a time span of more than 240 days show that water flow in the soil-plant system gradually reaches a quasi steady state situation, where water is taken up from the lower part of the soil column and redistributed to the upper soil during the night (Fig. 7b), whereas during the day mainly water from the upper soil layer is taken up by the roots to accommodate transpiration demand (Fig. 7c). Moreover, compared to a simulation where we assume no water flux from the roots to the soil, the daily transpiration with hydraulic lift is about 5% higher during the first 100 days.

Fig. 7
figure 7

Hydraulic lift scenario: Root water uptake in the soil profile: a after 120 d within the whole day at night times: 0.2 d (solid line), 0.4 d (long dashed line), 1.0 d (dotted line) and at day times: 0.6 d (dashed dotted line) and 0.8 d (short dashed line); b during night time and c during day time in each case after 1 d (solid line), 60 d (long dashed line), 120 d (short dashed line), 180 d (dashed dotted line) and 240 d (dotted line). d Depth distribution of the root area index RAI (m2 root surface area m  − 2 soil surface)

Lysimeter scenario

A more realistic scenario is given by the lysimeter scenario. Based on the measured daily precipitation, the estimation of daily potential evapotranspiration by the Penman–Monteith dual crop coefficient method (Allen 2000; Loos et al. 2007; Gayler et al. 2009) and the partitioning of potential evaporation and potential transpiration according to Allen (2000), the daily actual evapotranspiration, the daily changes in water storage and the daily percolation rates were simulated and compared to the corresponding measurements (Fig. 8). Results show that the model is able to reproduce the actual evapotranspiration and the changes in water storage quite well, with a model efficiency ME of 0.901 and a root mean square error RMSE of 0.298 for the evapotranspiration, and with the values ME = 0.977 and RMSE = 0.586 for the water storage change. The simulated dates and sizes of the percolation peaks differ more strongly from the measured data giving the values ME = 0.845 and RMSE = 1.812, in particular, the simulated percolation peaks seem to occur always somewhat later than the observed peaks.

Fig. 8
figure 8

Lysimeter scenario: a Daily precipitation, P (mm d − 1, grey bars), leaf area index LAI (m2 leaf surface area m  − 2 soil surface, dotted line); b estimated daily potential evapotranspiration PET (mm d − 1, open circles), together with simulated (solid lines) and measured (grey bars) daily actual evapotranspiration ET (mm d − 1), c daily changes in water storage ΔW (mm d − 1), and d daily percolation rates L (mm d − 1) of lysimeter S3 (Gayler et al. 2009)

Figure 9 shows the simulated diurnal course of transpiration, stem diameter, water flux at the root collar and the rate of change of xylem water content during three consecutive days of the lysimeter experiment. Transpiration and water flux at the root collar strictly follow the prescribed sinusoidal potential evapotranspiration. After the onset of transpiration, the rate of change of tree water content decreases due to the onset of transpiration but soon increases because root water uptake begins to refill plant water capacity, before it decreases again indicating the end of water replenishment of the xylem water by soil water taken up during the night. Correspondingly, the xylem diameter changes with a maximal diameter difference of 0.194 mm between day and night.

Fig. 9
figure 9

Lysimeter scenario: Simulated diurnal course of transpiration (mm d − 1, solid line), of root water uptake (mm d − 1, dotted line), of change rate of xylem water volume per unit soil surface area (mm d − 1) of the total tree (dashed line) and of xylem diameter (mm, grey line) during three days (24–27 July)

Discussion

A new model approach is introduced to describe water flow in the soil-plant system. The model extends previous approaches that were developed to describe hydro-dynamics in above-ground plants (Früh and Kurth 1999; Aumann and Ford 2002; Bohrer et al. 2005) by additionally including the hydraulic root system architecture in combination with a soil water flow model. In this way a complete water flow model for the whole soil-plant system is obtained, which then vice versa extends existing root water uptake modelling approaches (Doussan et al. 1998; Javaux et al. 2008) by including the above-ground water dynamics of the plant. Similar to the soil water flow model, which is based on the continuum approach as described by the Richards’ equation, also the plant water flow is described by this porous media equation following Bohrer et al. (2005). This description is in contrast to more discrete approaches that consider the plant as built up of storage compartments and describe the water flow by water exchange processes between these compartments (Perämäki et al. 2001; Steppe and Lemeur 2007). But, besides the advantages of a continuous mathematical formulation representing mass conservation (Früh and Kurth 1999), due to the similarity between the porous media approaches to describe water flow in the soil and the xylem, it seems also to be conceptually easier to derive a consistent model of hydro-dynamics of soil-plant systems based on the continuum approach.

In a first analysis, the new model was applied to investigate model behaviour and model sensitivity by scenarios with young European beech trees that have already been analysed in two previous studies (Javaux et al. 2008; Gayler et al. 2009). Model tests and numerical results show, that the model is successfully implemented and can efficiently simulate water flow in the soil-plant system based on the a simple representation of plant architecture. Mass balance errors are low and in an acceptable range. But since the global mass balance is a necessary but not a sufficient criterion for having a correct solution (Celia et al. 1990), the model needs further analysis, which may be based on comparison with exact solutions (Lehmann and Ackerer 1998) or as in our case by evaluating simulations of practical examples using experimental data.

In the case of the 3D simulation by Javaux et al. (2008) the effect of the steeper slope of the loam hydraulic properties generates a heterogeneous distribution of soil water potentials in the horizontal direction with drier soil near the roots. This distribution leads to a larger conductivity drop and produces an earlier water stress (Javaux et al. 2008). In the case of our 1D simulation this effect is mimicked by amplifying the direct impact of the lower xylem potential on the soil matric potentials directly near the roots resulting in a lower radial hydraulic conductivity between root axis and soil, which is modelled by Eq. 14 based on the definition of the radial rhizosphere conductivity k rs . If the distribution of radial conductivities of the roots is adapted in this way, the main conclusions about effects of soil type and conductivity values remain, although the root architecture was different in the drying scenario simulations we compared. But amplitudes or local values of considered variables are different due to the 1D simulation of soil water flow. Therefore, a full data comparison cannot be performed, even if locally the water uptake by the root system may be of similar magnitude.

Besides the model sensitivity of the radial root or rhizosphere conductivity, the most significant model sensitivity on simulated soil water uptake and transpiration is due to the stomatal hydraulic conductivity. This model behaviour is mainly caused by eventually rather low conductivity values, reflecting the significant role of both conductivities in controlling tree water flow which has been shown in several experimental studies (Steudle 2000; Lemoine et al. 2002; Sperry et al. 2003). Moreover, the water has to flow through nonxylary pathways in the root across the root cortex and endodermis and in the leaf through the mesophyll. These flow paths have low hydraulic conductivities and can account to a high proportion of the overall root respectively shoot hydraulic resistance (Sperry et al. 2003). In particular, the radial root conductivity may be strongly changed by aquaporin activity (Ehlert et al. 2009), which itself may be controlled by metabolism or be triggered by environmental factors (Steudle 2000). Also changes of leaf hydraulic conductivity may be mediated by plasma membrane aquaporins (Cochard et al. 2007).

Hydraulic lift has been shown to occur in more than 60 species (Espeleta et al. 2004), “but there is no fundamental reason why it should not be more common as long as active root systems are spanning a gradient in soil water potential and that the resistance to water loss from roots is low” (Caldwell et al. 1998). However, for European beech trees it seems to be an open question whether the radial hydraulic root conductivity from inside the roots towards the soil is large enough to permit considerable water flow to dry soil (Rewald 2008). Therefore our simulations of hydraulic lift only show that in principle the model is able to describe this kind of water redistribution from wet to dry soil under almost optimal conditions for hydraulic lift, i.e. assuming water saturation at the bottom of the soil profile with a resulting simulated leaf xylem matric potential ranging from −40,000 to −70,000 mm. Nevertheless, the hydraulic lift could be significant in mitigating plant water stress and need to be accounted for in models describing root water uptake and the onset of plant water stress in a water limited ecosystem (Siqueira et al. 2008).

The lysimeter scenario simulations confirm the applicability of the new model to describe daily water balance dynamics in soil-plant systems. Also the diurnal course of transpiration seems to be realistically described. In particular, the simulated diurnal changes of xylem diameter up to 0.194 mm (Fig. 9) are within the range of xylem or stem diameter changes that were observed for young European beech and Norway spruce (Pinus abies L.) trees (Steppe and Lemeur 2007; Perämäki et al. 2001). This agreement results from the fact that based on the parameterisation of the xylem water retention curve obtained from the experimental curve of Oertli (1993), we get a value for the elastic modulus E which is similar to values found by Steppe and Lemeur (2007) and Perämäki et al. (2001). However, this value might have been underestimated, since we assume that changes of xylem water potential correspond only to changes in the xylem water content neglecting related changes in water contents of phloem, cambium and heartwood.

Since the new plant water flow model focuses on the water flow in the sapwood tissue in the longitudinal direction of the stem and branch axes, the model still needs to be extended to describe water flow in the phloem and to account for the radial exchanges of water between xylem, heartwood, phloem and cambium. In particular the active uploading and unloading of sugars according to Münch’s hypothesis (Münch 1928) needs to be considered, because near the sources of sugars considerable water flow from the xylem to the phloem can occur and vice versa close to sugar sinks significant amounts of water can flow from the phloem to the xylem. Locally this water exchange can strongly influence the xylem water flow dynamics which also can lead to xylem diameter changes (Sevanto et al. 2002).

At the present state of development the model is supposed to be already useful to simulate the impact of differences in plant architecture and soil properties on transpiration and water uptake under various climatic conditions, in particular considering differences of leaf area and root area distribution or of root and shoot branching. The model may be also helpful to describe the impact of different distributions of conductivity values along the tree architecture in particular of larger trees to study their role in protecting the xylem against cavitation, while maintaining transpirational water flow. This description may lead to a better understanding of observed differences in xylem structures and functions between deep and shallow roots, roots, stem and branches (Nadezhdina 2010; North 2004) and of the relevance of nonxylary pathways for water uptake and transpiration (Sperry et al. 2003). The model may also be applied to better describe differences of stand water budgets between mixed forest stands of different species composition.

Conclusion

The newly proposed model of water flow in soil-plant systems presented in this study combines approaches to simulate water flow in the above-ground plant based on the porous media equation with root water uptake and soil water flow models also based on the porous media equation. In this way the model concept exploits the similarity between water flow in the soil and the xylem of plants. It could be shown that the model can simulate water balance dynamics in soil-plant systems with young European beech trees including some main known features of water flow in plants and soil such as the diurnal cycle of plant water content and xylem diameter change, stomatal control of transpiration, cavitation due to water stress, sensitivity of root water uptake to radial root conductivity, hydraulic lift, interaction of hydraulic soil properties, root distribution and plant water availability. Due to the simplified representations of plant architecture and distribution of soil properties that lead to 1D representations of water flow domains, the model is numerically very efficient and may therefore be also applicable to simulate water flow in canopies accounting for spatial and temporal variability of canopy structures based on rather simple representations of individual plant architecture.