Simulation of Water Flow in the Branched Tree Architecture

The model HYDRA, which simulates water flow in the branched tree architecture, is characterized. Empirical studies of the last decades give strong evidence for a close structure-function linkage in the case of tree water flow. Like stomatal regulation, spatial patterns of leaf specific conductivity can be regarded as a strategy counteracting conductivity losses which may arise under drought. Branching-oriented water flow simulation may help to understand how damaging and compensating mechanisms interact within the hydraulic network of trees. Furthermore, a coupling of hydraulic to morphological modelling is a prerequisite if water flow shall be linked to other processes. Basic assumptions of the tree water flow model HYDRA are mass conservation, Darcy's law and the spatial homogeneity of capacitance and axial conductivity. Soil water potential is given as a one-sided border condition. Water flow is driven by transpiration. For unbranched regions these principles are condensed to a nonlinear diffusion equation, which serves as a continuous reference for the discrete method tailored to the specific features of the hydraulic network. The mathematical derivation and model tests indicate that the realization of the basic assumptions is reproducible and sufficiently exact. Moreover, structure and function are coupled in a flexible and computationally efficient manner. Thus, HYDRA may serve as a tool for the comparative study of different tree architectures in terms of hydraulic function.


Introduction
According to network theory, the transmission characteristics of a network, e.g. an electric circuit, are determined by the network structure and the properties of the components (Oppenheim andWillsky 1983, Varju 1977).The branched structures of trees, which determine the flow paths of water, are open networks.Hence, we can understand causally the response dynamics of the tree hydraulic system to external conditions only if we consider the branched tree architecture.
The connection of hydraulic and morphological simulation offers particular chances.Morphological models, like the softwares AMAP (De Reffye and Blaise 1993) and GROGRA (Kurth 1994a), enable the systematic generation of architectural variants, which then could be compared in terms of their hydraulic behaviour.A morphological model depicts the structural matrix jointly underlying different processes of trees (light interception, carbon allocation, water flow, etc.).Thus a coupling of hydraulic to morphological modelling is a prerequisite if tree water flow shall be linked to other processes, e.g.carbon flow (Friih 1995, Kurth 1994b).
A further argument for the branching-oriented approach are spatial patterns of leaf specific conductivity within the tree crown, which are associated with branching laws (Ewers and Zimmermann 1984, Tyree et al. 1983, Zimmermann 1978).Like stomatal regulation, they can be regarded as a strategy counteracting xylem embolism, which may be triggered off by drought (Jones andSutherland 1991, Zimmermann 1983).Enhanced by a positive feedback, this process endangers the functioning of the hydraulic system (Tyree and Ewers 1991).Branching-oriented water flow simulation may help to understand how damaging and compensating mechanisms interact within the hydraulic network of trees.

