WATER RESOURCES RESEARCH, VOL. 37, NO. 10, PAGES 2457–2470, OCTOBER 2001

One-, two-, and three-dimensional root water uptake functions
for transient modeling
J. A. Vrugt and M. T. van Wijk
Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam, Amsterdam, Netherlands

J. W. Hopmans
Hydrology Program, Department of Land, Air and Water Resources (LAWR), University of California
Davis, California, USA

J. Šimunek
U.S. Department of Agriculture Salinity Laboratory, University of California, Riverside, California, USA

Abstract. Although solutions of multidimensional transient water flow can be obtainedby numerical modeling, their application may be limited as root water uptake is generally
considered to be one- or two-dimensional only. This is especially the case for trees. The
first objective of this paper is to test the suitability of a three-dimensional root water
uptake model for the simultaneous simulation of transient soil water flow around analmond tree. The soil hydraulic and root water uptake parameters were optimized byminimizing the residuals between measured and simulated water content data. Water
content was measured in a three-dimensional grid around a sprinkler-irrigated almondtree for a 16 day period, following irrigation. A second objective was to compare theperformance and results of the three-dimensional flow model with one- and twodimensional root water uptake models. For this purpose, measured water contents were
aggregated in the x and y direction in the one-dimensional case and in the radial direction
for the two-dimensional uptake model. For the estimation of root water uptake model
parameters a genetic algorithm was used to estimate the approximate global minimum of
the parameter space, whereas final parameters were determined using the Simplex
optimization algorithm. With the optimized root water uptake parameters, simulated andmeasured water contents during the 16-day period were in excellent agreement for all root
water uptake models. Most significantly, the spatial variation in flux density below therooting zone decreased when reducing multidimensional root water uptake to fewer
dimensions, thereby justifying the proposed multidimensional approach.

1.

Introduction

From a hydrological perspective, water uptake by root systems and their spatial distribution may exert a large degree ofcontrol on the water fluxes to the atmosphere and the groundwater [Canadell et al., 1996]. For an improved understanding of
the magnitude of these fluxes, accurate estimates of the temporal and spatial root water uptake patterns are needed.
Clearly, quantification of root water extraction rates also contributes to an understanding of chemical fluxes in the vadosezone in both ecological and hydrological studies [Somma et al.,1998] as well as their control by vegetation. Interactions between roots and soil in the rhizosphere influence the quality
and quantity of water transported in and exported from thevadose zone. An understanding of the interactions between theroots and surrounding soil and solutes under a variety ofchanging environmental conditions has large implications since
it will lead to a decrease in contamination of subsurface andsurface water resources by reducing the loss of fertilizers andother agrochemicals below the root zone [Clausnitzer and Hopmans, 1994; Clothier and Green, 1994]. Moreover, the rhizosCopyright 2001 by the American Geophysical Union.
Paper number 2000WR000027.
0043-1397/01/2000WR000027$09.00

phere might be responsible for accelerated breakdown of organic chemicals by biodegradation [Walton and Anderson,1990] and phytoremediation [Nyer and Gatliff, 1996]; hence athorough understanding of root function regarding uptake ofwater and associated solutes is warranted.
Actual root water uptake spatially depends not only on theroot density distribution but also on its temporal functioning as
determined by soil water availability and soil salinity. In addition to water stress in periods of low water availability, root
water uptake is also reduced when concentrations of soluble
salts exceed plant-specific threshold values [Homaee, 1999]. In
irrigated soils, particularly in arid and semiarid regions, plants
are generally subjected to both salinity and water stress. In
these regions, soil and water management practices are basedon maintaining a favorable soil water content and salinity status in the root zone, thereby minimizing periods of water stresswhile controlling leaching to minimize salinity stress.
The influence of plant root systems on water and chemical
movement can be better understood using soil water simulation models, provided that accurate spatial and temporal root
water uptake distributions are included [Musters and Bouten,
1999]. One of the earliest detailed quantitative studies of water
extraction by a plant root was presented by Gardner [1960].
This microscopic model considered a single root to be equiv-

2457

2458

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

alent to an infinitely long cylinder of uniform radius with water-absorbing properties. The steady state soil water flow equation was solved analytically assuming radial flow and imposedroot water uptake rates. Soil water matric head distributions
around the idealized root were calculated. This concept was
extended in later papers [Gardner, 1964, 1965; Gardner andEhlig, 1962] and proved to be very insightful but lacked practical applicability since the detailed geometry of the rooting
system is difficult to measure and is time dependent. Consequently, most root water extraction terms have been developed
using a macroscopic rather than a microscopic approach.
Nonetheless, all macroscopic models preserve the essence of
Gardner’s [1960] insight.
In the macroscopic approach of Richards’ equation a sinkterm representing the water extraction of the entire root system is included to describe transient multidimensional waterflow [Whisler et al., 1968; Molz and Remson, 1970; Clausnitzer
and Hopmans, 1994], according to⭸␪⫽ ⵜ ⴢ 关Kⵜ共h ⫺ z兲兴 ⫺ S共 x, y, z, t兲,⭸t(1)

where ␪ is the volumetric water content (L 3 L ⫺3 ), K is theunsaturated hydraulic conductivity tensor (L T 1 ), h (L) is thesoil water matric head, z (L) is the depth which is included forvertical flow only, and S is the volumetric sink term (L 3 L ⫺3
T ⫺1 ), representing root water uptake as a function of both
space and time. The benefits of such an approach are evident
since it allows direct integration of root water uptake withtransient soil water flow and provides natural interactions between transpiration and root water extraction, that is, as longas we know the mechanisms with which to describe S.
Though transient soil water flow in the vadose zone is often
simulated in one, two, and three spatial dimension; root water
uptake is generally considered simply to be a function of thevertical dimension only. For uniform crops with a spatially
uniform water uptake pattern, one-dimensional models may
suffice. However, for row crops and tree lines, for example, a
two-dimensional representation would be better. For isolated
trees, such as apples and almonds in large monocultures, theprocess of water uptake is complex, and a three-dimensional
representation would therefore seem appropriate [Green andClothier, 1999]. Flow models such as HYDRUS-3D, being an
update of the SWMS_3D code of Šimunek et al. [1995], allow
multidimensional root water uptake. However, spatial characterization of root water uptake data is generally lacking tosupport multidimensional root water uptake parameters.
Moreover, available uptake models are largely limited to one
dimension only [Feddes et al., 1976; Molz, 1981; Jarvis, 1989]
and describe variations in water uptake with soil depth whileallowing for reduction in uptake by soil water stress. Exceptions are the two-dimensional models proposed by Neuman etal. [1975], Warrick et al. [1980], Coelho and Or [1996], and mostrecently, Vrugt et al. [2001b].
In the past few years, computing capabilities have significantly improved the effectiveness of multidimensional soil water flow models to study spatial and temporal patterns of root
water uptake. Such a multidimensional approach in root water
uptake is needed if uptake is varying in space, thereby allowing
a more accurate quantification of spatial variability of the soil
water regime, including the water and solute flux densities
below the rooting zone. The objective of this study is threefold.
First, we test the suitability of a three-dimensional model forthe simultaneous, dynamic simulation of soil water flow androot water uptake. A three-dimensional finite element grid
over the considered soil domain serves to define the spatial
distribution of soil physical properties and root characteristics
and acts as a framework for the transient water flow model.
Soil physical and root parameters are subsequently estimatedusing inverse modeling while using the measured spatial distribution of water contents around a sprinkler-irrigated almond tree during a 16 day period. A second objective was tocompare the results of the three-dimensional analysis withnumerical models describing soil water flow and root water
uptake in one and two dimensions, with uptake parametersoptimized using the same field data set. Finally, the third goalof the simulation study is to evaluate the improved predictionof the spatial variability of soil water flux taking into consideration the multidimensionality of root water uptake. Whereas
verification of the presented multidimensional uptake model islimited to a single data set, the main goal of this study is toemphasize the importance of multidimensional root water uptake modeling in root zone domains with spatially distributed
root water uptake.

