Benchmarking test of empirical root water uptake models

. Detailed physical models describing root water uptake (RWU) are an important tool for the prediction of RWU and crop transpiration, but the hydraulic parameters involved are hardly ever available, making them less attractive for many studies. Empirical models are more readily used because of their simplicity and the associated lower data requirements. The purpose of this study is to evaluate the capability of some empirical models to mimic the RWU distribution under varying environmental conditions predicted from numerical simulations with a detailed physical model. A review of some empirical models used as sub-models in ecohydrological models is presented, and alternative empirical RWU models are proposed. All these empirical models are analogous to the standard Feddes model, but differ in how RWU is partitioned over depth or how the transpiration reduction function is deﬁned. The parameters of the empirical models are determined by inverse modelling of simulated depth-dependent RWU. The performance of the empirical models and their optimized empirical parameters depends on the scenario. The standard empirical Feddes model only performs well in scenarios with low root length density R


Introduction
The rate at which a crop transpires depends on atmospheric conditions, the shape and properties of the boundary between crop and atmosphere, root system geometry, and crop and soil hydraulic properties. The study and modelling of the involved interactions is motivated by the importance of transpiration for global climate and crop growth (Chahine, 1992) as well as by the role of root water uptake (RWU) in soil water distribution (Yu et al., 2007). The common modelling approach introduced by Gardner (1960), referred to as microscopic or mesoscopic (Raats, 2007), is not readily applicable to practical problems due to the difficulty in describing the complex geometrical and operational function of the root system and its complex interactions with soil (Passioura, 1988). However, it gives insight into the process and allows development of upscaled physical macroscopic models (De Willigen and van Noordwijk, 1987;Heinen, 2001;Raats, 2007;De Jong van Lier et al., 2008. In many one-and two-dimensional problems, macroscopic RWU is modelled as a sink term in the Richards equation, whose dependency on water content or pressure head is usually represented by simple empirical functions (e.g. Feddes et al., 1976;Feddes et al., 1978;Lai and Katul, 2000;Li et al., 2001;Li et al., 2006;Vrugt et al., 2001). Most of these models are derived from the Feddes et al. (1978) model, which consists of partitioning potential transpiration over depth according to root length density and applying a stress reduction function of piecewise linear shape -defined by five threshold empirical parameters -to account for local uptake reduction. Results of experimental studies (Arya et al., 1975b;Clothier, 1995, 1999;Vandoorne et al., 2012) and the development of physically based models (De Jong van Lier et al., 2008;Javaux et al., 2008) increased the understanding of the mechanism of RWU as a non-local process affected by non-uniform soil water distribution in the rhizosphere (Javaux et al., 2013). Accordingly, a plant can increase water uptake in wetter soil layers in order to compensate for uptake reductions in dryer layers to keep transpiration rate at the potential rate or mitigate transpiration reduction. Several empirical approaches have been developed over the years to account for this so-called compensation mechanism (Jarvis, 1989;Li et al., 2002Li et al., , 2006Lai and Katul, 2000). These models have been incorporated into larger hydrological models and tested at site-specific environments, showing an improvement in predictive quality for e.g. soil water content and crop transpiration (e.g. Braud et al., 2005;Yadav et al., 2009;Dong et al., 2010). Comparisons with physically based models (Jarvis, 2011;de Willigen et al., 2012) that implicitly account for compensation showed that models not including compensation, like Feddes et al. (1978), are less accurate with respect to crop transpiration and soil water content predictions under some circumstances, e.g. at a high root length density.
Recently, De Jong van Lier et al. (2013) developed a mechanistic model for predicting water potentials along the soilroot-leaf pathway, allowing the prediction of RWU and crop transpiration. This model was incorporated into the SWAP eco-hydrological model  by employing a piecewise function between leaf pressure head and relative transpiration, reducing the number of empirical parameters needed when compared to other relations (e.g. Fisher et al., 1981). Besides parameters describing soil hydraulic properties and root geometry, this new model requires information about root radial hydraulic conductivity, xylem axial conductance, and a limiting leaf water potential. Although conceptually interesting, the difficulty in obtaining the required input parameters makes the model less attractive for routine applications.
Other physical RWU models include Couvreur et al. (2012), comparable to De Jong van Lier et al. (2013), as well as more complex three-dimensional models (e.g. Javaux et al., 2013), which account for the full root architecture, requiring more input parameters and a higher computational effort. Specifically, the De Jong van Lier et al. (2013) model differs from the previously mentioned models by the fact that the RWU is based on matric flux potential with an equation derived from the microscopic RWU approach (De Jong van Lier et al., 2008), whereas in other models RWU is based on water pressure head. Scenarios in-cluding an osmotic potential can also be simulated by the model .
Empirical RWU models are more readily used because of their relative simplicity and lower data requirements. On the other hand, their empirical parameters do not have a clear physical meaning and cannot be independently measured. Their limitations under varying environmental conditions are not well established. For the case of the Feddes et al. (1978) transpiration reduction function, threshold values are available in the literature (Taylor and Ashcroft, 1972;Doorenbos and Kassam, 1986) for some crops and some levels of transpiration demand. Nevertheless, experimental (Denmead and Shaw, 1962;Zur et al., 1982) and theoretical (Gardner, 1960;De Jong Van Lier et al., 2006) studies indicate that these parameters should not depend only on crop type and atmospheric demand, but are also determined by root system parameters and soil hydraulic properties. Furthermore, threshold values are hardly ever validated, and they cannot be used for other models (like the Jarvis (1989) model) due to conceptual differences. Accurate values for crops accounting for more environmental factors are necessary to apply these models to a wider range of scenarios. Due to the great number of models developed over the years, it is paramount to investigate some of these before attempting to determine their parameters.
The general purpose of this study was to evaluate the ability of some empirical models to mimic the dynamics of RWU distribution under varying environmental conditions performed in numerical experiments with a detailed physical model proposed by De Jong van Lier et al. (2013). This detailed physical model accounts for hydraulic resistances from the soil to the leaf. We first review some empirical RWU models that have been employed in ecohydrological models and suggest some alternatives. By determining the parameters of the empirical models using inverse modelling of simulated depth-dependent RWU, it becomes clear to what extent the empirical models can mimic the dynamic patterns of RWU.

Theory
RWU and crop transpiration are linked through the principle of mass conservation for water flow in the soil-plantatmosphere pathway: where T a (L T −1 ) is the crop transpiration and S (L 3 L −3 T −1 ) is the root water uptake, dependent on crop properties and soil hydraulic conditions, a function of soil depth z (L), and z m (L) the maximum rooting depth. Equation (1) neglects the change in water storage in the plant, which is justified for daily scale predictions, assuming that α m Feddes et al. (1978) root water uptake reduction function. h 2 and h 3 are the threshold parameters for reduction in root water uptake due to oxygen deficit and water deficit, respectively. The subscripts l and h stand for low and high potential transpiration T p . h 1 and h 4 are the soil pressure head values above and below which root water uptake is zero due to oxygen and water deficit, respectively. (b) Root water uptake reduction function α m as a function of matric flux potential M; M ch and M cl are the critical values of M for high and low T p , respectively, below which the uptake is reduced and M max is the maximum value of M, dependent on soil type.
plants rehydrate to the same early morning water potentials on successive days (Taylor and Klepper, 1978). In a macroscopic modelling approach, RWU is calculated as a sink term S in the Richards equation, which for the vertical coordinate is given by where θ (L 3 L −3 ) is the soil water content, h (L) the soil water pressure head, K (L T −1 ) the soil hydraulic conductivity, t (T) the time, and z (L) the vertical coordinate (positive upward). To apply Eq.
(2), a functional expression for S is needed. Physical equations in analogy to Ohm's law have been suggested (see the review of Molz, 1981, for examples) as well as expressions derived by upscaling microscopic models (De Willigen and van Noordwijk, 1987;Feddes and Raats, 2004;De Jong van Lier et al., 2008. Alternatively, simple empirical models requiring less information about plant and soil hydraulic properties have also been proposed and are commonly used. Most of these models use the Feddes approach (Feddes et al., 1976(Feddes et al., , 1978, formulated as where α(h) is the RWU reduction function, defined by Feddes et al. (1978) as a piecewise linear function of h (Fig. 1). According to this approach, a reduction in S due to α(h[z]) < 1 directly implies a transpiration reduction, causing α(h) to be called a transpiration reduction function. S p is the potential RWU, which is determined by partitioning potential transpiration T p over depth. Several ways to estimate S p as a function of depth have been proposed (Prasad, 1988;Li et al., 2001Li et al., , 2006Raats, 1974), the most common being to distribute T p over depth according to the normalized root length density β (L −1 ) defined as a fraction of root length density R (L L −3 ): Different expressions for α have been suggested, normally considering α a function of θ (e.g. Lai and Katul, 2000;Jarvis, 1989), of h (e.g. Feddes et al., 1978), or of a combination of both (Li et al., 2006). Compared to θ , h seems to be more feasible because of its relation to soil water energy and the fact that obtained parameters of such a function would be more likely applicable to different soils. Some reduction functions, generally associated with reservoir models for soil water balance, correlate RWU with the effective saturation. Regarding the shape of the reduction curve, smooth non-linear functions constrained between wilting point and saturation have been used, as well as piecewise linear functions, but they are all described by at least two empirical parameters. The parameters of the smooth non-linear functions allow easy curve fitting, whereas in the piecewise functions they stand for the threshold at which RWU (or crop transpiration) is reduced due to drought stress, which has been an important parameter in crop water management. Metselaar and De Jong van Lier (2007) showed for a vertically homogeneous root system that the shape of α is not linearly related to soil water content or to pressure head. A linear relation between α and matric flux potential, a composite soil hydraulic function defined in Eq. (5), is physically more plausible and was experimentally corroborated by Casaroli et al. (2010). Matric flux potential is defined as where h w is the soil pressure head at the wilting point. Accordingly, a more suitable expression for α would be a piecewise linear function of M (Fig. 1). RWU can then be calculated by the Feddes model (Eq. 3), replacing its reduction function at the dry side with the alternative illustrated in Fig. 1 where M s is the bulk soil matric flux potential, M 0 the value of M at the root surface, and ρ(z) (L −2 ) a composite parameter, depending on R and root radius r 0 : where r m (= √ 1/π R) (L) is the rhizosphere radius -defined as the half distance between neighbouring roots -and a the distance relative to r m -r 0 where water content equals the average soil water content. In De Jong van Lier et al. (2013), this model is extended by taking into account the hydraulic resistances to water flow within the plant. Dividing water transport within the plant into two physical domains (from root surface to root xylem to leaf), assuming no water changes within the plant tissue and by coupling Eq. (6) for water flow within the rhizosphere, they derived the following expression relating water potentials and T a : where h 0 and h l (L) are the pressure heads at the root surface and leaf, respectively, L l (T −1 ) is the overall conductance of the root xylem-to-leaf pathway, and ϕ (T L −1 ) is defined as where K root (L T −1 ) is the radial root tissue conductivity (referring to the pathway from the root surface to the root xylem), and r x (L) is the xylem radius. An analytical solution of Eq. (8) for h 0 or M 0 depends on an expression for M 0 (h 0 ). For a particular case of Brooks and Corey (1964) where T r (= T a /T p ) is the relative crop transpiration. Because T a and h l are unknowns and T a is undefined when h l = h wl , the equation system cannot be solved analytically. An iterative solution was provided in De Jong van Lier et al. (2013) by defining a maximum transpiration rate T p,max , corresponding to T a (Eq. 8) for h l = h wl . The system of equations is then solved by defining plant stress in terms of T p,max , according to the following boundary conditions: unstressed conditions : T p,max ≥ T p : T a = T p , h l ≥ h wl , stressed conditions : T p,max < T p : h l = h wl , T a < T p .
In the De Jong van Lier et al. (2013) model, crop water stress, a condition for which T a < T p , is defined at the crop level (Tardieu, 1996) and begins when h l = h w . S can be calculated using Eq. (6) by solving Eq. (8), with h 0 (so M 0 ) variable over the root zone and controlled by plant hydraulic properties and soil hydraulic conditions. Figure 2 shows RWU for several values of h l and helps to understand how RWU is distributed over depth. h l can be regarded as a measure of water deficit stress over the whole root zone at crop level: as soil water is depleted, h l is reduced, thus increasing the driving force for RWU. As soil pressure head h s decreases, high uptakes are only achieved by lower h l . For a given h l value, RWU is substantially reduced as h s decreases. If h l is not reduced while h s decreases, S becomes negative (although not shown in Fig. 2, negative S is part of the extension of each curve) and water will flow from root to soil, a phenomenon called hydraulic lift or hydraulic redistribution (Jarvis, 2011). This situation occurs when parts of the root zone are wetter and RWU from these parts satisfies transpiration demand, hence h l is not reduced. Figure 2 also shows that RWU is sensitive to both R and h s and that it can be locally balanced by R and soil water content. Under homogeneous soil water distribution, RWU is partitioned proportionally to R. For heterogeneous conditions, RWU for lower R and higher R may be the same, depending on the stress level (indicated by h l ) and the h s (see Fig. 2). This is in agreement with experimental results reported by several authors (Arya et al., 1975a, b;Green and Clothier, 1995;Verma et al., 2014), who found less densely rooted but wetter parts of the root zone to correspond to a significant portion of RWU when more densely rooted parts of the soil were drier, allowing the crop to maintain transpiration at potential rates. Empirical model concepts that only use R for predicting RWU distribution over depth (under non-stressed conditions) are most common, and therefore these results have been interpreted as being due to a mechanism labelled "compensation" by which uptake is "increased" from wetter layers to compensate for the "reduction" in the drier layers (Jarvis, 1989;Šimůnek and Hopmans, 2009). It is clear, however, that this compensation concept is found merely on a reference RWU distribution based on R, and it only needs to be explicitly addressed in empirical models. In physical models, distinguishing compensation is not necessary since in such models "compensation" follows implicitly from the RWU mechanism.
In order to account for RWU pattern changes due to heterogeneous soil water distribution (the so-called "compensation"), several empirical models have been developed over the years. These models follow the general framework of the Feddes et al. (1978) model given by Eq. (3). Below we review these models and present a new empirical alternative.
2.2 Empirical root water uptake models accounting for compensation 2.2.1 The Jarvis (1989) model Jarvis (1989) defined a weighted-stress index ω (0 ≤ ω ≤ 1) as where, differently from Feddes et al. (1978), α was defined as a function of the effective saturation. Whereas Feddes et al. (1978) assume the RWU reduction directly to be reflected in crop transpiration reduction, the Jarvis (1989) approach employs a so-called "whole-plant stress function" given by where ω c is a threshold value of ω for the transpiration reduction. Substituting Eqs.
(3) and (4) into Eq. (1) (the mass conservation principle) and combining with Eq. (12) results in where α 2 is called the compensation factor of RWU, distinct from Feddes' α (Eq. 3), and which can be derived by defining T a by Eq. (12). In the Jarvis (1989) model, α accounts for local reduction of RWU and transpiration reduction is computed by Eq. (12). When ω = 1, there is no RWU reduction (α = 1 throughout the root zone) and the model prediction is equal to the Feddes model. For ω c < ω < 1, uptake is reduced in some parts of the root zone (as computed by α < 1), but the plant can still achieve potential transpiration rates by increasing RWU over the whole root zone by the factor α 2 .
When ω < ω c , even though the uptake is increased by the factor α 2 , the potential transpiration rate cannot be met. The threshold value ω c places a limit on the plant's ability to deal with soil water stress. When ω c tends to 0, relative transpiration calculated by Eq. (12) tends to 1, and the plant can fully compensate for uptake and transpire at the potential rate if α > 0 at some position within the root zone. In principle, any definition of α is applicable in Eq. 11, and commonly the Feddes et al. (1978) reduction function is used instead of the original Jarvis (1989) reduction function, e.g. in the HYDRUS model (Simunek et al., 2009). This modified version of the Jarvis (1989) model, hereafter referred to as JMf, will be further analysed. Nevertheless, one should be careful in setting up and interpreting the threshold parameters of JMf. The Feddes et al. (1978) model does not account for compensation, and the threshold pressure head value below which RWU is reduced (h 3 ) also represents the value below which transpiration is reduced, making h 3 values from the literature refer to this interpretation. Instead, in the JMf, the transpiration reduction only takes place when ω < ω c , and the soil pressure head in some layers is already supposed to be more negative than h 3 . Therefore, h 3 in JMf is less negative than its namesake in the Feddes model. In that sense, h 3 for the JMf is hard to determine experimentally. An option to do so would be by inverse modelling, optimizing outcomes of soil water flow models with experimental data.

Comparison to the De Jong van Lier et al. (2008) model
The physical basis of Jarvis (1989), defined by Eqs. (11) to (13) with using any α, has been questioned (Skaggs et al., 2006;Javaux et al., 2013). However, the Jarvis model has, to some extent, a physical basis, and a comparison with the physically based model of De Jong van Lier et al. (2008) can be made, as demonstrated in Jarvis (2010Jarvis ( , 2011. This is discussed in the following. De where T p,max (L T −1 ), differently from the definition in the De Jong van Lier et al. (2013) model, is the maximum transpiration rate reached when the root surface pressure head is constant over depth and equal to a limiting value h w . For such a condition M 0 = 0 and T p,max is given by From Eq. (14) we see that drought stress occurs when T p,max < T p . At the onset of drought stress, T a = T p,max . Under this condition, M 0 = 0 and S(z) becomes When T p,max > T p , T a = T p (no drought stress), and M 0 (> 0) is given by Jarvis (2011) observed the similarities between Eqs. (14) and (12) of the models, as well as the algebraic similarity between ω (Eq. 11) and T p,max (Eq. 15). Thus, Jarvis (2010) showed that both models provide the same results under drought stress if α and β(z) are defined as follows: where M max is the maximum value of M (i.e. at h = 0). By substituting Eq. (18) and (19) into Eq. (15) and comparing Eq. (12) with Eq. (14), ω c is found to be equal to Substitution of Eqs. (18) to (20) into Eqs. (12) and (11) Equation (21) is different from Eq. (6) and, therefore, the models cannot be correlated for these conditions. The Jarvis (1989) model predicts RWU by a weighting factor between ρ and M throughout rooting depth. Defining α and β by Eqs. (18) and (19), respectively, allowed us to correlate both models only for stressed conditions. These definitions and the resulting model will be further analysed. Li et al. (2001) proposed to distribute potential transpiration over the root zone by a weighted stress index ζ , being a function of both root distribution and soil water availability:

The Li et al. (2001) model
where α (-) and R (L L −3 ) were previously defined and the exponent l m is an empirical factor modifying the shape of RWU distribution over depth. Originally, the l m values were based on experimental works. For 0 < l m < 1, the RWU in sparsely rooted soil layers is increased in the attempt to mimic compensation. For l m > 1, which has no maximum, the uptake from more densely rooted soil layers increases. Thus, S p is given by and RWU is calculated by substituting Eq. (23) into Eq. (3), following the Feddes approach.
As an alternative to the Jarvis (1989) model, S p can be defined as function of root length density and soil water availability distribution. Compensation is directly accounted for by the weighted stress index in Eq. (22). However, using α to represent soil water availability in Eq. (22) does not mimic properly the compensation mechanism. Compensation may take place before transpiration reduction. Using α in Eq. (22) means that compensation will only take place after the onset of transpiration reduction when α in one or more layers is smaller than 1. The l m parameter may also be interpreted as to account for compensation under non-stressed condition. However, compensation as well as the shape of the RWU distribution are likely to change when a soil becomes drier, and a constant l m cannot account for that. Molz and Remson (1970) and Selim and Iskandar (1978) models

The
Decades before Li et al. (2001), Molz and Remson (1970) and Selim and Iskandar (1978) already suggested to distribute potential transpiration over depth according to root length density and soil water availability. Instead of using α to account for soil water availability, they used soil hydraulic functions. The weighted stress index was defined as where is a soil hydraulic function to account for water availability. Molz and Remson (1970) used soil water diffusivity D (L 2 T −1 ), and Selim and Iskandar (1978) used soil hydraulic conductivity K (L T −1 ) for in Eq. (24). RWU is then calculated by substituting Eq. (24) into Eq. (23) and then into Eq. (3) following the Feddes approach.
These models may better represent RWU and compensation than the Li et al. (2001) model. The compensation is implicitly accounted for by means of in ζ . Since decreases as soil water is depleted, in a heterogeneous soil water distribution ζ in wetter layers is relatively increased because the overall Rdz is reduced due to the reduction of in drier, more densely rooted soil layers. Differently from the Li et al. (2001) model, this change in RWU distribution can occur before the onset of transpiration reduction. Heinen (2014) compared different types of in Eq. (24) such as the relative hydraulic conductivity (K r = K/K sat ) and relative matric flux potential (M r = M/M max ), among others. He found large differences in predicted RWU patterns for different forms of , but did not indicate a preference for a specific one.

Proposed empirical model
In describing soil water availability, the matric flux potential M may be a better choice than K or D, since it integrates K and h or D and θ (Raats, 1974;De Jong van Lier et al., 2013). We propose a new weighted stress index, defined as The exponent l m provides additional flexibility on distributing T P over depth, as also shown by Li et al. (2001). The proposed model differs from Li et al. (2006) only on the hydraulic property to account for soil water availability. The α function used in Li et al. (2006) can only alter RWU distribution after the onset of transpiration reduction, as commented earlier. Contrastingly, M affects RWU distribution before transpiration reduction, integrating the effect of both K and h. The RWU can then be obtained by inserting Eq. (25) into Eq. (23) (S p ) and multiplying by any reduction function, such as the Feddes et al. (1978) and proposed reduction functions. In other words, the model follows the Feddes approach, which computes RWU by the two mentioned steps, which differ only with respect to the way S p is obtained: Eq. (25) (multiplied by T p ) versus Eq. (4).

Material and methods
3.1 Applied models Table 1 summarizes the empirical RWU models evaluated in this study. They all follow the original Feddes model (Eq. 3), but differ in how RWU is partitioned over rooting depth or how α is defined. For each model, except for Jarvis (2010), we defined a modified version by substituting the Feddes reduction function by the proposed reduction function (Fig. 1b), and these modified versions were also evaluated. The threshold values of the Feddes et al. (1978) reduction function for anoxic conditions (h 1 and h 2 ) were set to zero. The value of the parameter h 4 was set to −150 m. The other parameters of the models were obtained by optimization as described in Sect. 3.2.
All these models were embedded as sub-models in the SWAP ecohydrological model , allowing us to solve Eq. (2) and to apply it to different scenarios of root length density, atmospheric demand, and soil type (described in Sect. 3.2) to analyse the behaviour and sensitivity of the models. Simulation results of SWAP in combination with each of the RWU models were compared to the SWAP predictions in combination with the physical RWU

Drying-out simulation
Boundary conditions for drying-out simulations were no rain/irrigation and a constant atmospheric demand (potential transpiration) over time. The simulation continued until simulated crop transpiration by the physical RWU model approached zero. Soil evaporation was set to zero, making soil water depleted only due to RWU or bottom drainage. The free drainage (unit hydraulic gradient) at the maximum rooting depth was the bottom boundary condition. The soil was initially at hydrostatic equilibrium with a water table located at 1 m depth. We performed simulations for two levels of atmospheric demand given by potential transpiration (T p ) of 1 and 5 mm d −1 . We also considered three soil types and three levels of root length density, as described in the following.

Soil type
Soil data for three top soils from the Dutch Staring series (Wösten et al., 1999) were used. The physical properties of these soils are described by the Mualem-van Genuchten functions (Mualem, 1976); (Van Genuchten, 1980) for the K − θ − h relations: where = (θ − θ r )/(θ s − θ s ); θ, θ r and θ s are water content, residual water content and saturated water content (L 3 L −3 ), respectively; h is pressure head (L); K and K sat are hydraulic conductivity and saturated hydraulic conductivity, respectively (L T −1 ); and α (L −1 ), n, and λ are empirical parameters. The parameter values for the three soils are listed Table 1. Summary of empirical models used in this study. α f and α m are the Feddes et al. (1978) (Fig. 1a) and proposed reduction functions (Fig. 1b), S p (Eq. 4) is the potential root water uptake, ω (Eq. 11) and ω c are the weighted stress index and threshold value in the Jarvis (1989) model, and ζ m (Eq. 25) is the weighted stress index in the proposed models.

Model
Acronym Equation Feddes et al. (1978) (2010) Table 3. These soils are identified in this text as clay, loam and sand.

Root length density distribution
Three levels of root length density were used, according to the range of values normally found in the literature. We considered low, medium and high root length density for average crop values equal to 0.01, 0.1 and 1.0 cm cm −3 , respectively. For all cases, we set the maximum rooting depth z max equal to 0.5 m. Root length density over depth z was described by the exponential function: where R 0 (L L −3 ) is the root length density at the soil surface, b (-) is a shape-factor parameter, and z r (= z/z max ) is the relative soil root depth. The term (1−z r ) in Eq. (28) guarantees that root length density is zero at the maximum rooting depth. The parameter R 0 is hardly ever determined, whereas the average root length density of crops R avg is usually reported in the literature. Assuming R of such a crop given by Eq. (28), it can be shown that  Solving Eq. (29) for R 0 and substituting into Eq. (28) gives Figure 3 shows R(z r ) calculated from Eq. (30) for different values of b and R avg = 1 cm cm −3 . As b approaches zero, Eq. (30) tends to become linear; however, it is not defined for b = 0. In our simulations b was arbitrarily set equal to 2.0.

Optimization
The parameters of the empirical RWU models were estimated by solving the following constrained optimization Table 3. Mualem-van Genuchten parameters for three soils of the Dutch Staring series (Wösten et al., 1999) used in simulations. θ s and θ r are the saturated and residual water content, respectively; K s is saturated hydraulic conductivity and α, λ, and n are fitting parameters.  Table 4. Parameters of the root water uptake models estimated by optimization and their respective constraints .

Model Parameter Unit
where ( Table 1. p is the model parameter vector to be optimized, constrained in the domain . Both p and vary depending on the empirical RWU model used. Table 4 shows the parameters of each empirical RWU model that were optimized and their respective constraints . m is the number of soil layers (50 soil layers of 1 cm thickness) and n is the duration, in days, of the simulation. The Jarvis (2010) model has no empirical parameters and therefore requires no optimization. Equation (31) was solved using the PEST (Parameter ES-Timation) tool (Doherty et al., 2005) coupled to the adapted version of SWAP. PEST is a non-linear parameter estimation program that solves Eq. (31) by the Gauss-Levenberg-Marquardt (GLM) algorithm, searching for the deviation, initially along the steepest gradient of the objective function and switching the search gradually to the Gauss-Newton algorithm as the minimum of the objective function is approached. Upon setting PEST parameters, we made reference runs of SWAP with each empirical model using random values of p aiming to assess the ability of PEST to retrieve p. These reference runs allowed us to properly set up PEST for our case. For highly non-linear problems as in Eq. (31), the optimized parameter set depends on the initial values of b. We used five random sets of initial values for p in order to guarantee that GLM encountered the global minimum and also to check the uniqueness of the solution. Runs led to the same minimum in most cases, but if not, the minimum was compared and a fit run was made again.
The optimizations were performed for the drying-out simulation only. This guaranteed that RWU predictions from SWAP corresponded to the best fit of each empirical model to the De Jong van Lier et al. (2013) model. This analysis aimed to investigate the capacity of the empirical RWU models to mimic the RWU pattern predicted by the De Jong van Lier et al. (2013) model. The optimized parameters were subsequently used to evaluate the models in an independent growing season scenario.

Growing season simulations
In the growing season simulation, all models were evaluated by simulating the transpiration of grass with weather data from the De Bilt weather station, the Netherlands (52 • 06 N; 5 • 11 E), for the year 2006. The same root system distribution as in the drying-out simulations was used, i.e. a crop with roots exponentially distributed over depth as Eq. (30) (b = 2.0) down to 50 cm below the soil surface. We also performed simulations for the same three types of soils and root length densities. In all cases the crop fully covered the soil with a leaf area index of 3.0. Daily reference evapotranspiration ET 0 was calculated by SWAP using the FAO Penman-Monteith method (Allen et al., 1998). In the SWAP model, a potential crop evapotranspiration ET p is obtained by multiplying ET 0 by a crop factor, which for the grass vegetation was set to 1 . ET p was partitioned over potential evaporation E p and T p using parameter values for common crops given in the SWAP model (see Van Dam et al., 2008, for details).
The values of the empirical parameters of each RWU model corresponding to the type of soil and root length density were taken from the optimizations performed in the drying-out experiment. Each parameter was estimated for two levels of T p (1 and 5 mm d −1 ) and was linearly inter-  polated for intermediate levels of T p . For T p higher than 5 mm d −1 and T p lower than 1 mm d −1 , the values estimated for 5 and 1 mm d −1 , respectively, were used. As in the drying-out simulations, the bottom boundary condition was free drainage. Initial pressure heads were obtained by iteratively running SWAP starting with the final pressure heads of the previous simulation until convergence.  Figure 4 shows the RWU patterns for the case of the clay soil for the three evaluated root length densities R and the two levels of potential transpiration T p . It can be seen how R and T p affect RWU distribution and transpiration reduction as the soil dries out. The onset and shape of transpiration reduction is affected by the RWU pattern. For low R, the low number of roots in deeper layers is not sufficient to supply high RWU rates. When the upper layers become drier, transpiration reduction follows immediately. Under medium and high R, the RWU front moves gradually downward as water from the upper layers is depleted. Comparing from high to medium R, the RWU front goes even deeper, and transpiration is maintained at potential rates for a longer time (Fig. 4). Accordingly, the plant exploits the whole root zone and little water is left when transpiration reduction onsets, causing an abrupt drop in transpiration. Regarding T p , the RWU patterns are very similar for both evaluated rates, differing only in time scale: for high T p the onset of transpiration reduction and the shift in RWU front occur earlier. The uptake patterns for the sand and loam soil (not shown here) are very similar. However, for the sand soil potential transpiration is maintained a little longer and more water is extracted from deeper layers. For the loam soil, the onset of transpiration reduction occurred earlier.
The leaf pressure head h l over time shown in Fig. 4 illustrates how the model adapts h l to R and T p levels in a drying soil. Initially all scenarios have the same water content distribution and lower h l values are required for low R or high T p scenarios to supply potential transpiration rates. As soil becomes drier, h l is decreased to increase the pressure head gradient between bulk soil and root surface, thus maintaining RWU corresponding to the demand. Therefore, uptake in wetter layers becomes more important. Transpiration reduction only onsets when h l reaches the limiting leaf pressure head h wl (= −200 m), after significant changes in the RWU patterns, characterized by increased uptake from deeper layers.
For the high T p -low R scenarios, transpiration reduction starts on the first day of simulation, although the soil is relatively wet. This is a case of transpiration reduction under non-limiting soil hydraulic conditions due to high atmospheric demand (Cowan, 1965). For such conditions, the high water flow within the plant required to meet the atmospheric demand cannot be supported by the root system with a low R and hydraulic parameters given in Table 2. Higher atmospheric demand (here represented by T p ) leads to faster reduction of h l caused by the hydraulic resistance to water flow within the plant, and the transpiration rate and RWU are a function of h l . The physical model assumes a parsimonious relationship (Eq. 10) between transpiration and h l : transpiration rate is only reduced when h l reaches a limiting value h wl , which corresponds to a maximum transpiration rate T p,max allowed by the plant for the current soil hydraulic and atmospheric conditions. Under non-limiting soil hydraulic conditions, root system properties and plant hydraulic parameters (Table 2) are the major determining factors for T p,max , whereas soil hydraulic conditions play a minor role. Figure 5 shows T p,max as a function of K root for some values of L l with a constant soil pressure head of −1 m in the root zone for low R in the sandy soil. In this scenario, K root limits the crop transpiration and L l becomes important only when K root increases. The potential transpiration can be achieved by rais-  ing K root to about 10 −7 m d −1 . This can also be achieved by decreasing h wl (not shown in Fig. 5).
In the field, transpiration rate and root length density are related to each other: a high transpiration rate only occurs in a high leaf area, and a high leaf area implies a high root length density. Thus, even under very dry and hot weather conditions, a crop with a low R may not be able to realize a high transpiration rate. Furthermore, crop transpiration depends on the stomatal conductance. In the De Jong van Lier et al. (2013) model, this is implicitly taken into account by the simple relationship between h l and T a . However, stomatal conductance is relatively complex and depends on several environmental factors such as air temperature, solar radiation, and CO 2 concentration. Therefore, high potential transpiration rates may not be achieved because of the stomatal conductance reduction due to temperature or solar radiation. This behaviour could be simulated by the cou-

Root water uptake pattern predicted by the empirical models
In this section, we evaluate the empirical RWU models (models and their abbreviations are listed in Table 1) Table 5 and are discussed in Sect. 4.1.4, and therefore represent the best fit with VLM. The RWU patterns simulated by VLM and the empirical models for the sandy soil and high R scenario are shown in Figs. 6 and 7 for low and high T p , respectively. Both versions of the Feddes model (FM and FMm) predicted enhanced RWU from the upper soil layers. When the soil pressure head (h s for FM) or soil matric flux potential [M s for FMm] is greater than the threshold value for uptake reduction, these uptake patterns are equivalent to the vertical R distribution. For conditions drier than the threshold value (when α f and α m are less than 1), the predicted RWU patterns by the models become different (Figs. 6 and 7).
When reducing RWU for a period depending on R, T p , and h 3 , RWU from the upper soil layers predicted by FM rapidly decreases to zero. This zero-uptake zone expands downward as soil dries out. On the other hand, the uptake predicted by FMm is substantially reduced right after the onset of transpiration reduction, proceeding at lower rates and for a much longer time until approaching zero. These features become evident by comparing the shapes of both reduction functions (Fig. 8). α m is linear with M after M > M c , but it is concavely shaped as a function of h -as also shown by Metselaar and De Jong van Lier (2007) and De . This makes α m decrease abruptly for M > M c , causing a substantial decrease in RWU even when h is slightly below the threshold value. Therefore, RWU proceeds at low rates for a longer time. In contrast and due to the linear shape of α f , RWU predicted by FM remains higher for a longer time after h < h 3 . FM does not predict an abrupt change in RWU patterns, especially when T p is low (Fig. 6). When h approaches h 4 , α f is still relatively high and RWU makes h decrease rapidly. Another diverging feature between α f and α m , also shown in Fig. 8, is that the shape of α m varies with soil type (regardless of the value of its threshold parameter M c ), whereas α f does not. These different features of the reduction functions also affect the matching values of the parameters, as discussed below. Although the choice of the reduction function affects transpiration over time only slightly, RWU patterns are strongly affected (Figs. 6 and 7).
The RWU patterns predicted by the JMf and JMm models can be very different, as shown by Fig. 6 for the high R-low T p scenario. In this scenario, the JMf model did not predict any compensation because the optimal ω c equalled 1 (Table 5) -thus becoming identical to FM -and the optimal h 3 values for JMf and FM were similar. In Fig. 6, although h 3 values for FM and JMf (ω c = 1) are close to zero, the plant transpiration is near T p for a prolonged time due to a small reduction of α. These high R-low T p scenarios with a high R in deep soil layers allow RWU at higher rates when surface soil layers become drier (as predicted by VLM). Then, the reduction of ω c , an attempt to numerically predict compensation with JMf, makes the RWU pattern deviate even more from the VLM pattern. This is illustrated in Fig. 6 and by the optimal h 3 and ω c values shown in Table 5. In order to mimic the VLM uptake patterns, the value of h 3 for all soil types in this scenario was equal or close to zero. Decreasing h 3 or ω c to simulate compensation makes JMf predict higher uptake from upper layers, increasing the discrepancy between the models. The optimal ω c for all soil types was equal to 1 (in other words: there was no compensation). RWU in the upper layers predicted by VLM is substantially reduced within a few days, whereas reducing ω c in the JMf model to predict compensation has the side effect of causing an increase in uptake from upper layers. The model, therefore, is not able to adequately mimic the scenarios with compensation evaluated here. On the other hand, the JMm model was able to reproduce considerably well the VLM pattern for the evaluated scenarios due to the shape of α m as discussed above. As soon as M < M c in the upper layers, RWU decreased at a higher rate, compensated for by increasing uptake from the wetter, deeper layers. This agrees more closely with VLM predictions.
For high T p (Fig. 7), the JMf model can predict compensation (ω c < 1); however, its predicted RWU pattern is quite different from JMm and VLM. JMf predicts a higher longerlasting RWU near the soil surface than the other models that account for compensation. This makes soil water depletion more intense and RWU from these layers to cease sooner when h s becomes lower than h 4 . At this point, T a is predicted to continue to be equal to T p because of the low optimal ω c (= 0.19), which increases RWU from the deeper layers where h is close or equal to h 4 . JMm performed very differently, predicting uptake over the first few days (when M s > M c ) in accordance with R distribution. After M < M c in the upper soil layers, the RWU pattern started to change gradually and RWU increased at lower depths.
The proposed models (PM and PMm) are capable of predicting RWU patterns similar to VLM. For the low T p -high R scenario (Fig. 6), RWU is more uniformly distributed over depth than in the VLM model for the first days and uptake from upper layers is lower than that predicted by the VLM model. For high T p (Fig. 7), these models better represent RWU patterns and, in general, differences in predictions of RWU between the proposed models are small. The shape of the transpiration reduction over time, however, is smoother than predicted by the VLM model. Concerning the relative transpiration curve, the proposed models appear to be less precise than the other models that account for RWU compensation.
JMII does not mimic well the RWU pattern predicted by VLM for the high R-low T p scenarios. It overestimates uptake from surface layers during the first days. Before the onset of transpiration reduction, uptake from upper layers reaches zero, but it is compensated for by a higher uptake from deeper layers. The model is very sensitive to both R and M. For the high R-high T p scenarios, JMII provides better uptake pattern predictions (Fig. 7). However, the model does not perform well in the other scenarios with low and medium R (data not shown here). Comparing RWU predictions from JMf and JMII, the Jarvis-type models are affected by the definition of α. This becomes clear from Fig. 9, which shows the α of JMII (Eq. 18) as a function of h s and ω c (Eq. 20) for different soil types, expressed by M max . The α function shows that even though the soil resistance increases as the soil becomes drier, defining α by Eq. (18) does not seem plausible. In this case, α is suddenly reduced when the soil is still near saturation. When h s = 1 m, for instance, α is much lower than 0.5. Such behaviour is not reasonably compatible with for the α concept. The ω c values are also extremely low. The low α values are, however, balanced by high α 2 values (due to low ω and ω c values), leading to suitable values of RWU in a given soil layer. Nevertheless, the magnitudes of α and ω c are conceptually questionable. Therefore, we conclude that (i) the ω c value in Jarvis-type models, which sets the compensation level, depends on the definition of α. For instance, for the original Jarvis (1989) model, ω c = 0.5 corresponds to a moderate level of compensation. Surely, this would not hold if α were defined by Eq. (18); (ii) comparing the Jarvis (1989)  The fact that JMII is more sensitive to both R and M, as stated above, when compared to the other M-based models is attributed to the α function and the derived equations to express their parameters (Eqs. 19 and 20). It can be seen from Fig. 9c that β defined by Eq. (19) (β of JMII) tends to be higher when R increases and lower when R decreases compared to the β of JMf and JMm. Thereby, for the first days of simulations when the soil hydraulic conditions tend to be rather uniform over depth, JMII overestimates RWU compared to VLM predictions. This becomes more important for the high R-low T p scenarios. For such conditions, the RWU over depth predicted by the VLM tends to be more uniform, which seems reasonable as the low transpiration demand can be met by any small R in deeper soil depths. After some time, the discrepancies between VLM and JMII tend to increase, since the higher RWU in the upper layers reduces h; thus, because of the α shape of JMII, RWU in the upper layers is suddenly reduced towards zero. These are the main reasons for JMII not to predict well in the high R-low T p scenarios.  Figure 8. Feddes et al. (1978) (α f , grey lines) and proposed (α m , black lines) water uptake reduction functions as a function of soil pressure head h using their respective optimized parameters for the scenario of high root length density, three types of soil, and two potential transpiration levels.

Statistical indices
The performance of the empirical models was analysed by the coefficient of determination r 2 and the model efficiency coefficient E (Nash and Sutcliffe, 1970) calculated by comparing to the RWU and relative transpiration predicted by VLM. For the low R-high T p scenarios, the VLM predicts water stress (T a < T p ) from the beginning of the simulation as discussed in Sect. 4.1.1. The empirical models (except for JMf and JMm by setting ω c > 1) are not able to reproduce these results, thus these scenarios were not considered when analysing the performance of the models. Statistical indices for the evaluated scenarios of each model are concisely shown by the boxplots in Fig. 10. The width of whiskers indicates the range of the statistical indices for each model used in the evaluated scenarios. The outliers indicate whether a model had different performance at some scenarios than its overall performance. Focusing first on RWU, the figure shows that the proposed models performed better. The performance of PM was just a bit poorer  than PMm's, shown by the presence of an outlier and lower median. JMm performed as good as the proposed models, and only in two scenarios it had a bad performance as shown by the outliers in Fig. 10. The wider whiskers and presence Among the models that account for RWU compensation, JMf and JMII performed worst, especially in the high R-low T p scenarios. In general their performances were poorer for medium R scenarios, especially for low T p . Thus, the use of α m in Jarvis-type models promotes substantial improvements, especially from medium to high R scenarios. For low R scenarios all models performed well and the highest values of the boxes in Fig. 10 usually refer to this scenario.
In predicting transpiration, all models accounting for compensation performed well, except JMf. It can be noticed that JMII performed much better in predicting transpiration than RWU. As for the RWU, all models performed worse in high R scenarios than in low R scenarios.
As the evaluated models differ regarding the number of empirical parameters (from 0 to 2), it is important to use a statistical measure that accounts for this and penalizes the models with more parameters. The Akaike's information cri-teria (AIC) is a suitable measure for such a model comparison. The selection of the "best" model is determined by an AIC score, defined as (Burnham and Anderson, 2002): where K is the number of fitting parameters and L(θ |y) is the log-likelihood at its maximum point. The "best" model is the one with the lowest AIC score. Table 6 lists the best models for every scenario based on the AIC score. Overall, the AIC supports the above descriptive statistical analyses, indicating that the proposed models are the best models in predicting RWU estimated by VLM, especially from medium to high R scenarios. For the low R scenarios JMm is the best model. On predicting T r by VLM, the above analyses indicated that, in general, most models performed similarly. The AIC indicated comparable results, but overall JMm was the best model. The proposed models (PM or PMm) were the best models for high R-low T p scenarios.

Relation of the optimal empirical parameters to R and T p levels
The optimal values of the empirical parameters of all models (except JMII that has no empirical parameters) for all scenarios but the high T p -low R scenario are shown in Table 5. The threshold reduction transpiration parameters h 3 and M c (for FM and FMm, respectively) stand for the soil hydraulic conditions at which the crop cannot meet its potential transpiration rate. Conceptually, the higher R, the lower is h 3 or M c due to the larger root surface area for RWU, i.e. the crop can extract water in drier soil conditions. Similarly, lower h 3 and M c are expected for low T p . This can also be deduced from Figs. 6 and 7 by means of the predictions of relative transpiration and RWU by VLM. The optimal h 3 and M c values (Table 5) for FM and FMm, respectively, increase with R, contradicting their conceptual relation to R. For T p , there is no specific relationship for these parameters: whether they increase or decrease with T p depends on the value of R. In drying-out scenarios, soil water from top layers depletes rapidly due to the higher initial uptake. Thus, uptake from these layer starts to decrease, whereas RWU in deeper, wetter layers increases. This effect becomes stronger at higher R, as seen by the VLM predictions in Sect. 4.1.1. Because FM and FMm do not account for this mechanism, decreasing h 3 or M c in search of conceptually meaningful values would make these models predict higher RWU from upper layers (in accordance with the R distribution) for a longer period, increasing the discrepancy with VLM predictions. Therefore, their best fitted values are physically without meaning due to the model assumptions.
In order to interpret the parameters in Table 5 for JMf, one should first recall that α in JMf stands for the local RWU reduction due to soil hydraulic resistance. Thus, its h 3 parameter refers to the local soil pressure head at which RWU starts to decrease. It may be argued that RWU reduction occurs in drier soil conditions as R increases, i.e. h 3 is more negative for higher R (similarly as for FM and FMm). However, since JMf accounts for compensation, RWU is interpreted as a non-local process, and uptake from one layer depends on the water status and root properties from other layers (Javaux et al., 2013). Thus, the h 3 parameter from JM is affected by other parts of the root zone. Predictions by VLM show that RWU reduction from the upper layers starts at less negative pressure head values as R increases. Therefore, h 3 in JMf should increase with increasing R. The values of h 3 for JMf shown in Table 5 agree with this conceptual meaning. The M c parameter from JMm can be interpreted likewise.
Values for ω c from JMf for the high R-low T p scenarios equal 1, thus contradicting its conceptual meaning: as in these scenarios the compensation mechanism is more intense, ω c should be less than one for the medium and high R scenarios. The reason for ω c = 1 was discussed in Sect. 4.1.2. Conversely, ω c values for JMm follow the conceptual meaning.
The optimal parameters of the proposed models follow their logical relation to R and T p . The l m values for both models are very close. The optimal l m values are less sensitive to soil types and more sensitive to R.
High correlation parameters might result in uncertainties and a non-unique solution of the optimization problem. In general, the correlation parameter coefficients were low, except for some scenarios in which high correlation coefficients between ω c and h 3 (or M c ) were found. These high correlations may be due to model structure rather than to the data used for fitting the models, since the correlations for PM and PMm parameters were low (absolute correlation coefficient below 0.53).

Optimization using T r
The empirical models fitted only to RWU, since the primary interest is to evaluate the model's capability to predict the RWU patterns under different scenarios. RWU is not easily obtained in real conditions, making the use of physical RWU models a great advantage. On the other hand, plant transpiration, one of the main outputs in RWU models, is more easily measured. Thus, one might consider to fit the models to the temporal course of (relative) plant transpiration or to fit the models simultaneously to both plant transpiration and RWU, for which a rather complicated optimization scheme would be required.
We addressed this issue by fitting the models to the course of relative transpiration for some scenarios. The procedure was the same as explained in Sect. 3.2, but substituting S i,j in Eq. (31) by T r i . The results for some models in two contrasting scenarios of R are shown in Fig 11. Models that account for "compensation" can predict T r quite reasonably even when fitted to RWU only. The models that do not account for "compensation" do not mimic the T r course over time correctly for the high R scenario predicted by VLM, even when they are fitted to T r , and the predictive quality decreases when fitted to RWU. The most important aspect shown in Fig 11 is that fitting the models to T r can improve T r predictions but impairs their RWU predictions considerably, especially in high R scenarios. Conversely, if a model fits well to RWU, it can provide suitable transpiration predictions. This can also be seen by the analysis of Sect. 4.1.3, when the proposed models and JMm had good performance in predicting T r as well.

Growing season simulation
By evaluating the RWU models under real weather conditions during a relatively dry year and considering the same soil types and crop characteristics as for the drying-out experiment, it was possible to use the calibrated parameters for specific soil type and root length density. This evaluation is important to analyse whether our calibration of the empirical models with a single drying-out experiment results in consistent predictions for other circumstances. Models were not evaluated for the low R scenario because the empirical models (except JMf and JMm) were not able to mimic those conditions for high T p (Sect. 4.1.1). Figure 12 shows the time course of cumulative actual transpiration simulated by SWAP using all the RWU models, together with rainfall and T p throughout the growing season period. Following the first dry spell, T ac predicted by FM and FMm, not accounting for "compensation", starts to be lower than predictions from other models. Two or three more dry spells occur in the evaluated period. The magnitude of the underestimation, however, varies with soil type and R. For the medium R-loam soil scenario, for instance, the T ac for all models are similar. The T ac at the end of the evaluated period predicted by VLM for low R (not shown in Fig. 12) was much lower and approximately equal for the three soil types (40.45, 40.05 and 40.08 cm for clay, loam and sand soil, respectively). In fact, a higher R resulted in an increasing difference of cumulative transpiration between soil types. Most water is extracted from the clay soil, followed by sand and loam. Little difference of cumulative transpiration is found between medium and high R: for sand and clay soil, the cumulative transpiration was slightly higher for high R; for the loam soil it was and practically identical.
Comparing cumulative T a predicted by the empirical models with VLM predictions shows that the models that do not account for compensation underestimate cumulative T a from 2.0 % (medium R -sand soil scenario) to 13.9 % (high R -  clay soil scenario). Overall, the highest underestimates occurred for high R. All other models predict similar values. Therefore, for total actual transpiration prediction, any of the evaluated models accounting for compensation might be suitable after calibration.
An overall analysis of model performance is shown in Fig. 13 and a list of the "best" model for each scenario based on AIC is shown in Table 7. The best performances are from the models that account for compensation. An improvement of JMf by using the proposed reduction function can be observed. Among the models that account for compensation, JMf had the worst performance. JMII also was poor in predicting RWU, but showed good performance in estimating plant transpiration. Overall, the best performances were also obtained by the proposed models (PM and PMm) and by the modified Jarvis (1989) model (JMm) in predicting RWU. These results also indicate that the strategy of designing a single drying-out experiment to calibrate an empirical model is successful.
According to the AIC, PM, PMm and JMm are best in predicting RWU. Regarding T r predictions, Fig. 13 shows considerably high statistical indices (E and r 2 ) for all models that account for "compensation". However, the AIC, which penalizes the models with more parameters, indicates that JMII was the "best" model for most of the scenarios.
In general, the proposed models as well as JMm showed better performance than the other empirical models. It should be noted, however, that these models are based on M, making them closer to the physical De Jong van Lier et al. (2013) model. In this regard, it is important to separately compare JMf and JMm and PM and PMm. The only difference be-tween JMf and JMm is the α reduction, which resulted in considerable improvements as discussed. In the proposed models, M is included in S p (z) to distribute T p over depth. In PMm, α m is used instead of the Feddes reduction function (used in PM). These simple modifications were sufficient to allow these empirical models to be fitted too mimic the predictions made by the more complex physical model.

Conclusions
Several simple RWU models have been developed over the years, and we outlined some of these models and also proposed alternatives. Some of these models were embedded as sub-models in the SWAP eco-hydrological model    -Regarding the ability of the models in predicting plant transpiration, all models accounting for compensation have good performance. The AIC indicates that JMII is the "best model". This model is also more suitable for blind predictions, as no empirical parameters need to be estimated.
-The simulations of a growing season with grass confirmed these findings, suggesting that an experiment of soil drying-out for two levels of potential transpiration, as performed, is adequate for analysing the performance of RWU models and retrieving their empirical parameters by defining the objective function in terms of RWU.
It should be noticed that the predictions from the De Jong van Lier et al. (2013) physical model do not represent a real system. However, they show to be consistent with observed behaviour and have adequate sensitivity to variables and system boundaries. It is common practice to refer to the parameter compilation made by Taylor and Ashcroft (1972), which does not account for the dependence of the parameters on soil type. Moreover, these parameters depend on type of transpiration reduction function; although not explicit in the Taylor and Ashcroft (1972) compilation, it is usually considered to refer to the Feddes model. The best empirical models in predicting RWU, based on the comparison with the De Jong van Lier et al. (2013) physical model (the proposed models and JMm), contain one additional parameter, also dependent on soil type, root length density and potential transpiration. Although the parameters for three soil types, root length density and potential transpiration are provided in this study, a more robust and complete calibration may be necessary, mainly because general values of plant hydraulic resistances were used. Due to the dependence of the empirical parameters on soil type and potential transpiration, parameterizing the selected empirical models for a specific crop might require more effort than when using the physical model with parameters that can be determined independently. The use of the physical model predictions, as in this study, seems a good strategy to calibrate the empirical models. Ultimately, the option for the empirical or physical model will be based on the desired complexity and understanding of the system, and on the availability of parameter values.
Model, input data and optionally modeling results are available from the corresponding author upon request.