A simple three-dimensional macroscopic root water uptake model based on the hydraulic architecture approach

Many hydrological models including root water uptake (RWU) do not consider the dimension of root system hydraulic architecture (HA) because explicitly solving water flow in such a complex system is too time consuming. However, they might lack process understanding when basing RWU and plant water stress predictions on functions of variables such as the root length density distribution. On the basis of analytical solutions of water flow in a simple HA, we developed an “implicit” model of the root system HA for simulation of RWU distribution (sink term of Richards’ equation) and plant water stress in three-dimensional soil water flow models. The new model has three macroscopic parameters defined at the soil element scale, or at the plant scale, rather than for each segment of the root system architecture: the standard sink fraction distribution SSF , the root system equivalent conductance Krs and the compensatory RWU conductanceKcomp. It clearly decouples the process of water stress from compensatory RWU, and its structure is appropriate for hydraulic lift simulation. As compared to a model explicitly solving water flow in a realistic maize root system HA, the implicit model showed to be accurate for predicting RWU distribution and plant collar water potential, with one single set of parameters, in dissimilar water dynamics scenarios. For these scenarios, the computing time of the implicit model was a factor 28 to 214 shorter than that of the explicit one. We also provide a new expression for the effective soil water potential sensed by plants in soils with a heterogeneous water potential distribution, which emerged from the implicit model equations. With the proposed implicit model of the root system HA, new concepts are brought which open avenues towards simple and mechanistic RWU models and water stress functions operational for field scale water dynamics simulation.


Introduction
Plants impact the terrestrial water cycle, in particular, through evapotranspiration.Root water uptake (RWU) affects underground water dynamics, with consequences on plant water availability and groundwater recharge.However, even though hydrological and climate models are sensitive to RWU and plant water stress parameters (Desborough, 1997;Zeng et al., 1998), no consensus exists on the modelling of these two processes (Feddes et al., 2001;Skaggs et al., 2006;Raats, 2007).
From a conceptual point of view, two main approaches exist today, which differ in the way they predict the volumetric rate of RWU, or "sink term", of Richards' equation in volume elements of soil: where θ is the volumetric water content (L 3 L −3 ), t is the time (T), K is the unsaturated soil hydraulic conductivity (L 2 P −1 T −1 ), H s is the total soil water potential (P), which will be referred to as the "soil water potential", and S is the sink term (L 3 L −3 T −1 ).Note that the units of K and H s differ from standards of soil physics, but were chosen for consistency with those used in plant physiology.

V. Couvreur et al.: A macroscopic root water uptake model based on the hydraulic architecture approach
The first approach promotes a detailed, physically-based modelling of water flow, from the soil-root interfaces to the plant collar, inside the three-dimensional root system hydraulic architecture (HA) whose segments hydraulic properties can be defined individually.This approach originated from the proposition of Van Den Honert (1948) to express water transport in plants as a catenary process (i.e., the partial process that encounters the highest resistance governs the velocity of the whole), of which Ohm's law equation is an appropriate mathematical expression.This approach was later developed by Landsberg and Fowkes (1978), Lhomme (1998) or Doussan et al. (1998a) to describe water flow in root system architectures.Coupled with a threedimensional soil water flow model, it led to quite sophisticated RWU models at the plant scale (Doussan et al., 2006;Javaux et al., 2008;Schneider et al., 2010).Such models may predict compensatory RWU, hydraulic lift and RWU under water stress conditions without any additional feature than hydraulic principles.Moreover, each of their parameters (root system architecture, root segment radial conductivity and axial conductance to water flow) has a physical meaning and can be measured directly.However, the cost of characterising these parameters is a major drawback when using such models.In addition, this type of model is very demanding in terms of computational power and time, which explains why it cannot be used at crop management relevant scales (Schröder et al., 2009).
The second approach, generally favoured in crop management models and soil-vegetation-atmosphere transfer models (SVATs), relies on model parameters defined at the soil element scale and at the whole plant scale; these parameters will be referred to as "macroscopic parameters".This type of model usually predicts RWU as the product of the potential transpiration rate (T pot ) by a spatially distributed root parameter (e.g., the relative root length density -rRLD), a stress function depending on the local soil water/osmotic potential (Feddes et al., 1976), and sometimes a compensatory RWU function (Jarvis, 1989).Despite its simplicity and potential efficiency, this approach is also subject to criticism.First, most of the macroscopic parameters cannot be directly determined or measured and, thus, require a calibration.This calibration stage is subject to major limitations: low sensitivity and non-uniqueness of the model parameters, lack of extrapolation power and uncertainty on the measurements used for the calibration (Musters and Bouten, 2000;Hupet and Vanclooster, 2005;Vandoorne et al., 2012).Second, by using root length or mass density distributions, these models neglect the effect of root hydraulic properties and architecture while numerous authors show their significant influence on RWU (Pierret et al., 2006;Schneider et al., 2010).Finally, predicting RWU as the product of T pot by other factors forces all local RWU rates to nullify when T pot is null, which prevents the simulation of hydraulic lift while this process is proven to occur during low transpiration rate periods (Dawson, 1996;Song et al., 2000).
Far from trying to emphasise the best approach, Feddes et al. (2001) encouraged continuing the development of both modelling approaches by increasing the complexity and completeness of existing physically-based models (for accuracy and scientific understanding of the modelled process) while keeping macroscopic RWU models as simple as possible (so that appropriate computational weight could be paid to each modelled process, depending on its importance).
How to improve a RWU model while keeping it as simple as possible is a complex task.Recently, Raats (2007), De Jong Van Lier et al. (2008) and Jarvis (2011) attempted to do so with the following approach: deriving a macroscopic RWU model from an approximate analytical solution of a detailed RWU model.Yet, all of these models tended to neglect the effect of the root system HA.
In this paper, we developed a macroscopic RWU model based on analytical solutions of water flow in a simple HA and validated it for a more complex HA.A new spatially distributed root parameter instead of the rRLD and new stress and compensatory RWU functions emerged from this approach, which call for a complete revision of classical RWU models.