2.

Materials and Methods

2.1. One-, Two-, and Three-Dimensional Root Water
Uptake Model
Recently, Vrugt et al. [2001b] proposed a one-dimensionalroot water uptake model, which was shown to be very flexible.
It was based on the model by Raats [1974]␤ 共 z兲 ⫽

冋

1⫺

册

z
e ⫺共 pz/Zm兲兩z*⫺z兩,Zm(2)

where ␤ ( z) is a shape factor describing the spatial distributionof potential root water uptake with depth, Z m (L) is the
maximum rooting depth, and p z and z* (L) are empiricalparameters. These parameters are included to provide for zeroroot water uptake at z ⫽ Z m to account for asymmetrical root
water uptake with depth and also to allow for a maximum root
water uptake rate at any depth, Z 0 (0 ⱕ Z 0 ⱕ Z m ). The
asymmetry in root water uptake with soil depth is determinedby the ratio between p z for z ⱕ z* and the p z value for z ⬎ z*.
To reduce the number of parameters, p z is set to unity forvalues of z ⬎ z*, whereas it is a fitted value for z ⱕ z*. The
value of Z 0 , the depth of maximum uptake, can simply becalculated from the first derivative, or冏d ␤ 共 z兲⫽ 0.
dz z0Assuming axial symmetry in root water uptake while using
the same root water uptake model leads to the following twodimensional root water uptake model:␤ 共r, z兲 ⫽

冉

1⫺

z
Zm

冊冉

1⫺

冊

r
e ⫺关共 pz/Zm兲兩z*⫺z兩⫹共 pr/Rm兲兩r*⫺r兩兴,Rm(3)

where R m (L) is the maximum rooting length in the radial
direction, r (L) is the radial distance from the origin of thetree, and p r (dimensionless) and r* (L) are additional empirical parameters. Here ␤ (r, z) (dimensionless) denotes thetwo-dimensional spatial distribution of potential root water
uptake. Vrugt et al. [2001b] showed that the root water uptakein (3) is extremely flexible and allows spatial variations of water

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

uptake as influenced by nonuniform water application (e.g.,
drip irrigation) and root length density patterns.
Including an additional exponential term in (3) leads to athree-dimensional root water uptake model, which is expressed
as

冉
冉

x
Xm

␤共x, y, z兲 ⫽ 1 ⫺
䡠 1⫺

冊冉
冊

1⫺

y
Ym

冊

z ⫺关共 p x /X m 兲兩x*⫺x兩⫹共 p y /Y m 兲兩y*⫺y兩⫹共 p z /Z m 兲兩z*⫺z兩兴
e
,
Zm

(4)

where X m and Y m (L) are the maximum rooting length in thex and y direction; x and y (L) are the distances from the originof the tree in the x and y direction; p x (dimensionless), p y(dimensionless), x* (L), and y* (L) are empirical parameters;and ␤ ( x, y, z) (dimensionless) denotes the three-dimensional
spatial distribution of potential root water uptake. As in (1) we
set p x , p y , and p z to unity for z ⬎ z*, x ⬎ x* and y ⬎ y*,respectively.
Denoting the normalized root water uptake S m (L 3 L ⫺3⫺1T ) as the volume of water extracted per unit volume of soil
with time, it follows that
One-dimensional description
S m共 z兲 ⫽冕␤ 共 z兲T pot
Zm

,

(5a)

␤ 共 z兲 dz

0

Two-dimensional axial symmetry

␲ R 2␤ 共r, z兲T pot

S m共r, z兲 ⫽

冕冕
Zm

2␲

0

,

Rm

(5b)

r ␤ 共r, z兲 dr dz

0Three-dimensional description
S m共 x, y, z兲 ⫽X mY m␤ 共 x, y, z兲T pot

冕冕冕
Xm

0

Ym

0

.

Zm

(5c)

␤ 共 x, y, z兲 dx d y dz

0

So integration of any of the above expressions over the spatial
domain leads to the result that cumulative potential root water
uptake is equal to the potential transpiration T pot.
To provide for root water uptake under water-stressed conditions, a soil water stress response function was included [Van
Genuchten, 1987],

␥ 共h兲 ⫽冋 冉
1⫹1
h共 x, y, z, t兲
h 50

冊册
p

,

(6)

where h is the soil water matric head at a particular spatial
location, h 50 (L) is the soil water pressure head at which root
water uptake rate is reduced by 50%, and p (dimensionless) isa fitting parameter. The parameter p is usually assumed to be
3 [Van Genuchten and Gupta, 1993].
Finally, the actual root water uptake rate at any particularspatial location can be calculated from
S共h, x, y, z兲 ⫽ ␥ 共h兲S m共 x, y, z兲,
where for an almond tree,(7)

T pot ⫽ ETalmond ⫺ E s

2459

(8)

⫺1

where S(h, x, y, z) (T ) is the actual root water uptake andE s (L T ⫺1 ) denotes soil evaporation. ETalmond defines thepotential ET by the almond crop and is computed from theproduct of K c and ET0, where K c is the crop coefficient (dimensionless), and ET0 (L T ⫺1 ) is the reference evapotranspiration. Hence the actual transpiration rate T a can be computed from
One-dimensional (1-D)
Ta ⫽

冕

Zm

S共h, z兲 dz,

(9a)

0

Axial symmetrical (2-D)
Ta ⫽

2␲
␲R2

冕冕
Zm

0

Rm

rS共h, r, z兲 dr dz,(9b)

S共h, x, y, z兲 dx d y dz.

(9c)

0

Three-dimensional (3-D)
Ta ⫽
2.2.

1
X mY m

冕冕冕
Xm

0

Ym

0

Zm

0

Field Description and Measurements

The experimental plot of our almond orchard includes about
one quarter of the wetted area of a microsprinkler irrigating asingle almond tree [Koumanov et al., 1997; Vrugt et al., 2001b].
Hence, despite the spatial variability in irrigation amounts by
the microsprinkler, we assumed that this instrumented area of
2.0 ⫻ 2.0 m (Figure 1) was representative. Twenty-five PVC
neutron probe access tubes were installed in a square grid of
0.50-m spacing to a depth of 90 cm. The neutron probe wascalibrated from gravimetric measurements using soil samples
collected at soil depths of 15, 30, 45, 60, 75, and 90 cm during
and after access tube installation. Separate calibration curves
were used for the 0 –15 cm surface soil and the 30 –90 cm soil
depth interval. Standard errors of the estimate of volumetricwater content curves were ⬃0.01 m3 m⫺3 at the 15 cm depth
interval and 0.02 m3 m⫺3 at all other measurements depths.
The field is slightly undulating, and the soil is a shallow gravelly
loam [Andreu et al., 1997], overlaying a sloping high-density
restricting clay layer at about the 90 –120 cm soil depth. The
studies by Andreu et al. [1997] and Koumanov et al. [1997]
indicated that root water uptake during the growing season was
mainly limited to the top 40 – 60 cm and that drainage at the 90
cm depth occurred primarily by lateral flow along the sloping
restricting clay layer.
The measurements were conducted during the September
13–29 period in the summer of 1995, after the microsprinkler
system was used to wet up the whole soil profile. Neutronprobe measurements were taken on September 13, immediately after the irrigation at 1300, 1500, and 1800 LT; during theperiod September 14 through September 17, every 4 hours at
600, 1000, 1400, and 1800 LT; and during the period September 18 through September 29, daily at ⬃1000 LT. This resulted
in 31 water content measurements at each spatial location. To
test the one- and two-dimensional root water uptake model,
the three-dimensional local measurements of water content
were reduced to one and two dimensions, respectively. For one
dimension, all the water content measurements at a specificdepth for all x and y locations were arithmetically averaged.