Basic Assumptions of the Model HYDRA
The fundamental measure to describe water transport in trees is the water potential *F (Pa), which quantifies the free energy per unit volume of an aqueous phase (Slatyer 1967, Mohr andSchopfer 1993).Water flows spontaneously along falling gradients of x ¥.Within plants, *F is a sum of different components: (1 where P is the hydrostatic pressure (Pa); PH20 g z is the gravitational potential (Pa) with PH 2 O as the density of water, g as the gravitational acceleration, z as the height over a reference level; 7C* is the osmotic pressure (Pa); x* is the matric potential (Pa).
Assuming that temporal and spatial changes of the components 71* and T* are without importance for long-distance water transport, eq. ( 1) simplifies to (2) Darcy's law specifies the relation between the ^P-gradient and the mass flow F (kg • s" 1 ): where K is the hydraulic conductivity normalized by length (kg • s~] • MPa~' • m).Darcy's law provides a natural measuring instruction for K. Principally simple is the gravity flow method: by stacking a water column of defined height over a woody specimen, a H'-gradient is given and the flow driven by the gradient is registered (e.g.Zimmermann 1978).The method integrates over the total volume of the sample piece and thus yields a spatially effective value of K.
Hydraulic conductivity is not independent of water potential.Within transpiring plants, hydrostatic pressure normally is negative (Holbrook et al. 1995, Pockman et al. 1995), that is, water is in the state of a superheated fluid (Trevena 1987).
With declining pressure and thus declining potential, the probability of cavitation events increases (Oertli 1971, Tyree andSperry 1989).Vulnerability curves quantify the progressive loss of conductivity with pressure or water content (Cochard 1992, Edwards and Jarvis 1982, Sperry and Tyree 1988,1990, Tyree and Dixon 1986).
The hydraulic capacitance c (kg • dirr 3 • MPa~') represents the ability of tree organs to store or release water with isothermal ^-changes: where 0 is the water content (kg • dm" 3 ), c often is obtained from the slope of dehydration isotherms (Edwards andJarvis 1982, Tyree andYang 1990).As with the measurement of K, the methods yield spatially integrated results referring to the total volumes of sample pieces.
If we synthesize the hydraulic system by a model, then the theoretical synthesis should correspond to an appropriate empirical analysis: the decomposition of the tree to sample pieces, whose material properties are determined, has to ensure well-defined boundary conditions during measurement.In this sense, leaves or leaf blades are useful subunits for the measurement of capacitance.If the network of branching axes is decomposed by cross sections perpendicular to the axial directions, axis segments result, for which spatially effective values of capacitance and axial conductivity can be measured under conditions with geometrically clear boundaries and minimal boundary effects (impairment of tissues, water losses).Hence, it is possible to compose the hydraulic network from axis segments and to compute the system dynamics on the basis of measured hydraulic properties.The underlying assumption is that capacitance and axial conductivity are superimposed in a spatially homogeneous manner.
The subdivision into bark and wood, within the xylem cylinder the differentiation between inner and outer sapwood, parenchymatic and tracheidal/ tracheal tissue, late and early wood, rather indicate a spatial separation of capacitive and conductive functions.However, a quantitative derivation of system dynamics from structure could neither be empirically based nor checked, if the conductor-capacitor-separation is the premise.With ex-ception of the wood-bark relationship, the geometrical relations between capacitive and conductive elements are highly complex (Braun 1970, Wilson andWhite 1986), furthermore, cavitation and refilling of conduits may cause transitions between both elements (Tyree and Yang 1990, Waring et al. 1979, Borghetti et al. 1991).
On the other hand, modelling on the basis of the homogeneity assumption ensures empirical falsificability.The modelled system structure can be backed up with measured data, as it is possible to check system dynamics derived from structure by spatially integrating measurements (the measuring of flow by the heat balance method seems appropriate (Cermäk et al. 1976, Kucera et al. 1977, Schulze et al. 1985, Granier 1987)).
Now, if we further accept mass conservation, we arrive at a trinity of basic assumptions: (1) Mass conservation, (2) Darcy's law, (3) Spatial homogeneity of capacitance and axial conductivity.
Since we restrict ourselves to the axial flow direction, the mathematical expression of mass conservation is the 1-dimensional continuity equation where/is the flux density (kg • s" 1 • nr 2 ).With eq. ( 3) and (4), using K as hydraulic conductance (kg • s" 1 • MPa 1 • nr 1 ), we obtain dt Adding the sink term E(t)l(x) for the transpiration leads to is the transpiration flux density (kg • s" 1 • nr 2 ), £(x) the local leaf area (nr 2 • nr 3 ) with a formal reference to the axis volume which can be dropped in the discrete formulation.Eq. ( 5) is a member of the family of nonlinear diffusion equations.With initial and boundary conditions, an initial boundary value problem (IBVP) results, which defines a unique spatio-temporal dynamics.The soil water potential ^soii is given as boundary condition *¥(x = 0, t) = ^soj] with constant *Fgoii.Since the transpiration is represented as a sink term of eq. ( 5), the boundary condition is one-sided.If we choose as initial state thermodynamic equilibrium, then ^(x, t = 0) = ^Psoji is valid.
The IBVP involving eq. ( 5) simplifies the following aspects: (1) Leaves are not treated as compartments of their own, but leaf areas and capacities are attributed as parameters to the axis segments.( 2) The temporal course of transpiration flux density is a spatially uniform model input.
(3) Stomatal control is not included.( 4) The feedback of root water-uptake to soil water dynamics is not part of the model.

Numerical Method
For the above IBVP no closed-form analytical solution exists, it is to be approximated numerically by the solution of a discrete set of equations instead.Since hydraulic properties change discontinuously at branching nodes, a discrete formulation could be chosen a priori.However, with discrete formulations we may encounter unexpected difficulties.The solution may oscillate with extreme amplitude (as in the model of Tyree 1988), or significant deviations from mass balance may occur (Istok 1989).In general, the relation between a continuous equation and a discrete model is not easy (Ames 1977).To assure a transparent transition from the basic assumptions to the final computation results, we first condensed the assumptions to a continuous model and then selected a discrete method, whose quality thus can be judged with reference to the continuous equation.Therefore we maintain this reference in the limit case of constant cross section which excludes branchings.
Here only the essentials of the numerical method are outlined, more details will be given in a future publication.
For the IBVP involving eq. ( 5), the finite difference method of Douglas and Jones (1963) is chosen.This is a Crank-Nicolson scheme linearized by a predictor-corrector method.The method is unconditionally stable and converges with O[(Ax) 2 + (At) 3 ' 2 ].We were able to show that the simplifications of von Rosenberg (1969) also apply to the specific problem (Friih 1995).
After we had extended the method to variable cross sections and finally to the network of branching axes, we obtained a finite difference method for which the spatial discretization scheme in the vicinity of a branching node is shown in Fig. 1.If we omit for the sake of clarity the temporal discretization, the difference equation for the general branched case is: where L 7 is the leaf area (m 2 ) and A V} the volume (dm 3 ) of the discrete, branched cell shaded in Fig. 1.The cell is bordered by the midpoints Xj-i and xJ + ± with the superior index v to distinguish the grid points located downstream to the branching node.Eq. ( 6) mirrors mass conservation, since the difference between the flow Fj-i into the cell and the partial flows FJ + ± out of it, is balanced by the change of water content and the transpirational water loss.For Fj_± and FJ + ± difference expressions are to be substituted according to Darcy's law (eq.( 3)).Assuming irreversibility of conductivity losses and using vulnerability curves of Cochard (1992), the conductivities of tree segments are actualized, when local pressures decrease.
The parametrization of the discrete IBVP follows a principle of Tyree (1988) which derives the local initial conductivities and capacities from structural data (diameter and volume of axes, leaf area).This makes the adaption of functional information simple when the structural data are exchanged.
For the model tests presented here, Picea abies trees generated by GROGRA were used and parametrized as follows: For the initial conductivities empirical regressions of the form k = a • (ö?/[cm])P were taken from Cochard (1992) with a = 5.370 • 10" 5 kg • s-' • MPa-' • m and (3 = 2.44.The capacitance of shoot axes (0.1 kg/(dm 3 • MPa)) and of leaves (0.04 kg/(m 2 • MPa)) were taken from Tyree (1988) and provisionally treated as constants.We approximated the vulnerability curves of Cochard (1992)  The numerical method requires the solution of two sets of linear equations per time step, where the number J of equations equals the number of spatial grid points.For the solution of the linear sets a Cholesky algorithm is used.Its count of floating point operations normally increases with P (Press et al. 1992).However, we were able to tailor the algorithm to our specific problem, so that the operation count only increases with J.This specific adjustment is crucial for the practicability of network modelling, since we usually have to deal with a J between 10 4 and 10 5 .The branching points and the end points of branching axes determine a metric, the distances of which vary markedly and randomly.No correlation is known between this metric and the spatial pattern of water potential gradients, which may vary by up to 3 orders of magnitude (Tyree and Ewers 1991).Since the locations of these spatial grid points are more or less fixed, the room for the setting of points is markedly constricted in contrast to ordinary discretization problems (e.g.Trompert and Verwer 1990).The spatial discretization in the tree water flow models of Edwards et al. (1986) and Tyree (1988) has no orientation to numerical requirements.But if the dependency of hydraulic function on tree architecture shall be studied, a standardized discretization is necessary which guarantees accuracy independently of structural peculiarities.Hence, a discretization algorithm forms a natural interface between structural and functional modeling.
The multi-step discretization algorithm starts with a so-called base grid which contains no other grid points than branching and terminal nodes.The final product is called a hydraulic map following Tyree (1988).Both, base grid and hydraulic map, are stored as lists of nodes.To each node its own number and the number of the neighbours are assigned.The numbers are supplemented by the geometrical data of the segment which joins the node to its proximal neighbour.The base grid is generated by the software GROGRA.We can use two different sources for base grid generation: (1) Empirical information from the total mapping of single trees, (2) tree structures which are generated by GROGRA on the basis of growth grammars.
The multi-step discretization algorithm realiz-es a compromise between the numerical requirements of uniform accuracy and locally limited variation of step-widths (see Friih 1995).According to their local leaf specific conductivities, the base grid segments are subdivided in such a way that approximately the same steady-state V Fdifference A?¥f tat drops along all segments of the tree.We call A¥^f l , the steady-state accuracy level with given E. Finally, in the following procedure the spatial pattern of step-widths is smoothened.
For temporal discretization an algorithm was taken from Hornung and Messing (1984) which adjusts time steps such that the change of water content 0 per step approximates a desired optimum A9 f>/ ,,, where a halving of A9 op , on average leads to a halving of time steps.

Results of Model Tests
The convergence property of the finite difference method of Douglas and Jones (1963) is mathematically proven, but a number of specific adjustments were done, and about 12 000 lines of the computer program HYDRA are between the sets of difference equations and the final simulation results.Therefore model tests should give additional evidence that the physical and physiological assumptions are realized correctly by the simulation method.
In most of the model tests presented, steps or rectangles of transpiration flux density were given as input signals.The sudden changes between the zero level and 5.0 • 10" 6 kg/(s • m 2 ), which corresponds to a normal daily transpiration maximum of Picea abies, were broadened by HY-DRA to 60.0 s.This signal type is preferred for two reasons: it demands high standards of the numerical method, and the step response characterizes system behaviour quite comprehensive (Oppenheim and Willsky 1983).

Spatial Convergence
In the Picea abies growth grammar, sketched by Kurth and Lanwert (1995), stochastic variations were eliminated and the lengths of branching axes extended by a factor of 1.5.From this a 4-year-old model tree was generated and transformed to a base grid.It served as the common ancestor for the generation of 4 hydraulic maps differing from each other only in terms of spatial discretization, while the branching and terminal nodes are identical.For the 4 maps, the steady-state accuracy levels AYfa,, = 4.0 • 10~3,2.0• 10~3,1.0• 10" 3 and 5.0 • 10"* MPa were chosen, which in the mean corresponds to a progressive halving of spatial steps.The control parameter A0 (;/)/ for temporal discretization was set to 10" 4 kg/dm 3 , ^F^a was 0.0 MPa, thermodynamic equilibrium as initial condition, E was a rectangle signal with a breadth of 1 h located within a total time interval of 2 h.
We found that with each halving of space steps the relative ^-differences stay within 4 % at the half-times of responses (ca. 3 min after stepwise changes), and within 0.2 %, when the steady-state flow maxima are reached.The differences did not change significantly from the first to the third halving of steps.

Temporal Convergence
From the trees in the tests for spatial convergence, the hydraulic map with the steady-state accuracy level ATf^ = 10" 3 MPa for E = 5.0 • 10~6 kg/(s • m 2 ) was taken.4 simulations were run with AG opt = 10-3 , 5-10^, 10" 4 10" 5 kg/dm 3 .In the mean this approximately corresponds to a progressive halving of time steps, which however may change during an individual simulation by 3-4 orders of magnitude.Transpiration, total time interval, initial and boundary condition were set as in the tests for spatial convergence.
We observed that with each halving of time steps the relative ^F-differences stay within 2 % at the half-times of responses, and within 0.01 % when the steady-state flow maxima are reached.The differences did not change significantly from the first to the third halving of steps.

Routine Checking of Mass Balance
During each simulation, HYDRA automatically generates a protocol of mass balance in which for individual time steps the input-output bal-ance of the tree is compared to its internal mass change.Independently of tree architecture and temporal course of transpiration, the differences observed never exceeded the order of machine accuracy (2.2 • 10~1 6 ).The simulation protocols also confirm that in steady-state situations the input-output balance and the internal mass change equal 0 (within the scope of machine accuracy).

Sensitivity to Changes of Conductivity, Capacity and Leaf Area
By means of the Picea abies growth grammar (Kurth and Lanwert 1995) a 15-year-old model spruce was generated and transformed to a hydraulic map.It served as a control tree with leaf areas as given by the growth grammar and conductivities and capacites following the data mentioned in the section about parametrization.By means of a global scaling factor of 2.0 resp.0.5, 6 tree variants were generated which differed in the levels of one of the three properties conductivity, capacity and leaf area.The steady-state accuracy level AH'f^ of all test trees was within 4.0 • 10" 3 MPa, the control parameter A0 o/7/ was set to 1(H kg/dm 3 .¥ 8oa was 0.0 MPa, as initial condition thermodynamic equilibrium was chosen.At starting time, E step-wisely changed from 0 to 5 • 10~6 kg/(s • m 2 ), the total time interval was 7 h.The step response of water flow at the stem base was used to characterize the system behaviour of the test trees.We found that a doubling of one of the three quantities 1 / K, C,L led to raisings of half-times by factors < 3.0, while the halving of a quantity caused reductions of half-times by factors > 0.4.

Numerical Stability
The same hydraulic map was used as in the tests for temporal convergence.The control parameter AQ opt for temporal discretization was set to 10^, ^soii was 0.0 MPa.For the stability tests, perturbation terms were superimposed on the initial condition, which normally should assign ^Fsoii to all spatial grid points Xj.The perturbation terms had the form a • r,(-1.0,1.0), where the function r/-1.0,1.0) generates random numbers evenly distributed within the interval [-1.0, 1.0] and a is the perturbation amplitude.At starting time E step-wisely changed from 0 to 5 • 10~6 kg/ (s • m 2 ), the total time interval was 2 h.
The simulated potentials were compared to the results of a non-perturbed simulation.The total X Frange covered during simulations was about 4 • 10" 2 MPa.A perturbation amplitude a = 5 • 10~3 MPa was dampened to fractions between 6 and 1 % at the half-time of step response (ca. 3 min after start of simulation), and at the steady-state end of simulations to fractions between 2.5 and 0.025 %.Perturbations with a > 5 • 10" 3 MPa decayed only incompletely, which can be traced back to irreversible conductivity losses.

A-posteriori Error Analyses
With the test simulations for numerical stability, additionally a-posteriori error analyses were carried out (Watkins 1991) to evaluate how accurately the linear sets of difference equations are solved.The analyses revealed that the residua of the sets do not exceed machine accuracy (2.2 • 10" 16 ) by more than 1 order of magnitude.

Discussion
From the results of model tests for spatial and temporal convergence we conclude that A^f^ < 4.0 • 10-3 MPa and AQ opt < 10" 3 kg/dm 3 are scopes of the control parameters of step widths appropriate for simulation, because a further reduction of spacings leads to no significant change of the solution function.Mass conservation within the limits of machine accuracy is a consequence of our conservative difference equations and by no means a matter of course: finite element models of groundwater or soil hydrology may exhibit deviations from mass conservation up to several percents (Istok 1989).From the comparison of the hydraulic dynamics of test trees differing in their levels of conductivity, capacity or leaf area, we conclude a moderate sensitivity of the model to changes of these quantities.The rapid decay of perturbations, which were superimposed with high amplitude on the initial state, points to a high numerical stability of the model and thus reproducibility of the results.The tests for sensitivity and stability indicate that the matrices of the linear sets of equations in the discrete IBVP have a small condition number (Watkins 1991).Thus we can interprete the tiny residua of the aposteriori error analyses to the effect that the accuracy of the computed solutions is limited only by round-off errors due to computer hardware.Summarizing the test results, we have strong evidence that the physical and physiological assumptions, which we condensed in our way from the continuous starting equation up to the discrete IBVP, are realized accurately and reproducibly by the computational method.
HYDRA simplifies aspects concerning the integration of trees within the soil-plant-atmosphere-continuum, which cannot be mirrored by the tree water flow model alone.However, the approach is designed such that spatially variable transpiration flux densities and the stomatal control can be included without change of the numerical method.A realistic spatial distribution of transpiration is started by coupling HYDRA to a 3-dimensional model of microclimatics (Lanwertetal. 1987).
The reason why we confined the domain of water flow modelling to the above-ground tree parts is that, because of its strong feedback to the hydraulic variables of the soil, root water-uptake cannot be separated from the water movement within soil.In turn, there is increasing evidence that biochemical signalling mediates a dependency of stomatal conductance on fine root water status (Davies and Zang 1991, Gollan et al. 1985, Jones and Sutherland 1991, Turner et al. 1985).To meet these interdependences decisive for hydraulic integration, we prepared a coupling of HYDRA to SilVlow, a numerical soil water flow model (Hörmann and Schmidt 1995).The preconditions are favourable: both models reflect the flow paths of water, they are based on similar equations and discretization follows similar criteria.We are on the way to complete the theoretical synthesis of the hydraulic system of trees.Effective solution algorithms, the parametrization principle, standardized space discretization and the linkage to the morphological model GROGRA (Kurth 1994a), which can generate and vary tree architectures, predestine HYDRA as a tool to study the dependency of hydraulic function on structure.

Fig. 1 .
Fig.1.Spatial discretization scheme in the vicinity of a branching node.The shaded region, composed of segment halves joining at the grid point xj, denotes the cell for which mass conservation is fulfilled by eq.(6).Arrows symbolize flows across the borders of the cell and transpirational water loss.