Theory
We first considered a simple root system as the analogue of a circuit of resistances to water flow (like in the first, detailed, approach described in the introduction), and derived analytical solutions of its water relations by solving the system of equations of Ohm's and Kirchhoff's laws.Figure 1 shows the simple root system HA (discretised in point elements called root nodes, conceptually divided into a root xylem node (green circle) and the associated soil-root interface (brown circle); the root nodes being connected by root segments) together with the variables used to express water relations in this system: the water potential in xylem vessels at the plant collar H collar (P), which will be referred to as the "plant collar water potential", soil-root interfaces water potentials H sr (P), root xylem nodes water potentials H x (P), axial resistances to water flow R x (P T L −3 ) between two root xylem nodes, radial resistances to water flow R r (P T L −3 ) between a root xylem node and the associated soil-root interface (note that in the following developments, R x , R r and the root architecture are supposed not to change with time), root axial water flow rates Q x (L 3 T −1 ), root radial water flow rates Q r (L 3 T −1 ) and the actual transpiration rate T act (L 3 T −1 ).

Shape of the simple root system hydraulics model
A system of twelve equations can describe water flow rates in the HA shown in  H x,4 and isolating the RWU rates Q r,i , analytical solutions of Q r,i at each root node can be found, which have a common structure: ( Three terms appear in this structure: (i) the actual transpiration rate T act , (ii) ϕ i (L 3 T −1 ), which is the only term that depends on the soil water potential distribution, and (iii) f i (−), which is a fraction depending on the root axial and radial resistances and on the root architecture.
According to its properties, ϕ i represents the "compensatory RWU" process at the i-th root node (ability of the root system to adapt its uptake distribution in response to the soil water potential distribution).For any uniform water potential at the soil-root interfaces (uniform H sr vector), there is no compensatory RWU, i.e., ϕ i = 0 at each root node.In case of positive ϕ i , the uptake rate at the i-th root node is increased of a value that equals ϕ i f i , as compared with a situation where H sr is uniform.For a negative ϕ i , the uptake rate at the ith root node is reduced as compared with the uniform case.When ϕ i < −T act , the uptake rate is negative and water flows from the root into the soil, like when hydraulic lift occurs.
When there is no compensatory RWU, which will be considered as the "standard condition", f i equals the uptake fraction at the i-th root node.We, therefore, propose to call it the "Standard Uptake Fraction at the i-th root node", or shortly, "SUF i " (−).All together, the SUF i form a vector (SU F ) whose sum is 1.

Expression for the compensatory root water uptake
By gathering the variables H sr,i in a specific way in the analytical expression of ϕ 1 , it comes out that several coefficients can be factorised into the exact expressions of the SUF i divided by a constant term: which is equivalent to: H sr,j SUF j (4) where ρ (P T L −3 ) depends on root axial and radial resistances (see Eq. A2), N is the total number of root nodes, j ranges from the first to the last root node, and i is 1.While Eq. ( 4) also applies to ϕ 2 (see Eq. A6), this is not the case for ϕ 3 and ϕ 4 (see Eqs. A7 and 5): where However, when R x,i R r,j (resistances to root axial water flow much lower than resistances to root radial water flow), α ≈ 1 and Eq. ( 4) can be generalised to ϕ 3 and ϕ 4 .Under this condition, the coefficients of the relation linking ϕ i to H sr,i and to the other H sr is constant in time and uniform for all root nodes.Hereafter, we consider that Eq. (4) applies at any i-th root node.
If we consider Eq. ( 4) as an Ohm's law equation, where "the flux" is the compensatory RWU rate at the i-th root node (ϕ i ), and "the potential difference" is the difference between the i-th soil-root interface water potential (H sr,i ) and a weighted-mean water potential of all soil-root interfaces H sr,j SUF j ), then the factor 1 ρ can be considered as the

V. Couvreur et al.:
A macroscopic root water uptake model based on the hydraulic architecture approach effective conductance of the compensatory RWU process.Therefore, we propose to call this parameter "compensatory RWU conductance" and to denote it by the symbol "K comp " (L 3 P −1 T −1 ).Just like ρ, K comp only depends on root axial and radial resistances to water flow.By combining Eqs. ( 2) and (4), we, therefore, obtain the following complete equation for the simple architectural RWU model: It can be demonstrated that Eq. ( 7) is exact for RWU prediction in any root system with negligible root axial resistance.

Water stress function
In the previous paragraph, we derived the uptake distribution in the root system when T act is given.In hydrological models, when the soil water potential is not a limiting factor for RWU, T act is supposed equal to the so-called potential transpiration rate (T pot ), which depends on atmospheric conditions and leaf properties (Van Den Berg et al., 2002).When the plant roots cannot sustain the atmospheric demand for transpiration, isohydric plants control their stomatal conductance in order to keep the water potential in their leaves at a threshold value.Under these conditions, we assume the water potential in the leaves, and consequently at the plant collar (H collar ), to remain constant over time.This implies that T act needs to be calculated from a prescribed H collar .
The system of twelve equations that describes water flow in the HA represented in Fig. 1 also links T act , H collar , H sr and the root hydraulic properties (R x and R r ).After eliminating and H x4 , the following equality can be demonstrated: where K rs is the "equivalent conductance of the root system" (L 3 P −1 T −1 ), or inversed Thevenin equivalent resistance of the root resistance network linking the plant collar to the soil.For more detail, the analytical solution that led to Eq. ( 8) is given in Appendix (Eq.A8).
Two interesting features appear in Eq. ( 8).First, the factor multiplying each soil-root interface water potential H sr,j is the exact analytical expression of the corresponding SUF j .Second, the factor K rs can be calculated from Thevenin theorem (Thévenin, 1883), which is commonly used in electric circuits, to calculate the equivalent resistance of a resistance network.
The relation between T act and H collar allows the use of the plant collar water potential as stress indicator.As long as no stress occurs (i.e., H collar predicted from Eq. 8 is higher than a threshold value), T act = T pot .In case stress occurs (i.e., H collar is lower than the threshold value), H collar is fixed at the threshold value and T act is predicted from Eq. ( 8).
It is worth noting that under the simplifying hypothesis that allowed the generalisation of Eq. (4) (R x,i R r,j ), the analytical expression of K rs becomes equal to that of K comp , which would reduce the number of parameters of the simple RWU model and water stress function to two (SU F and K rs ).

Expression of the simple root system hydraulics model at the soil element scale
The simple RWU model developed in Eqs. ( 7) and ( 8) was set up for the root system architecture domain and, thus, cannot directly apply at the soil element scale at which the sink term of Richards' equation is defined.In this section, we focus on the conversion of Eqs. ( 7) and ( 8) into macroscopic equations (i.e., defined at the soil element scale).
The following developments first require a clear definition of what the soil nodes and elements are.Like the root system architecture, which is discretised in root nodes, the soil domain is discretised in soil nodes (point elements at which the soil water potential is defined).Adjacent soil nodes delimit a volume (for instance a cube), which is called a soil element.
Practically, when simulating RWU by an explicit HA, the soil domain is linked to the root architecture domain in two ways: (i) the RWU rate in a soil element is the sum of the RWU rates of the root nodes that it contains; (ii) the soil-root interface water potential is calculated as a spatial interpolation of the surrounding soil nodes water potentials.
Link #1 is mathematically expressed as follows: where S k (L 3 L −3 T −1 ) is the volumetric RWU rate in the kth soil element, V k (L 3 ) is the volume of the k-th soil element, ε ik (−) is a factor which is 0 except when the k-th soil element contains the i-th root node (then it equals 1), and Q r,i (L 3 T −1 ) is the RWU rate at the i-th root node.Note that for a given i, only one value of ε ik is different from zero (i.e., a root node is contained in only one soil element).Conversely, for a given k, a variable number of values of ε ik can be different from zero (i.e., a soil element can contain no, one, or several root nodes).
In order to obtain an analogue of Eq. ( 7) at the soil element scale, link #2 is approximated by the following assumption: "the water potentials of all soil-root interfaces located within the k-th soil element equal the averaged water potential of the soil nodes of the k-th soil element (H s,k )", which is mathematically expressed as follows: For clarification, Eq. ( 10) can be summarised as H sr,i = H s,k when the i-th root node is inside the k-th soil element (then, ε ik = 1).
Hydrol.Earth Syst.Sci., 16, 2957-2971, 2012 www.hydrol-earth-syst-sci.net/16/2957/2012/By combining Eqs. ( 7), ( 9) and ( 10), we obtain: where M is the total number of soil elements and SSF k (−) is the value of N i=1 ε ik SUF i , which we propose to call the "Standard Sink Fraction in the k-th soil element" for similar reasons as for SUF i .Since the standard sink fraction distribution is simply a different partition of the standard uptake fraction distribution, its properties are similar, but in the soil domain instead of the root architecture domain: (i) the SSF k depend on the root system architecture and its hydraulic properties, and (ii) all together, the SSF k form a vector (SSF ) whose sum is 1.It can also be demonstrated that, by using the same approximation of link #2, an analogue of Eq. ( 8) at the soil element scale is found: The model developed in Eqs. ( 11) and ( 12) allows the simulation of RWU from macroscopic parameters defined at the soil element scale (SSF ) and at the plant scale (K comp , K rs ).
In the next pages, it will be referred to as the "model implicitly accounting for the root system HA" (or shortly, "implicit model"), because its inputs and outputs are defined at the soil element scale, and no root system architecture has to be considered during a simulation using this model.However, the model structure, and its parameters, keep an "implicit" footprint of the root system HA, because they were deduced from solutions of water flow in a HA.It is notable that under uniform soil water potential conditions, Eqs. ( 11) and ( 12) simplify into the following equations:

Methodology
The validity of the implicit model for realistic root systems relies on three hypotheses: (i) the equations developed with the simple HA apply for more complex root system HA (i.e., the solutions of water flow equations have the same structure in complex HA, independent of the simple HA shown in Fig. 1), (ii) root axial resistances to water flow values are low enough to neglect their effect on compensatory RWU, (iii) all the root nodes located within a certain soil element have a soil-root interface water potential which can be approximated by the averaged water potential of the soil nodes of the soil element.
In order to explicitly simulate water flow in more complex root system HA, we used the Doussan model (which will be referred to as the "explicit model"), which allows solving numerically the system of equations of Ohm's and Kirchhoff's laws in any HA by inverting a set of linear equations.
As explained in more detail in the following sections, the explicit model was used to parameterise the implicit model.After that, in order to validate the implicit model, two series of tests comparing its predictions with those of the explicit model were carried out.First, "instantaneous tests" were used to check the existence and properties of the macroscopic parameters for a maize root system HA.Second, "long term tests" were used to quantify the accuracy of the implicit model, when coupled to Richards' 3-D soil water flow equations, over a longer time period.
3.1 Description of the complex maize root system and of the simulation domain

Root system architecture and hydraulic properties
Our objective was to generate a maize root system whose hydraulic and geometric properties were as close as possible to reality.An 80-days old maize root system of 35 000 nodes was generated with RootTyp (Pages et al., 2004).This code generates root systems by taking into account plant-specific genetic properties like insertion angles of the different root types, their trajectories, average growth speed and distances between lateral roots.The corresponding values were parameterised based on information from Tardieu and Pellerin (1990) and Girardin (1999).The environmental Root-Typ parameters were optimised in order to fit measured root length density profiles from a maize field (Tardieu, 1988).
Figure 2a shows the optimised root system architecture.
Variable maize root hydraulic properties evolving with root segment age and type were obtained from Doussan et al. (1998b).Based on root average growth speed (Girardin, 1999), principal root hydraulic properties could be expressed as a function of segment age instead of distance from root tip (see Fig. 3).The chosen threshold plant collar water potential at which stress occurs was −15 000 hPa.

Simulation domain properties
In order to represent a field root distribution while limiting the computational needs, the generated root was located in a periodic soil domain of 15 cm (direction of the maize rows) on 75 cm (direction perpendicular to the maize rows).This domain was periodic at its vertical boundaries for soil water fluxes, root system architecture and root water fluxes.No flux boundary conditions were imposed at the top and bottom of the soil domain.The depth of the soil domain was 124.5 cm and the spatial discretisation 1.5 cm.
Figure 2b shows the boundaries of the periodic domain (in red), the central root architecture (in black) and the root branches that crossed the vertical boundaries of the periodic domain (in green).

Existence and properties of the macroscopic parameters for the complex root system
According to their analytical expressions in the simple HA, the macroscopic parameters SSF , K rs and K comp should be constant as long as the root system architecture and hydraulic properties do not change with time.In this section, we further checked whether such parameters exist in the complex HA by calculating them under different conditions supposed not to affect their values.

Standard sink fraction distribution and root system equivalent conductance
The SSF and K rs were obtained by respectively using Eqs.( 13) and ( 14) in which the sink terms and plant collar water potential were numerical solutions of Doussan equation for uniform soil water potential conditions.
Both parameters were first calculated for the following conditions: uniform soil water potential of −150 hPa and actual transpiration rate of 1200 cm 3 d −1 .Whether these reference parameters also apply for other uniform soil water potentials (−50, −500, −1500 and −5000 hPa) and actual transpiration rates (1, 600, 1800 and 2400 cm 3 d −1 ) was subsequently checked.

Compensatory root water uptake conductance
If Eq. (11) applies to the complex HA, at a given time, the vector of compensatory RWU ϕ should be a linear function of the vector of soil water potentials H s : Note that the second term of the difference on the right-hand side of Eq. ( 15) is a scalar value.slope, whose value would be K comp , but might have different intercepts.
The linearity of the H s − ϕ relation and the constancy of its slope, i.e., K comp , was evaluated at several times during a scenario with a dynamic water content (scenario "Equil" in Table 1).We used the R 2 of ϕ k versus H s,k plots as linearity criterion.The ϕ k were obtained by (i) solving the Doussan equation, which gave us S k under non-uniform soil water potential conditions, and (ii) using the following equation: In addition, the influence of the order of magnitude of the root axial conductances on the linearity of the H s − ϕ relation was investigated.We also checked for which root axial conductances K comp could be approximated by K rs .

Validation of the implicit model
The "long term test-scenarios" of water uptake by the maize root system were run with both explicit and implicit models.
To quantify the impact of the compensatory RWU, the implicit model was also run with K comp = 0, representing a simulation without compensatory RWU.To evaluate the accuracy of the implicit model as compared to the explicit model, the time evolutions of the mean absolute differences (MAD) between sink terms and water contents simulated by both models were calculated for each scenario.In parallel, the coefficients of determination (R 2 ) between plant collar water potentials and transpiration rates simulated by both models were calculated.
Five one-week scenarios, with different soil hydraulic properties (from Carsel and Parrish, 1988), initial soil matric potential profiles and root hydraulic properties, were selected (cf.Table 1).
Note that neither root architecture nor root hydraulic properties changed with time.Root hydraulic properties were defined on the basis of root segments ages at the beginning of the scenarios.
A time series of 600 cm 3 day −1 plant −1 sinusoidal day/night T pot was chosen as root boundary condition.This corresponds to a reference evapo-transpiration (ET 0 ) of 4.5 mm day −1 , which is typical for a warm Belgian summer day (Baguis et al., 2010), under a well-developed maize crop (K c = 1.2) where the surface per plant (Surf) is 15 × 75 cm.The relation used to predict the daily potential transpiration rate (L 3 T −1 ) was: T daily = ET 0 K c Surf.

Existence and properties of the macroscopic
parameters for the complex hydraulic architecture

Standard sink fraction distribution and root system equivalent conductance
The macroscopic parameters SSF and K rs characterised under the reference conditions matched those calculated for all of the tested conditions respectively with R 2 of 1.0000 and absolute difference percentages lower than 0.014 %.These results support the existence of the macroscopic parameters SSF and K rs for the complex root system HA.In Fig. 4, we can see that decreasing the root axial conductances of principal roots leads to a shift of SSF to regions of the root system that are closer to the plant collar, whereas the opposite is the case when the radial conductances of young roots are increased.It is notable that all the shown SSF differ from the relative root length density vector (rRLD).
Figure 5 shows that decreasing the root radial or axial conductances leads to a nonlinear decrease in K rs and its value is more affected by a decrease in root radial conductances.
The observed properties for both SSF and K rs of the complex root system are in agreement with those deduced from their analytical expressions for the simple root system: they are sensitive to root hydraulic properties and not to T act and H s .

Compensatory root water uptake conductance
Figure 6 illustrates the spatial relation between ϕ (predicted by the explicit model) and both H s and root hydraulic properties in a vertical slice of a soil profile with vertical gradient of soil water potential.It shows that the gradient of compensatory RWU spatially follows the gradient in soil water potential.At locations where H s,k is lower, i.e., at the top of the profile in this example, ϕ k is negative, meaning that the uptake rate is reduced as compared to the uniform soil water potential case.The reduced uptake at the top of the profile is Table 2. Comparison of the values of K comp and K rs for different root axial conductance levels of magnitude (as compared to the reference root hydraulic properties, Fig. 3).
Factor multiplying the reference axial conductances (Fig. 3 compensated by an increased uptake at the bottom of the profile (i.e., ϕ k is positive), where H s,k is higher than average.Equation ( 15) is in agreement with this observation, together with the fact that the compensatory RWU is reduced in intensity when the root radial conductance values (and, thus, K comp ) are lower.
Figure 7 shows the ϕ k versus H s,k plots obtained from the "Equil" scenario for different days (coloured crosses).It is observed that the slope of these plots (i.e., K comp ) does not change significantly at different time steps.In addition, the straight lines predicted with the implicit model fit the ϕ k versus H s,k plots.
In Fig. 8, we checked the linearity of this relation by computing the determination coefficient (R 2 ) of the ϕ k versus H s,k plots obtained from the "Equil" scenario.For the reference root hydraulic properties (yellow), R 2 is always higher than 0.945.Figure 8 also shows a significant loss in R 2 (which can be explained by scattering or nonlinearity) when the root axial conductances decrease while the gain in R 2 is slight for the same increase in root axial conductances.
These two tests support the validity of Eq. ( 15) for the complex root system architecture and reference root hydraulic properties by (i) attesting the linearity of the H s − ϕ relation and (ii) showing the invariability of the macroscopic parameter K comp .Yet, the user should avoid applying Eq. ( 15) for root systems with low root axial conductance values.
Table 2 shows the relative difference between K comp and K rs for a broad range of root axial conductances values (the root radial conductances being kept constant).It confirms that they can be considered as one single parameter for relatively high root axial conductances (100 times these of Fig. 3), but not in the other cases.Therefore, K rs and K comp were considered as two independent parameters when using the implicit model in the soil water dynamics scenarios.

Sink term and water content distribution predictions
Figure 9a shows the evolution of the MAD of sink terms predicted by both explicit and implicit models for the "CL" scenario, with compensatory RWU (green line) and without compensatory RWU (red line).The mean absolute sink term (blue line) is given as reference for comparison.The same is shown in Fig. 9b, for water content distributions, with the mean loss in water content from the beginning of the scenario (blue line) as reference for comparison.
The MAD on the sink term prediction is globally lower than 2 % of the mean absolute sink term while it reaches 20 % of the mean absolute sink term when the compensatory RWU process is neglected.
As shown in Fig. 9b, these differences on the sink term prediction imply differences on the water content: the MAD on the water content is lower than 1 % of the mean water content loss (which can be defined as the ratio between the cumulated transpiration and the volume of the soil domain) while it reaches 10 % of the mean water content loss when the compensatory RWU process is neglected.
Globally, all MAD and maximum differences results (not shown) are close to those of the "CL" scenario with best performance for the "Inv" scenario and worst performance for the "SCL" scenario.We also notice that the maximum differences on the water content predictions are all lower than 1 % in absolute water content when the compensatory RWU process is considered while they vary between 3 and 10 % in absolute water content when this process is neglected.

Plant collar water potential and actual transpiration rate predictions
Determination coefficients of H collar and T act predicted with the implicit model versus the explicit model equal 1.0000 in all scenarios.This result confirms the accuracy of the relation described by Eq. ( 12).
Figure 10 shows that reductions in both T act and H collar are stronger when no compensatory RWU occurs.

Computing time
For each of the five scenarios, the computing time with the implicit model was a factor 28 to 214 shorter than with the explicit model.This time reduction results from the fact that the Doussan matrix inversion, which is computationally heavy, particularly for root system architectures with high number of root nodes, is not needed anymore.
This time consumption reduction is an additional step towards high precision water dynamics modelling at larger scales.

Shape of the implicit root water uptake model
The shape of the implicit RWU model given in Eq. ( 2) was directly derived from analytical solutions of water flow in a root system HA, conferring a strong physical basis to our model.This shape is different from the common "product of factors" in which the sink term is proportional to the transpiration rate, to the root length density, and to factors that express the effect of local conditions (e.g., local water potential) and compensation mechanisms on RWU (Raats, 1974;Feddes et al., 1976;Molz, 1981;Prasad, 1988;Jarvis, 1989;Lai and Katul, 2000;Simunek and Hopmans, 2009).In the proposed model, the RWU rate (S k V k ) is the superimposing (i.e., sum) of two processes.First, the "standard RWU" term (T act SSF k ), which accounts for RWU in a soil with uniform water potential and is driven by the Hydrol.Earth Syst.Sci., 16, 2957-2971, 2012 www.hydrol-earth-syst-sci.net/16/2957/2012/ transpiration rate.Second, the "compensatory RWU" term (ϕ k SSF k ), which accounts for the internal adjustments of uptake rate distribution due to spatial variations in soil water potentials.This model split appears to be beneficial, not only in terms of accuracy when compared with the explicit model (see Sect. 4.2.1), but also since it may predict hydraulic lift as evidenced in Fig. 7: even though the transpiration is null at the selected times, we can see that it is not the case for most of the ϕ k .Water actually flows, from high water potential zones to low water potential zones, through the root system, in both explicit (crosses) and implicit (solid lines) models.Indeed, as the sink term is not proportional to the transpiration rate, and because the compensatory RWU term nullifies only when the soil water is at hydrostatic equilibrium, the implicit model allows water flow in the macroscopic root system when the transpiration rate is null.It is notable that, like in the macroscopic RWU model of Jarvis (2011), hydraulic lift can be considered as an extreme form of compensatory RWU.
A simple and exact expression of the compensatory RWU term could not be extracted from the analytical solutions of water flow in the simple HA.A single linear relation (Eq.4) between local compensatory RWU and local soil water potential, with constant slope, was chosen as an approximation of the exact solution.It was proven that this assumption is more accurate for higher than for lower root axial conductance values (Fig. 8), but shows satisfying results for realistic maize root hydraulic properties.As long as the main barrier to water uptake by plant roots is the radial transport path from root epidermis to xylem, rather than the axial path along xylem conduits (Frensch and Steudle, 1989;Steudle and Peterson, 1998), which is the case for a wide range of other plants, this assumption seems to be valid and the model applicable.

Shape of the implicit water stress function
In contrast to current RWU models in which compensatory RWU is closely linked to local water stress and local soil water potential (e.g., Feddes et al., 1976;Jarvis, 2011), in our model water stress is only determined by the water potential at the plant collar and is predicted independently of compensatory RWU.Compensatory RWU is driven by gradients in soil water potential close to roots and controlled by the conductance of the root system, but is independent of T act and of water stress and, therefore, is a phenomenon that is simulated continuously.Conversely, stress only occurs when T pot is high enough so that the soil water potential leads, for a given K rs , to a threshold plant collar water potential, which triggers water stress.The use of H collar rather than local soil water potentials as a regulator of transpiration is in agreement with its widely accepted role in transpiration regulation (Comstock, 2002;Christmann et al., 2007;Rodrigues et al., 2008).Even though active control by the plant on root radial conductances was not considered in our model, neither were osmotic stress or stress due to anoxia, the model structure allows including these effects in the model making use of biophysical process understanding.For instance, osmotic potentials in the root and soil water could be included in the total water potential.

Emerging macroscopic parameters
Three macroscopic parameters, which describe the plant macroscopic behaviour, emerged from our model.These parameters help understand the specific impact of root properties on RWU.
The standard sink fraction distribution (SSF ) is the normalised distribution of the sink term in the soil domain when the soil water potential is uniform.Its role is similar to that of the rRLD in Feddes RWU model.Both of them are proportional to RWU when there is no compensatory RWU and are independent of the soil water potential and transpiration rate (see Sect. 4.1.1 for the SSF ).However, the rRLD represents neither root architecture nor root hydraulic properties, while the SSF does (see Fig. 4 for sensitivity of the SSF to root hydraulic properties for the same root architecture and, therefore, same rRLD).Yet, root architecture and hydraulic properties are more difficult to derive from in situ measurements of plant roots than rRLD.However, this type of information can be obtained from root growth models calibrated on root profile measurements (as done in this study) or from noninvasive methods like MRI (Pohlmeier et al., 2008), X-ray or neutron tomography (Menon et al., 2007).The sensitivity of SSF to the root architecture type should be characterised in the future.
The root system equivalent conductance (K rs ) is the inversed Thevenin equivalent of the root resistances to water flow network linking the plant collar to the soil.In Sect.4.1.1,we showed that this macroscopic parameter still exists and keeps its explanatory power for the complex root system HA.The fact that K rs accurately links both T act and H collar to the soil water potential (Eq.12) makes it a key variable in the prediction of water stress for the plant.Figure 5 illustrates the sensitivity of this parameter to root hydraulic properties and evidences the fact that, in a realistic maize root system HA, K rs is more sensitive to radial than to axial root conductances to water flow.In addition, according to Thevenin theorem, K rs also depends on the root system architecture.
The compensatory RWU conductance (K comp ), whose value is close to that of K rs , but not equal to it for realistic root axial conductances (see Table 2), accounts for how intense the compensatory RWU, due to spatial variability of water potentials in the soil, will be.It was shown in Figs.7 and 8 that considering a constant K comp is a satisfying assumption for the wide range of soil water potential distributions covered by the "Equil" scenario, validating its existence and explanatory power for the complex maize root HA.We, thus, obtained a one-parameter compensatory RWU function, which is less than the number of parameters that are used in other compensatory RWU functions that have been proposed in the literature (generally two parameters whose value depend on the transpiration rate; Jarvis, 1989;Simunek and Hopmans, 2009).
Finally, the macroscopic parameters (SSF , K rs and K comp ) can be directly calculated from physical parameters, such as root hydraulic parameters and root architecture.For that, Doussan equation needs to be solved to obtain (i) S k and (ii) H collar under uniform soil water potential conditions, and (iii) S k under non-uniform soil water potential conditions.Then, by respectively using Eqs.( 13)-( 15), the macroscopic parameters can be found.If this physical information was not available, these parameters could as well be determined using inverse modelling.When using the common assumption of approximating the standard sink fraction distribution by the rRLD (or a function of depth), only two parameters (K comp and K rs ) need to be defined to model the standard RWU, the compensatory RWU, and the onset of transpiration reduction due to soil water limitations.Therefore, this model concept may lend itself better to a parameterisation by inverse modelling for root hydraulic property assessment (Draye et al., 2010), than in other RWU models of which the parameters depend, in addition, on the transpiration rate.

The soil equivalent water potential sensed by plants
Although numerous publications referred to Ohm's law for describing soil-plant-atmosphere water flows (Gardner and Ehlig, 1963;Feddes and Rijtema, 1972;Landsberg and Fowkes, 1978;Lhomme, 1998;Guswa, 2005), how to define an "equivalent soil water potential" when the soil water potential varies within the root zone was still a pending question.
Equation ( 12) governs water flow between the soil compartment and the plant collar, and its shape is that of an Ohm's law equation: the water flow equals the root system equivalent conductance times a difference in water potential.Since the first term in this difference is the plant collar water potential, the second term should be the water potential of the soil compartment.As a consequence, the expression of the equivalent soil water potential for soil-plant water flow should be the following: where H s,eq is the equivalent soil water potential sensed by the plant and equals the SSF -weighted mean soil water potential.Yet, H s,eq does not only drive plant water stress, but is also used to define the compensatory RWU in the k-th soil element, which is proportional to the difference between the water potential in the k-th soil element and H s,eq (see Eq. 11).

Accuracy of the implicit model
Globally, both implicit RWU model and stress function showed good performances in the tested scenarios (see Sects. 4.1.1 and 4.1.2) and the differences between the explicit and implicit models sink term predictions induced differences always lower than 1 % in absolute water content, which is lower than the uncertainty on water content measurements in cropped fields with currently used nondestructive methods (around 5 % for both ERT and GPR, and 2.5 % for TDR (Walker et al., 2004;Weihermüller et al., 2007;Brunet et al., 2009)), confirming that the implicit model can confidently be coupled to Richards' soil water flow equation for soil water dynamics simulations.
The fact that MAD between explicit and implicit models predictions were ten times higher when the compensatory RWU was neglected, evidenced the high significance of this process in soil water dynamics modelling, which is in agreement with several studies (Li et al., 2001;Dong et al., 2010;Jarvis, 2011).
Yet, the accuracy of the compensatory RWU prediction tended to decrease with increasing soil water potential heterogeneity (see the increasing scatter in Fig. 7), explaining the slight increase in the MAD at the scenario ends (see Fig. 9a).We also expect such a decrease in the implicit model accuracy when working with larger soil elements since the hypothesis on the uniformity of the soil-root interfaces water potential inside a soil element is expected to be less satisfying under these conditions, but this effect still has to be quantified.

Mathematical evidence of plant strategies against water stress
The simple water stress equation (Eq.12) and its macroscopic parameters (SSF and K rs ) can be used to postulate simple and intuitive strategies for plants to avoid water stress.Even though some of the suggested strategies may seem obvious, we think that mathematical evidence brings its own interest.From the root system point of view, maintaining high collar water flux at the limit of water stress (H collar fixed at −15 000 hPa) can be achieved by modifying two plant macroscopic parameters: (i) increasing K rs or (ii) increasing H seq by increasing SSF k at places where H s,k is high.For instance, the first can be achieved by globally increasing the rooting density (strategy observed by Li et al., 2010) or the root system conductance by the production of aquaporines (in agreement with Parent et al., 2009), and the second by favouring root exploration in wetter zones, or creating "tap roots", with high xylem conductance values, whose tips are established in wetter zones (in agreement with Ong et al., 1998).
In addition, Fig. 10 highlights the fact that the compensatory RWU process also plays a major role in delaying water stress, the absence of compensatory RWU process being translated in strong reduction in T act .

Conclusions and outlook
In this study, we propose a new model for RWU and water stress predictions for 3-D soil water flow models.The model parameters have the distinctive feature to be defined at the soil element scale, but to take implicitly into account the full root system HA.Based on analytical solutions of water flow in a simple root system, assuming high values of root axial conductances and local homogeneity of the soil-root interfaces water potentials, we derived two equations representing the general hydraulic behaviour of the root system (Eqs.11 and 12).In contrast to current RWU models, we propose a decoupling of water stress and of compensatory RWU processes and a model structure which is more appropriate for hydraulic lift simulation than the common "product of factors" structure.The suitability of these expressions was numerically proved for a more complex and realistic root system and the relevance of the compensatory RWU process in soil water dynamics modelling and in plant water stress prediction was emphasised.
We identified three macroscopic parameters controlling the macroscopic uptake process: (i) the standard sink fraction distribution (SSF ), (ii) the root system equivalent conductance (K rs ) and (iii) the compensatory RWU conductance (K comp ).These three characteristic properties can be obtained a priori based on Doussan equation for a given root architecture, root hydraulic properties and soil grid geometry or could be estimated by inverse modelling of water dynamics data.
In addition, a new definition was given to the effective soil water potential sensed by the plant in a soil with spatially heterogeneous water potential distribution.
As compared to current 3-D RWU models, the high computing speed of the proposed model opens new avenues to use mechanistic RWU models for inverse modelling and for large-scale water dynamics simulations.
Drawbacks of the proposed model are the fact that it was not made for predicting biological stress like that due to anoxia, its current limitation to plants having relatively large xylem conductance, and the fact that it should be coupled to 3-D Richards' equation resolution for soil water dynamics predictions.Note also that our simulations did not consider explicitly local hydraulic gradients in the rhizosphere zone around roots but this process can be accounted for by grid refinement (Schroder et al., 2009) and by choosing adequate local rhizosphere hydraulic properties (Carminati et al., 2011).

2 Figure 1 . 9 Fig. 1 .
Figure 1.Scheme of the simple HA and positions of the variables.Soil-root interface water 3 potentials H sr , root xylem nodes water potentials H x and the plant collar water potential H collar 4 are represented at the circles positions.(a) Axial resistances R x between two root xylem nodes 5 and radial resistances R r between a root xylem node and the associated soil-root interface.(b) 6Root axial water flow rates Q x , root radial water flow rates Q r , the actual transpiration rate T act 7 and their positive directions (arrows).8 9

Table 1 .
Scenario names and characteristics.