2460

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

Figure 1. Schematic view of the experimental plot. Location of the neutron probe access tubes is indicatedby solid dots. Circles approximate neutron probe measurement volumes, and curved lines indicate theaveraging volumes for the two-dimensional simulations.

Hence the one-dimensional data set was determined from 31
water content values at each depth and measurement time
during the September 13–29 monitoring period.
To simplify testing of the two-dimensional root water uptake
model, the three-dimensional grid measurements of water content needed to be reduced to two dimensions (r and z). For
this we assumed that (1) the root water uptake around the tree
was axisymmetrical and (2) the measurement volume of theneutron probe water content measurements was a sphere witha constant radius of ⬃0.25 m.
For each depth interval the rectangular measurement grid of
Figure 1 was partitioned into five concentric 0.6 m wide circular strips with their origin determined by the neutron accesspipe location closest to the tree trunk. Second, a radial averagewater content value was computed for each of the five soil
areas (0.2– 0.8, 0.8 –1.4, 1.4 –2.0, 2.0 –2.6, and 2.6 –3.2 m concentric circles) using weighting factors for each neutron probe
location with values equal to the fraction of the measurement
volume fitting within the respective concentric soil area. We
used 0.6 m wide strips for each of the five soil areas to ensurethat enough water content measurements were contained
within the respective strip. Moreover, the averaging using the
0.6 m wide strips gave the best agreement in total water depletion of the reduced two-dimensional domain as comparedto the original three-dimensional grid of water content measurements. Since the averaging procedure was applied to depth
intervals of 0 – 0.15, 0.15– 0.3, 0.3– 0.45, and 0.45– 0.6 m during
the September 13–29 period, the final two-dimensional map
included 20 average water content values at each measurement
time, namely, four depth intervals and five radial distanceincrements.
In Figure 2 we present the measured three-dimensional
volumetric water content distributions at three different times
during the monitoring period, as illustrated by depth interval–
averaged water content values. At the beginning of the period,
volumetric water content values are ⬃0.24 m3 m⫺3 at the soil
surface and increase downward to ⬃0.28 m3 m⫺3. Clearly, as
time proceeds, the water content at the 60 cm depth is still

relatively high, whereas a water uptake pattern is becoming
apparent in the surface layers, with volumetric water content
values ranging between 0.04 and 0.08 m3 m⫺3 in the top soil
layers.
2.3.

Soil Water Flow Modeling

The adapted models that solved for the one-, two-, andthree-dimensional solutions of (1) were HYDRUS-1D [Šimunek et al., 1998], HYDRUS-2D [Šimunek et al., 1999], andHYDRUS-3D (update of the SWMS_3D code of Šimunek et
al. [1995]), respectively. All three models use the Galerkinfinite element method based on the mass conservative iterative
scheme proposed by Celia et al. [1990].
The unsaturated hydraulic properties for all three models
are defined by [Van Genuchten, 1980; Mualem, 1976]Se ⫽␪ 共h兲 ⫺ ␪ r⫽ 关1 ⫹ 共 ␣ 兩h兩兲 n兴 ⫺m␪s ⫺ ␪rm⫽

n⫺1
n

1/m m 2
K共 ␪ 兲 ⫽ K sS 0.5
e 关1 ⫺ 共1 ⫺ S e 兲 兴 ,

(10)
(11)

where ␪ s (L 3 L ⫺3 ) is the saturated water content, (L 3 L ⫺3 )is the residual water content, ␣ (L ⫺1 ) and n (dimensionless)are curve shape parameters, and K s (L T ⫺1 ) denotes thesaturated hydraulic conductivity.
In the preprocessing phase the soil domain was discretized
into a rectangular grid of finite elements with similar size
elements defining the spatial resolution to avoid differences intruncation errors between the three simulation models. For all
three simulation models the vertical domain was 0.6 m deep. In
HYDRUS-1D the vertical domain was discretized into 40
equidistant finite elements. For the two-dimensional model thesimulated flow domain was 3 m long in the radial direction,using a grid spacing of 0.05 m in the radial and 0.015 m in thevertical direction. For the three-dimensional simulations thedimension of the soil domain was 0.15 by 2.65 m in the x-yplane. The flow domain was discretized into a structured meshof 13,500 blocks corresponding with 14,400 nodes, with a variable node spacing between 0.01 and 0.10 m. A structured mesh

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

2461

Figure 2. Measured three-dimensional volumetric water content distributions at three different times during the monitoring period.
was used to avoid large errors between measured and computed water content values at nodes with overlapping neutronprobe measurement volumes and was a compromise between
available computing time and nodal density.
Figure 3 presents the daily estimated boundary conditions asfunction of time during the monitoring period. Numericalmodeling requires estimates of potential transpiration T pot andsoil evaporation E s . As no direct measurement of the transpiration of the almond tree was available, we used the following
approach. Daily reference evapotranspiration ET0 were provided by a nearby weather station of the California Irrigation
Management Information System (CIMIS). Almond potential
ETalmond was calculated from ET0 and the appropriate crop
coefficients K c . Snyder and Pruitt [1988] recommended a value
for K c of 0.91, corresponding to conditions of 60% canopy soil
surface coverage for drip-irrigated trees in the Sacramento
Valley. Ritchie’s [1972] equation was used to estimate soil
evaporation. The radiation interception was calculated using
the empirical function for maize [Snyder et al., 1985], while weused an upper limit of stage 1 cumulative soil water evapora-

Figure 3. Soil surface boundary conditions during simulation
period (Time 0 corresponds with September 13).

tion of 6 mm and a partitioning factor of 0.4 between stage 1
and stage 2 evaporation [Ritchie, 1972]. The potential transpiration of almond trees T pot was obtained by subtracting the soil
evaporation from ETalmond (equation (7b)). Although we
agree that this approach for calculating T pot is quite simplistic,
it is to be expected that errors in estimated daily T pot amounts
are smaller than 10%. We should also note that the focus of
this study was to demonstrate the development and functionality of spatially distributed root water uptake models.
Owing to the lack of accurate flux information, we assumed
a unit hydraulic gradient at the lower boundary (gravity flow).
This approach seemed most appropriate since water balancecalculations using the estimated ET and measured infiltration
data indicated that a drainage term was required to matchmeasured soil water storage data changes. The water content
measurements immediately after the irrigation of September
13 were used as initial condition for all numerical simulations.
2.4.

Parameter Optimization by Inverse ModelingIn the inverse modeling stage of this study a total of 6 (e.g.,
Z m , p z , z*, n, K s , and h 50 ), 9 (e.g., Z m , p z , z*, R m , p r , r*,n, K s , and h 50 ), and 12 (e.g., X m , p x , x*, Y m , p y , y*, Z m , p z ,z*, n, K s , and h 50 ) root water uptake and soil hydraulicparameters were identified simultaneously using the one-,two-, and three-dimensional HYDRUS models, respectively.
Despite measurement of the soil hydraulic properties of anearby location in the same almond orchard [Andreu et al.,1997], the soil heterogeneity within the orchard led us to also
optimize some of the soil hydraulic parameters simultaneouslywith the respective root water uptake model parameters. While
fixing the parameters ␪ s and ␣ to reported values [Andreu et al.,1997] of 0.30 m3 m⫺3 and 9.4 m⫺1, respectively, the soil hydraulic properties were assumed to be characterized by thefitting parameters n and K s of (10) and (11). Since some of themeasured water content values were smaller than the residual
water content values reported by Andreu et al. [1997], theresidual water content ␪ r was fixed to 0.0 m3 m⫺3.
Since optimization algorithms such as Levenberg-Marquardt
or Simplex method are generally only applicable to identify alimited number of unique parameters, an alternative was

2462

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

Table 1. Range of Parameter Values Used With Genetic Algorithm (GA) With HYDRUS-1D, HYDRUS-2D, andHYDRUS-3D Simulations
HYDRUS-1D
Z m, mz*, mpz

h 50 , mn

K s , cm d⫺1

0.00
1.00

1.00
5.00

0.10
15.0

⫺0.20
⫺10.0

1.20
3.00

1 ⫻ 10⫺2
100.00

Minimum
Maximum

HYDRUS-2D

Minimum
Maximum

Z m, m

R m, m

x*, m

r*, m

pz

pr

h 50 , mn

K s , cm d⫺1

0.00
1.00

1.00
5.00

0.00
1.00

1.00
5.00

0.10
15.0

0.10
15.0

⫺0.20
⫺10.0

1.20
3.00

1 ⫻ 10⫺2
100.00HYDRUS-3D

Minimum
Maximum

X m, mY m, m

Z m, m

x*, m

y*, m

z*, m

px

py

pz

h 50 , mn

K s , cm d⫺1

1.00
5.00

1.00
5.00

0.00
1.00

1.00
5.00

1.00
5.00

0.00
1.00

0.10
15.0

0.10
15.0

0.10
15.0

⫺0.20
⫺10.0

1.20
3.00

1 ⫻ 10⫺2
100.00needed to optimize the larger set of parameters of this study.
Recently, it has been shown that genetic algorithms (GA) area powerful tool for parameter identification, if the number offitted parameters is large [Bäck, 1996; Wang, 1991; Holland,1975]. The genetic algorithm is a search procedure based on
the mechanics of natural selection and natural genetics thatcombines an artificial survival of the fittest with genetic operators [Holland, 1975]. The GA differs from other search methods as it searches among a population of parameter sets rather
than the parameter values themselves using probabilistic transition rules. We applied the GA presented by Penny and Lindfield [1995] with the small adaptation that the best performing
parameter combination is not mutated in the next generation.
We used a crossover percentage of 85 to ensure a relatively fast
convergence to the global optimum, whereas a mutation factorof 0.15 was used to avoid optimized solutions in local minima.
The population size, representing the number of first generation parameter combinations, was set to 120, whereas the final
optimized parameter combination was selected after 200 generations. Assuming that the residuals, representing errors between measured and optimized volumetric water content values, are Gaussian distributed, independent, and homocedastic(constant variance), the fitness of a chromosome was calculated by the following objective function (OF):冘N

OF(b) ⫽关 ␪ *共t i兲 ⫺ ␪ 共t i, b兲兴 2,

(12)

i⫽1

where N is the number of observations and ␪ *(t i ) and ␪ (t i , b)
denote the measured and predicted water content values, respectively, at time t i . The parameter vector b characterizes thechromosome with the genes representing the fitting parameters. The allowable ranges of the parameters included in b foreach numerical model are presented in Table 1.
Although GAs are an effective means of reaching the global
minimum region, they are not necessarily efficient in finding
the exact optimum location. Therefore the results of the genetic algorithm were used as initial values for a subsequent
Simplex algorithm (SA) to determine the local minimum of OF
within the global minimum region as determined by the GA.
Using a sensitivity analysis in which each parameter was variedwith 10% of its final optimized value, while keeping the additional parameters fixed at their value determined by the GA,
only those six parameters that were most sensitive to model
output were fine-tuned. Both the GA and Simplex optimization were carried out using MATLAB, version 5.3 [TheMathWorks, 1999].
The uncertainty of each optimized parameter bj , j ⫽ 1, 䡠 䡠 䡠 ,m, was determined from the diagonal elements of the parameter covariance matrix C [Kool and Parker, 1988; Šimunek andHopmans, 2001], representing the estimate of the standard
deviation s j ,s j ⫽ 冑C jj,(13)

whereas model performance was evaluated by the root meansquared error (RMSE), computed fromRMSE ⫽

冑冘N关 ␪ *共t i兲 ⫺ ␪ 共t i, b兲兴 2

i⫽1

N⫺m

,

(14)

where N and m denote the number of measurements and total
number of parameters, respectively.
A single forward simulation of the HYDRUS-3D model
required between 5 and 60 min on a PIII 466 MHz computer,depending on the parameter combination provided. So thecomputational time for one computer to perform all 24,000
model runs for the GA optimization iterations was extremely
long. Instead, we used 40 PIII 400 MHz slave computers connected with one master computer to perform the optimizations.

3.

Results and Discussion

3.1.

Three-Dimensional Simulations

The parameter vector of the best performing “chromosome”
of the final population after 200 generations of the threedimensional root water uptake parameters using GA in combination with the HYDRUS-3D flow model is presented inTable 2. Also included are the final results after final tuning of
the selected parameters using the SA with their confidenceVRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

2463

Table 2. Optimized Parameter Values and Their 95% Confidence Intervals After Genetic Algorithms (GA) and Simplex
Algorithm (SA) for the HYDRUS-3D Modela
HYDRUS-3D
Parameter
X m, mY m, m
Z m, m
px
py
pz
x*, m
y*, m
z*, m
n
K s , cm d⫺1
h 50 , m
RMSE, m3 m⫺3R2
X 0, m
Y 0, m
Z 0, m
a

GA
3.31
2.65
0.43
1.86
2.62
2.57
1.92
1.92
0.35
1.72
1.82
⫺0.85
0.0183
0.91
1.53
1.64
0.26

95% Confidence IntervalSA
Root Model
2.60
0.45

2.21
0.40
Soil Hydraulic Model
1.74
1.60
0.0180
0.92

Lower

Upper

CVb

2.84
2.43
0.32
1.67
2.46
2.24
1.61
1.64
⫺0.06

3.77
2.78
0.57
2.06
2.77
2.89
2.22
2.42
0.86

7.10
3.46
13.33
5.38
2.86
6.23
7.81
4.75
57.50

1.71
0.10
⫺1.05

1.78
3.10
⫺0.64

1.15
46.88
12.35

Derived Parameters
1.53
1.61
0.27

Open space means parameter held constant to value found by genetic algorithm.
CV is coefficient of variation.

bintervals, values for the derived parameters X 0 , Y 0 , and Z 0 ,and RMSE and R 2 values. The optimum maximum rooting
depth Z m (0.43 ⬍ Z m ⬍ 0.45 m) is in excellent agreementwith the results obtained by Koumanov et al. [1997] for thesame experimental plot, confirming that active root water uptake was limited to the top 40 cm only. The position of maximum root water uptake of the almond tree under nonstressed
conditions (X 0 ⫽ 1.53, Y 0 ⫽ 1.61, Z 0 ⫽ 0.27 m) agrees wellwith the surface area of maximum irrigation application by
microsprinkling [Koumanov et al., 1997]. This was so despite
the location of the microsprinkler at the far corner along thetree row (see Figure 1), as caused by nonuniform water applications during the growing season. The optimized n value ofthe soil hydraulic functions (n ⫽ 1.74) agrees well with thereported n value (n ⫽ 1.44 ⫺ 1.99) for this soil obtainedwith the instantaneous profile method [Andreu et al., 1997].
The relatively high value of the h 50 parameter (h 50 ⫽ ⫺0.85m) is an indication of the small water holding capacity of thiscoarse-textured soil. The optimized saturated conductivity(K s ⫽ 1.60 cm d⫺1) is lower than the reported range of
34.1– 62.4 cm d⫺1 measured under saturated conditions for the0 – 60 cm depth interval [Andreu et al., 1997]. The lower optimized K s value is to be expected since the measured experimental conditions were such that the soil rooting zone was lessthan saturated, thereby eliminating the influence of the macropores on the estimated K s . Using the same data set for their
two-dimensional analysis, Vrugt et al. [2001b] showed that the
HYDRUS-2D model was well able to predict water content
dynamics for spatial locations not included in the calibration
period.
The measured water content values for all depths and measurement locations are correlated with simulated water content values using the final optimized parameter values in Figure 4. Measured values match the simulated values with an R 2
value of 0.92. The overall RMSE value of 0.018 m3 m⫺3 is low,

considering that the standard errors of the neutron probe water content measurements are already 0.01 m3 m⫺3 (15 cm
depth) and 0.02 m3 m⫺3 (all other soil depths).
Figure 5 presents three-dimensional maps of simulated water content and root water uptake intensity, S m ( x, y, z),averaged for the indicated soil compartments (0 – 0.15, 0.15–
0.30, 0.30 – 0.45, and 0.45– 0.60 m) at the three different timesof September 18, 23, and 29. These maps were obtained fromarithmetic averaging over all nodal values within each depth
interval followed by interpolation using SURFER [Golden
Software, 1996]. At the beginning of the period, after irrigation
on September 13, maximum actual water uptake rates approached 8 ⫻ 10⫺4 m3 m⫺3 hr⫺1. As the soil becomes depletedFigure 4. Measured versus simulated soil water contents
around the almond tree obtained using the calibrated
HYDRUS-3D model.

2464

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

Figure 5. Simulated three-dimensional volumetric water content and potential root water uptake distributions at three times during the monitoring period.
in water, regions of maximum root water uptake shifted toother locations within the rooting zone where soil water was
most readily available [Green and Clothier, 1999]. For example,
close observation of Figure 5 shows that the general root water
uptake pattern changes with time from maximum uptake
around (X 0 , Y 0 , Z 0 ) toward the outside perimeter of therooting volume, as caused by changes in soil water stress withtime. Although there are differences in the spatial pattern ofsoil water content between simulated and measured water content values (compare Figures 2 and 5), the magnitude of simulated water contents at the different depth intervals agreed
well with the measured water contents.
In Figure 6a we present contour plots of the time-averaged
RMSE of water content at the four depth intervals. These
contour plots were computed using arithmetic averaging of all
nodal values within the respective soil volume of the corresponding neutron probe measurement and thus resulted in 25
RMSE values for each depth interval. Although, in general,
RMSE values are small, relatively large error values arepresent at the bottom corner in the spatial domain near x ⫽2.4 m and y ⫽ 0.4 m. Differences between measured andsimulated water content values are likely because of model
errors as caused by restrictive assumptions regarding the geometry of the rooting system, homogeneity of soil hydraulicproperties within the spatial domain, and the prescribed root
water uptake model. For example, the water uptake model
assumes a single region of maximum uptake, whereas in reality
more regions within the rooting zone may show local maximum
uptake as caused by water application nonuniformities and soil
environmental factors affecting root growth. Figure 6b presents a box plot of the time averaged RMSE values for all 100

spatial locations (25 tubes ⫻ 4 depth intervals). The box plotshows the single outlier with RMSE ⫽ 0.04 at 1.9 ⬍ x ⬍ 2.4
m in the bottom right-hand corner of the spatial domain as well
as a clear clustering of RMSE values between 0.01 and 0.02 m3
m⫺3, with their magnitude about equal to the standard error of
the water content measurements with the neutron probe. The
various horizontal lines represent the minimum, maximum,
25th percentile, mean, and 75th percentile of the RMSE values, respectively.
3.2. Dimensional Effects on Parameter Optimization
Results
Final optimized parameter values after using the GA andfine-tuning with the SA for the HYDRUS-1D, HYDRUS-2D,and HYDRUS-3D models and their 95% confidence intervals
are presented in Table 3. Also included are the derived parameter values of X 0 , Y 0 , and Z 0 and the fitting results as
expressed by the RMSE and R 2 values. Since many errors may
occur, including measurement, model, and numerical errors,
an uncertainty analysis of the optimized parameters makes up
an important part of parameter estimation. Therefore we included 95% confidence intervals for the optimized parameters,
calculated using the Jacobian matrices and residuals for the
final optimized solution.
As Table 3 shows, optimized parameter values for the maximum rooting depth Z m (0.41– 0.49 m), n (1.74–1.91), h 50(⫺0.53 to ⫺0.85 m), and Z 0 (0.27– 0.28 m) are close between
the three cases. It indicates that the information content of thewater content measurements, whether aggregated or not, is themost robust for these parameters [Vrugt et al., 2001a]. To
obtain convergence of the inverse solution it was essential toVRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

2465

Figure 6. (a) Three-dimensional spatial distribution of root mean squared errors of water content. (b) Box plot of RMSE for
all measured locations and times.
select initial values for the other parameters, close to their
expected values [Inoue et al., 1998]. Problems with nonuniqueness of these parameters are caused by the presence of numerous local minima, as can occur when this many parameters are

optimized simultaneously [Duan et al., 1992]. Whereas all our
optimization results for the numerical models indicate thedepth of maximum root water uptake to be at ⬃0.25 m, thestudy by Andreu et al. [1997] concluded that maximum uptake

Table 3. Optimized Parameter Values and Their 95% Confidence Regions After Simplex Algorithm (SA) for the
HYDRUS-1D and HYDRUS-2D Modelsa
HYDRUS-1D
Parameter
X m, mY m, m
Z m, m
R m, m
px
py
pz
pr
x*, m
y*, m
z*, m
r*, m
n
K s , cm d⫺1
h 50 , m
RMSE, m3 m⫺3R2
X 0, m
Y 0, m
Z 0, m
R 0, ma

Optimized Value(Confidence Interval)

HYDRUS-2D
CV

Optimized Value(Confidence Interval)

HYDRUS-3D
CV

Root Model
0.49 ⫾ 0.045

2.98 ⫾ 0.74

4.21

12.49

0.27 ⫾ 0.04

7.09

1.91 ⫾ 0.22
0.72 ⫾ 0.19
⫺0.76 ⫾ ⫺0.13
0.0068
0.98

13.56
8.76
5.89

0.41 ⫾ 0.49
3.99 ⫾ 0.49

10.76
6.22

3.32 ⫾ 4.14
3.44 ⫾ 0.55

12.64
8.18

0.31 ⫾ 0.13
1.96 ⫾ 0.18

20.79
4.53

Soil Hydraulic Model
1.75 ⫾ 0.11
0.75 ⫾ 0.16
⫺0.53 ⫾ ⫺0.16
0.01540.91
Derived Parameters

0.27

For completeness we also report the HYDRUS-3D model results.

0.28
1.96

3.26
10.93
15.10

Optimized Value(Confidence Interval)

CV3.31 ⫾ 0.46
2.60 ⫾ 0.18
0.45 ⫾ 0.13

7.10
3.46
13.33

1.86 ⫾ 0.20
2.62 ⫾ 0.15
2.57 ⫾ 0.32

5.38
2.86
6.23

1.92 ⫾ 0.30
2.21 ⫾ 0.21
0.40 ⫾ 0.46

7.81
4.75
57.50

1.74 ⫾ 0.04
1.60 ⫾ 1.50
⫺0.85 ⫾ ⫺0.21
0.0180
0.92

1.15
46.88
12.35

1.53
1.61
0.27

2466

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

Figure 7. Optimized soil water retention curves using HYDRUS-1D, HYDRUS-2D, and HYDRUS-3D
flow models with corresponding multidimensional root water uptake models. Also included are independently
measured ( ␪ , h) points [Andreu et al., 1997].

occurred at the soil surface (0 –15 cm) and decreased further
down the soil profile. However, their study did not include soil
evaporation as a possible mechanism of soil water depletion
near the soil surface. The coefficient of variation (CV) indicates that for most parameters the confidence intervals at theminimum are typically small (CV ⬍ 10%).
The optimized soil water retention and unsaturated soil hydraulic functions of the HYDRUS-1D, HYDRUS-2D, andHYDRUS-3D models are presented in Figure 7. Also included
are the measured ( ␪ , h) data using the multistep outflow
method from soil cores taken at the 30 cm soil depth for a
nearby location. Both the measured ( ␪ , h) points and the
optimized retention curves clearly show the small water holding capacity of this shallow gravely soil. Whereas the optimized
retention functions of the HYDRUS-2D and HYDRUS-3D
model match extremely well, the HYDRUS-1D optimization
shifts the water retention to an even more coarser-textured
soil, increasing drainage as compared to the multidimensional
flow models. Differences between independently measured ( ␪ ,h) points and the optimized curves mostly occur outside therange of experimental water content values (␪ ⬍ 0.08) and aremost likely caused by the sensitivity of the optimized retention
curve to the fixed residual water content value ( ␪ r ⫽ 0).
Clearly, parameters obtained with parameter estimation are to
be used only within the measurement range for which they
were determined [Inoue et al., 2000; Vrugt et al., 2001a]. Additionally, the optimized unsaturated hydraulic conductivity andmeasured (K, ␪ ) points show the rapid decrease of the hydraulic conductivity with decreasing water content. The optimized saturated conductivity of 0.46 cm d⫺1 of the SA optimization is much lower than the reported range of 34.1– 62.4 cm
d⫺1 measured under saturated conditions for the 0 – 60 cm
depth interval by Andreu et al. [1997]. However, one should
realize that the saturated hydraulic conductivity in this study ismuch more a water balance parameter, controlling the magnitude of the lower boundary flux, than it is a soil physicalparameter affecting soil water flow in the soil domain. Moreover, the saturated conductivity determined by Andreu et al.
[1997] was measured under saturated conditions, when macropores play a major role.
Although not presented, the parameter correlation matrix
for the one-, two-, and three-dimensional models using theoff-diagonal terms of the covariance matrix showed that correlations were typically low for the three-dimensional (3-D)

optimizations but that parameter correlations increased whendecreasing spatial dimensions (2-D and 1-D). For example, thecorrelation between the parameters n and h 50 was high (R ⫽0.99) for the 2-D model. Low parameter correlations areimportant as they increase the likelihood of uniqueness of the
final solution. The decrease in parameter correlation with increasing flow dimensions is likely caused by the corresponding
increase in number of observations used in the optimization,
relative to the increasing number of fitting parameters. As
shown in Figure 8, the simulated water contents with the final
parameter estimates obtained with the SA compare favorablywith the corresponding measured water contents for both the
HYDRUS-1D and the HYDRUS-2D model. The increase in
RMSE with increasing model dimension is caused by the significantly increased number of water content observations (3-Dversus 2-D and 1-D) included in the objective function (equation (12)) with increasing spatial dimension of the optimization
problem. This is especially the case as the number of fitting
parameters increased only slightly with spatial dimension.
Moreover, the averaging of the water contents used in the oneand two-dimensional simulations reduced the general watercontent variability, thereby decreasing the final RMSE values.
Figure 9 presents a comparison between the optimized spatial distributions of potential root water uptake over the spatial
domain as obtained using the one-, two-, and three-dimensional root water uptake and flow models. Starting with thethree-dimensional model results, spatial values of the optimized potential root water uptake function, ␤ ( x, y, z), werearithmetically averaged and subsequently normalized in theradial direction (equation (5b)) to obtain an average S m (r, z).
Subsequently, a similar averaging procedure was carried out inthe radial direction to show the average S m ( z). Using thetwo-dimensional flow and root uptake model, spatially averaging and normalization (equation (5a)) of S m (r, z) resulted inanother average S m ( z). The close match between the three
one-dimensional potential root water uptake distributions inFigure 10 indicates that the suggested aggregation of measured
three-dimensional water content measurements to arrive at
two- and one-dimensional root water uptake models is valid.
In Table 4 we present the different water balance components with their spatial variations (standard deviations) as simulated by the HYDRUS-1D, HYDRUS-2D, andHYDRUS-3D models. Spatial variability in drainage flux wascomputed from simulated flux density values at the 55 cm soil

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

2467

Figure 8. Measured versus simulated soil water contents around the almond tree after parameter optimization using HYDRUS-1D and HYDRUS-2D flow models.

depth. The agreement between the various listed water balancecomponents for the different numerical models is satisfactory.
As different numerical models are used with differences inaggregation of water content values between spatial dimensions, perfect agreement is unlikely. Cumulative soil evaporation and drainage components are typically small compared tototal cumulative root water uptake. Differences between actualand potential transpiration are caused by water stress (differences between T pot and T a ). As soil water storage is used in
the objective function of (12), its value must remain approximately equal between simulations. Consequently, changing
drainage amounts compensates for differences in root water
uptake between 1-, 2-, and 3-D simulations. Most importantly,
the standard deviation results in Table 4 shows that the spatial
variation in drainage rate and root water uptake decreases

when reducing multidimensional soil water flow and root water
uptake to decreasing spatial dimensions. This may have large
implications for chemical transport in root zones, as drainage
rates and corresponding chemical transport rates will vary according to root water uptake distribution.
Figure 10 presents a detailed two-dimensional contour plotof the spatial variability of cumulative flux density (mm) during
the September 18 –29 monitoring period at the 55 cm soil
depth, as computed from the HYDRUS-3D model. Cumulative net soil water flow is downward (positive values), except
for a small portion of the rooting zone domain at ( x, y) ⫽(1.4, 2.5). Although the results in Figure 10 are influenced by
the choice of the lower boundary condition, Figure 10 clearly
demonstrates that the spatial variability of the drainage rate
below the rooting zone is large, with values increasing as cor-Figure 9. Comparison of optimized spatial distributions of potential root water uptake in one, two, andthree dimensions.

2468

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING

Table 5. Components of the Soil Water Balance (andStandard Deviation) Obtained Using One-Dimensional Root
Water Uptake Model in Multidimensional Flow Model

T pot, mm
E, mm
T a , mm
Drainage, mm⌬ storage, mm
RMSE, m3 m⫺3
R2

Figure 10. Two-dimensional contour plot of spatial variability in cumulative drainage at the 0.55 m soil depth during
monitoring period.

responding root water uptake values decrease. The increasing
accurate spatial description of root water uptake and soil waterflow with increasing spatial dimension is essential to improvemodel predictions of water fluxes and contaminant transport
through the vadose zone. Moreover, total chemical load to thegroundwater will depend on local concentration and fluxes andtheir spatial variability. Specifically, although the averagechemical load can be small, as computed from average flux andconcentration values in one-dimensional simulations, the actual chemical load can be much larger. For example, this is thecase if local regions of high drainage rates, as controlled by low
root water uptake, correspond with high concentration values.
3.3. One-Dimensional Root Water Uptake Model in
Multidimensional Flow Modeling
The final analysis investigates the need for a multidimensional root water uptake approach in multidimensional flow
modeling. For this purpose, we compare simulation results
when including a one-dimensional root water uptake model inthe multidimensional water flow modeling. For this purpose
the soil hydraulic parameters (h 50 , n, and K s ) and root parameters (Z m , p z , and z*) were taken from the optimized
HYDRUS-1D model. The final fitting results are expressed by
the RMSE and R 2 values between measured and simulated
water content values. These values, combined with the waterbalance components, including their spatial variability as determined by standard deviation values, are presented in Table
5. As expected, correlations are low (R 2 values of ⬃0.6), andTable 4. Components of the Soil Water Balance Using
HYDRUS-1D, HYDRUS-2D, and HYDRUS-3Da

T pot, mm
E, mm
T a , mm
Drainage, mm⌬ storage, mm

HYDRUS-1DHYDRUS-2DHYDRUS-3D

48.70
6.99
42.04
7.78
56.81

48.70
6.99
36.53 (35.18)
8.62 (5.33)
52.18

48.70
6.49
39.45 (39.60)
13.06 (8.11)
59.00

aParentheses give standard deviation. Variables are cumulative potential almond transpiration T pot, soil evaporation E, actual root water
uptake T a , drainage, and change in storage.

HYDRUS-2DHYDRUS-3D

48.70
6.99
42.46 (37.48)
10.80 (4.58)
60.25
0.0553
0.61

48.70
6.25
41.97 (38.02)
11.34 (6.44)
59.57
0.0627
0.57

RMSE values are large (0.05– 0.06 m3 m⫺3), indicating that a
one-dimensional root water uptake model is unable to capturesoil moisture variations caused by spatially variable root water
uptake. Moreover, whereas total root water uptake and its
variation are almost identical to the comparable HYDRUS-1D
values in Table 4, total drainage was increased and spatial
variability decreased. The inability to predict within root zone
soil moisture variability using a one-dimensional root water
uptake model in multidimensional rooting systems clearly favors the need for multidimensional root water uptake and flow
models, if detail at this scale is required.

4.

Conclusions

In this paper, we have tested the suitability of a threedimensional model for the simultaneous, dynamic simulationof soil water flow and root water uptake. After parameter
optimization of the selected root water uptake model and soil
hydraulic parameters, the agreement between simulated andmeasured water contents values during the 16 day period wasgood, with an overall time-averaged root mean squared error
value of 0.018 m3 m⫺3. These results are excellent, bearing in
mind that the standard error of the water content measurements was between 0.01 and 0.02 m⫺3 m⫺3. Subsequently,
using the same field data set of multidimensional volumetricwater content values, the results of the three-dimensional root
water uptake model were compared with inverse modelingdata, describing root water uptake and soil water flow in two
and one dimensions. Independently measured soil water retention data agreed favorably with the optimized retention curves
using either one-, two-, or three-dimensional root water uptake
with corresponding multidimensional water flow models. The
high value of the optimized water stress parameter agreed withthe low water holding capacity of the sandy field soil.
Optimized root water uptake distributions between one-,two-, and three-dimensional flow models with corresponding
root water uptake models were almost identical. These results
provide evidence that the presented spatial aggregation of soil
moisture data is adequate for calibration purposes to arrive ateffective root water uptake parameters. Also, when comparing
water balance components between the three models, all models were in approximate agreement. However, major differences occurred for the spatial variation in root water uptake
and drainage rates between one-dimensional and multidimensional models. This loss of information regarding variability ofdrainage rates and root water uptake clearly justifies the needfor multidimensional root water uptake and flow models, especially when the fate and transport of chemicals below therooting zone for single trees is of concern.

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELING
Acknowledgments. We acknowledge K. Koumanov for providing
the spatially distributed water content data. We thank the Land andWater Resources Research and Development Corporation(LWRRDC) and CSIRO Land and Water, Australia, making possible
the sabbatical leave of J.W.H. Specifically, we thank K. L. Bristow of
the CSIRO Davies Laboratories in Townsville, Australia, for his efforts in this regard as well as for his valued comments and discussions.
The Earth Life Sciences and Research Council (ALW) partly supported the investigations of J.A.V. with financial aid from the Netherlands Organization for Scientific Research (NWO). The authors
acknowledge the constructive review by B. Clothier, the associate editor, and one anonymous reviewer.

References
Andreu, L., J. W. Hopmans, and L. J. Schwankl, Spatial and temporal
distribution of soil water balance for a drip-irrigated almond tree,
Agric. Water Manage., 35, 123–146, 1997.
Bäck, T., Evolutionary Algorithms in Theory and Practice: Evolution
Strategies, Evolutionary Programming, Genetic Algorithms, Oxford
Univ. Press, New York, 1996.
Canadell, J., R. B. Jackson, J. R. Ehleringer, H. A. Mooney, O. E. Sala,
and E. D. Schulze, Maximum rooting depth of vegetation types at
the global scale, Oecologia, 108, 583–595, 1996.
Celia, M. A., E. T. Bouloutas, and R. L. Zarba, A general massconservative numerical solution for the unsaturated flow equation,
Water Resour. Res., 26, 1483–1496, 1990.
Clausnitzer, V., and J. W. Hopmans, Simultaneous modeling of transient three-dimensional root growth and soil water flow, Plant Soil,
164, 299 –314, 1994.
Clothier, B. E., and S. R. Green, Rootzone processes and the efficientuse of irrigation water, Agric. Water Manage., 25, 1–12, 1994.
Coelho, F. E., and D. Or, A parametric model for two-dimensional
water uptake intensity by corn roots under drip irrigation, Soil Sci.
Soc. Am. J., 60, 1039 –1049, 1996.
Duan, Q., S. Sorooshian, and V. Gupta, Effective and efficient global
optimization for conceptual rainfall-runoff models, Water Resour.
Res., 28, 1015–1031, 1992.
Feddes, R. A., P. Kowalik, K. Kolinska-Malinka, and H. Zaradny,Simulation of field water uptake by plants using a soil water dependent root extraction function, J. Hydrol., 31, 13–26, 1976.
Gardner, W. R., Dynamic aspects of water availability to plants, Soil
Sci., 89, 63–73, 1960.
Gardner, W. R., Relation of root distribution to water uptake andavailability, Agron. J., 56, 41– 45, 1964.
Gardner, W. R., Dynamic aspects of soil-water availability to plants,
Annu. Rev. Plant Physiol., 16, 323–342, 1965.
Gardner, W. R., and C. F. Ehlig, Some observations on water movement to plant roots, Agron. J., 54, 453– 456, 1962.
Golden Software, version 6.04c, Golden, Colorado, 1996.
Green, S., and B. Clothier, The root zone dynamics of water uptake bya mature apple tree, Plant Soil, 206, 61–77, 1999.
Holland, J. H., Adaptation in Natural and Artificial Systems: An Introductory Analysis With Applications to Biology, Control, and Artificial
Intelligence, Univ. of Mich. Press, Ann Arbor, 1975.
Homaee, M., Root water uptake under non-uniform transient salinity
and water stress, Ph.D. thesis, pp. 169, Agric. Univ. Wageningen,Wageningen, Netherlands, 1999.
Inoue, M., J. Šimùnek, J. W. Hopmans, and V. Clausnitzer, In situ
estimation of soil hydraulic functions using a multistep soil-water
extraction technique, Water Resour. Res., 34, 1035–1050, 1998.
Inoue, M., J. Šimùnek, J. S. Shiozawa, and J. W. Hopmans, Simultaneous estimation of soil hydraulic and solute transport parametersfrom transient infiltration experiments, Adv. Water Resour., 23(7),677– 688, 2000.
Jarvis, N. J., A simple empirical model of root water uptake, J. Hydrol.,
107, 57–72, 1989.
Kool, J. B., and J. C. Parker, Analysis of the inverse problem fortransient unsaturated flow, Water Resour. Res., 24, 817– 830, 1988.
Koumanov, K. S., J. W. Hopmans, L. J. Schwankl, L. Andreu, and A.
Tuli, Application efficiency of micro-sprinkler irrigation of almondtrees, Agric. Water Manage., 34, 247–263, 1997.
Molz, F. J., Models of water transport in the soil-plant system: A
review, Water Resour. Res., 17, 1245–1260, 1981.
Molz, F. J., and I. Remson, Extraction term models of soil moistureuse by transpiring plants, Water Resour. Res., 6, 1346 –1356, 1970.

2469

Mualem, Y., A new model for predicting the hydraulic conductivity ofunsaturated porous media, Water Resour. Res., 12, 513–522, 1976.
Musters, P. A. D., and W. Bouten, Assessing rooting depths of anAustrian pine stand by inverse modeling soil water content maps,
Water Resour. Res., 35, 3041–3048, 1999.
Neuman, S. P., R. A. Feddes, and E. Bresler, Finite element analysis of
two-dimensional flow in soil considering water uptake by roots, I,
Theory, Soil Sci. Soc. Am. J., 35, 224 –230, 1975.
Nyer, E. K., and E. G. Gatliff, Phytoremediation, Ground Water Monit.
Rem., 16(1), 58 – 62, 1996.
Penny, J., and G. Lindfield, Numerical Methods Using MATLAB, 328
pp., Ellis Horwood, Chichester, England, 1995.
Raats, P. A. C., Steady flows of water and salt in uniform soil profiles
with plant roots, Soil Sci. Soc. Am. Proc., 38, 717–722, 1974.
Ritchie, J. T., Model for predicting evaporation from a row crop withincomplete cover, Water Resour. Res., 8, 1204 –1213, 1972.
Šimùnek, J., and J. W. Hopmans, Parameter optimization and nonlinear fitting, in Methods of Soil Analysis, Part I, 3rd ed., Monogr. 9,edited by G. C. Topp and J. H. Dane, Am. Soc. of Agron., Madison,Wis., in press, 2001.
Šimùnek, J., K. Huang, and M. T. van Genuchten, The SWMS_3D
code for simulating water flow and solute transport in threedimensional variably saturated media, Version 1.2, Res. Rep. 139,
155 pp., U.S. Salinity Lab., U.S. Dep. of Agric., Agric. Res. Serv.,
Riverside, Calif., 1995.
Šimùnek, J., M. Šejna, and M. T. van Genuchten, The HYDRUS-1D
software package for simulating the one-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media, Version 2.0, 202 pp., Rep. IGWMC-TPS-70, Int. Ground Water Model.
Cent., Colo. Sch. of Mines, Golden, 1998.
Šimùnek, J., M. Šejna, and M. T. van Genuchten, The HYDRUS-2D
software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media, Version
2.0, Rep. IGWMC-TPS-53, 251 pp., Int. Ground Water Model. Cent.,
Colo. Sch. of Mines, Golden, 1999.
Snyder, R. L., and W. O. Pruitt, Crop coefficients, in Irrigation Scheduling—A Guide for Efficient On-Farm Water Management, edited by
D. A. Goldhamer and R. L. Snyder, Pub. 21454, 67 pp., Univ. of
Calif. Div. of Agric. and Nat. Resour., Davis, 1988.
Snyder, R. L., D. W. Henderson, W. O. Pruitt, and A. Dong, California
irrigation management information system final report, vol. I, Dep.
of Land, Air, and Water Resour., Univ. of Calif., Davis, 1985.
Somma, F., J. W. Hopmans, and V. Clausnitzer, Transient threedimensional modeling of soil water and solute transport with simultaneous root growth, root water, and nutrient uptake, Plant Soil, 202,281–293, 1998.
The MathWorks, MATLAB, version 5.3, MathWorks, Natick, Mass.,
1999.
Van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
Van Genuchten, M. T., A numerical model for water and solute movement in and below the root zone, Res. Rep. 121, U.S. Salinity Lab.,
Agric. Res. Serv., U.S. Dep. of Agric., Riverside, Calif., 1987.
Van Genuchten, M. T., and S. K. Gupta, A reassessment of the crop
tolerance response function, Bull. Indian Soc. Soil Sci., 4, 730 –737,
1993.
Vrugt, J. A., W. Bouten, and A. H. Weerts, Information content ofdata for identifying soil hydraulic parameters from outflow experiments, Soil Sci. Soc. Am. J., 65, 19 –27, 2001a.
Vrugt, J. A., J. W. Hopmans, and J. Šimunek, Calibration of a twodimensional root water uptake model for a sprinkler-irrigated almond tree, Soil Sci. Soc. Am. J., in press, 2001b.
Walton, B. T., and T. A. Anderson, Microbial degradation of trichloroethylene in the rhizosphere: Potential application to biologicalremediation of waste sites, Appl. Environ. Microbiol., 56, 1012–1016,
1990.
Wang, Q. J., The genetic algorithm and its application to calibrating
conceptual rainfall-runoff models, Water Resour. Res., 27, 2467–
2471, 1991.
Warrick, A. W., D. O. Lomen, and A. A. Fard, Linearized moistureflow with root extraction for three dimensional, steady conditions,
Soil Sci. Soc. Am. J., 44, 911–914, 1980.
Whisler, F. D., A. Klute, and R. J. Millington, Analysis of steady stateevapotranspiration from a soil column, Soil Sci. Soc. Am. Proc., 32,
167–174, 1968.

2470

VRUGT ET AL.: ROOT WATER UPTAKE FUNCTIONS FOR TRANSIENT MODELINGJ. W. Hopmans, Hydrology Department, Department of Land, Air
and Water, University of California, 1 Shields Avenue, 123 Veihmeyer
Hall, Davis, CA 95616, USA. (jwhopmans@ucdavis.edu)
J. Šimunek, USDA Salinity Laboratory, University of California,
Riverside, CA 95207, USA.
M. T. van Wijk and J. A. Vrugt, Institute for Biodiversity andEcosystem Dynamics, University of Amsterdam, Nieuwe Achtergracht
166, Amsterdam 1018 WV, Netherlands.

(Received October 23, 2000; revised April 18, 2001;
accepted May 8, 2001.)