Vol. 65

JULY–AUGUST 2001

No. 4

DIVISION S-1—SOIL PHYSICS
Calibration of a Two-Dimensional Root Water Uptake Model
J. A. Vrugt, J. W. Hopmans,* and J. Šimunek
ABSTRACT

mates of the temporal and spatial root water uptake
patterns are needed. Quantification of root water extraction rates also contributes to an improved understanding of chemical fluxes in the vadose zone in both
ecological and hydrological studies (Somma et al., 1998),
as well as their control by vegetation. Water uptake byrooting systems can control the timing and the amountof chemical pollutant loadings to the groundwater
through elimination of preferential flow patterns of water and chemicals, or by regulation of absorption ofnutrients or trace elements, thereby reducing their concentration levels in the deep vadose zone or groundwater (Clothier and Green, 1994). Moreover, the rhizosphere might be responsible for accelerated breakdownof organic chemicals by biodegradation (Walton andAnderson, 1990).
Actual root water uptake not only depends on theroot distribution and its functioning, but also on soil
water availability and salinity. In addition to water stressin 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 andwater stress. In these regions, soil and water management practices are based on maintaining a favorable
soil water content and salinity status in the root zone,
thereby minimizing periods of water stress while controlling leaching to minimize salinity stress.
The influence of plant–root systems on water andchemical 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; Musters, 1998). The

Although solutions of multidimensional transient water flow canbe obtained by numerical modeling, their application may be limited
in part as root water uptake is generally considered to be one-dimensional only. The objective of this study was to develop and test a twodimensional root water uptake model, which can be incorporated
into numerical multidimensional flow models. The two-dimensional
uptake model is based on the model by Raats, but is extended witha radial component. Subsequently, the root water uptake model wasincorporated into a two-dimensional flow model, and root water uptake parameters were optimized, minimizing the residuals between
measured and simulated water content data. Water content was measured around a sprinkler-irrigated almond tree (Prunus laurocerasus
M.J.Roem) for a 16-d period at 25 locations, following irrigation. To
calibrate the flow and root water uptake model, a genetic algorithm(GA) was used to find the approximate global minimum of the optimized parameter space. The final fitting parameters were determinedusing the Simplex algorithm (SA). With the optimized root water
uptake parameters, simulated and measured water contents during
the 16-d period were in excellent agreement, with R2 values generally
ranging between 0.94 and 0.99 and a root mean squared error (RMSE)
of 0.015 m3 m⫺3. The developed root water uptake model is extremely
flexible and allows spatial variations of water uptake as influenced by
nonuniform (drip irrigation) and uniform water application patterns.

F

rom a hydrological perspective, water uptake byroot systems and their spatial distribution can largely
control water fluxes to the atmosphere and the groundwater (Canadell et al., 1996). For an improved understanding of the magnitude of these fluxes, accurate esti-

J.A. Vrugt, Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam, The Netherlands, Nieuwe Prinsengracht 130,
Amsterdam, 1018 VZ; J.W. Hopmans, Hydrology Program, Dept.
Land, Air and Water Resources (LAWR), University of California,
Davis, CA 95616, USA; J. Šimunek, USDA Salinity Laboratory, University of California, Riverside, CA 95207, USA. Received 18 Nov.
1999 *Corresponding author (jwhopmans@ucdavis.edu).

Abbreviations: CIMIS, California Irrigation Management Information System; CV, coefficient of variation; GA, genetic algorithm;
RMSE, root mean squared error; SA, simplex algorithm.

Published in Soil Sci. Soc. Am. J. 65:1027–1037 (2001).

1027

1028

SOIL SCI. SOC. AM. J., VOL. 65, JULY–AUGUST 2001

most common approach for modeling root water uptake
in unsaturated flow is based on introducing a sink term,
S, in the Richards equation (Whisler et al., 1968; Molz,
1981; Clausnitzer and Hopmans, 1994) describing transient multidimensional water flow:

water contents around a sprinkler-irrigated almond tree
during a 16-d monitoring period.

⳵␪⫽ ⵜ[Kⵜ(h ⫺ z)] ⫺ S
⳵tRoot Water Uptake Model[1]

where ␪ is the volumetric water content (L3 L⫺3 ); K is
the unsaturated hydraulic conductivity (L T⫺1 ); h (L)is the soil water pressure head; z (positive downwards)
denotes the gravitational head (L) (to be included forthe vertical flow component only); and S is the volumetric sink term (L3 L⫺3 T⫺1 ), being a function of both
space and time. However, application of multidimensional flow models requires the spatial characterizationof root water uptake as well. Available uptake models
are largely limited to one dimension only (Feddes etal., 1974; Molz, 1981; Jarvis, 1989), describing variations
in water uptake with soil depth while allowing for reduction in uptake by soil water stress. Exceptions are thetwo-dimensional models proposed by Neuman et al.
(1975) and Warrick et al. (1980). However, their application to describe different types of root distributions isfairly limited. Most recently, Coelho and Or (1996) proposed bivariate gaussian root distribution density functions (normal, semilognormal and lognormal) as parametric models for the description of root water uptake
patterns of drip-irrigated row crops.
It is the objective of this study to develop a flexible
multidimensional root water uptake model to be integrated in the two-dimensional HYDRUS-2D flow code
(Šimunek et al., 1999). Furthermore, the root watermodel was calibrated using the spatial distribution of

MATERIALS AND METHODSAs basis of the proposed root water uptake model, we usedthe exponential model by Raats (1974),1 ⫺z␤(z ) ⫽ e ␭␭[2]where ␤(z) is the spatial root water uptake distribution withdepth (L⫺1 ); ␭ (L) is selected such that at depth ␭ the cumulative root water uptake is 63% of total uptake over the wholeroot zone; and z is the depth in the soil profile (z ⱖ 0).
The proposed model excludes ␭, but includes three additionalparameters

冤

冥

pz

z ⫺z |z*⫺z|e m;
␤(z ) ⫽ 1 ⫺
zm

zⱖ0

[3]

where ␤(z) denotes the dimensionless spatial root distribution
with depth; zm is the maximum rooting depth (L); and pz (-)
and z* (L) are empirical parameters. These parameters were
included to provide for zero root water uptake at z ⱖ zm, toaccount for nonsymmetrical root water uptake with depth andto allow for maximum root water uptake at any depth, zmax
(0 ⱕ zmax ⱕ zm ). The nonsymmetry in root water uptake withsoil depth is determined by the ratio between pz for z ⱕ z*, andthe pz value for z ⬎ z*. To reduce the number of parameters, pz
is set to unity for values of z ⬎ z*, whereas it is a fitted value
for z ⱕ z*. The value of z ⫽ zmax can be calculated from thefirst derivative, i.e., d␤/dz ⫽ 0.
As the potential cumulative root water uptake must equal
the potential transpiration rate (Tpot ), the normalized root
water uptake distribution, Sm (T⫺1 ), with depth is computed

Fig. 1. Representation of different root water uptake models, Sm(z), with depth.

1029

VRUGT ET AL.: CALIBRATION OF A TWO-DIMENSIONAL ROOT WATER UPTAKE MODEL

Table 1. Parameter values for different water uptake distributions under nonstress conditions.

from

␤(z)Tpot
Sm(z) ⫽ z⫽zm冮 ␤(z)dz[4]z⫽0

Since in most field studies, zm is known a priori, the uptake
model of Eq. [3] contains only two unknown parameters (pzand z*). Figure 1 shows six different possible configurationsof normalized root water uptake distribution, using zm equals
1.0 m. The corresponding parameter values for the different
models are listed in Table 1. The first four root water uptake
models (A, B, C, and D) have maximum root water uptake at
the soil surface (zmax ⫽ 0), whereas the other two uptake
distributions (E and F) simulate maximum uptake at zmax values of 0.2 and 0.5 m, respectively, as may be caused by subsurface drip irrigation.
For the characterization of the uptake intensity along theradial direction we used a similar expression as Eq. [3];冤

␤(r ) ⫽ 1 ⫺冥

prr ⫺r |r*⫺r
e m ;
rm

rⱖ0

[5]

where ␤(r) characterizes the dimensionless spatial distributionof unstressed root water uptake in the r direction; r m is the
maximum rooting length in the radial direction (L); r is theradial distance from the origin of the plant (L) and pr and r *are empirical parameters with units (-) and (L), respectively.
As in Eq. [3] we set pr to unity for r ⬎ r *.
Combining the uptake intensity along the z-direction (Eq.
[3]) with the uptake intensity along the radial direction leads

Model
A
B
C
D
E
F

Literature (adapted from)Raats, 1974
Prasad, 1988
Hoffman and van Genuchten, 1983

z*

pz

m
0.00
1.00
0.20
0.20
1.00
1.00

-†
0.01
1.00
10.0
1.00
2.00

† Model does not include parameter similar to pz.

to a two-dimensional root water uptake model, which can be
expressed as

冤冢

冣冥冤冢

z
␤(r,z ) ⫽ 1 ⫺
zm

冣冥

pz

pr

r ⫺ z |z*⫺z|⫹r |r*⫺r|
冣
m
1⫺
e 冢m
rm

[6]

where ␤(r,z) denotes the dimensionless two-dimensional spatial distribution of root water uptake.
Both the presented water uptake model and the bivariate
gaussian density functions presented by Coelho and Or (1996)
contain six parameters. However, the proposed model includes at least two parameters (zm and r m ) with a physicalmeaning. In contrast to the model of Coelho and Or (1996),
Eq. [6] can be directly evaluated in the limit as z →0 and r →0. The single expression of Eq. [6] can simulate a wide varietyof root water uptake patterns, whereas Coelho and Or (1996)
introduce different distribution functions that need to be evaluated between uptake patterns.
Denoting the normalized root water uptake, Sm, as the vol-Fig. 2. Four different configurations of two-dimensional spatial distribution of potential root water uptake, ␤(r,z ).

1030

SOIL SCI. SOC. AM. J., VOL. 65, JULY–AUGUST 2001

Table 2. Parameter values for the two-dimensional root water
uptake configurations in Fig. 2.
fitting parameters
Figure

zm

rm

z*

r*

pz

1.00
1.00
1.00
1.00

2.00
2.00
2.00
2.00

pr

zmax

r max

(-)
0.00
0.20
0.00
0.20

0.00
0.00
1.00
1.00

1.00
1.00
1.00
5.00

1.00
1.00
4.00
2.00

0.00
0.20
0.00
0.20

0.00
0.00
1.00
1.00

␲R 2␤(r,z) Tpotzm rm

冤 冢 冣冥
h p
h50

[8]

where h is the soil water pressure head at spatial location (r,z ); h50 is the soil water pressure head (L) at which root water
uptake rate is reduced by 50%; and p is a fitting parameter(-). The parameter p is usually assumed to be 3 (van Genuchten
and Gupta, 1993).
Finally, the actual root water uptake rate at (r,z ) can becalculated fromume of water extracted per unit time and volume of soil whileassuming axial symmetry, it follows that (Šimunek et al., 1999)

Sm(r,z) ⫽1
1⫹

derivedmA
BCD

␥(r,z,h) ⫽S(h,r,z) ⫽ ␥(r,z,h)Sm(r,z)[9a]Tpot ⫽ KcET0 ⫺ Es

[9b]

and

[7]

2␲ 冮 冮 r␤(r,z)drdz
0 0

where Sm(r,z ) denotes the normalized root water uptake rate(T⫺1) and R is the size of the flow domain in the r-direction.
If r m ⬍ R then the value of R in Eq. [7] automatically equals r m.
Four different types of two-dimensional root water uptake
patterns for unstressed conditions are presented in Fig. 2,
which were simulated with the proposed model in Eq. [6].
The corresponding parameters of the different water uptake
distributions are listed in Table 2. A root water uptake intensity pattern as presented in Fig. 2a,b is expected for uniform
water application. Much different uptake patterns are shownin Fig. 2c (surface irrigation) and Fig. 2d (subsurface irrigation), where the location of maximum uptake intensity shifts
to locations with maximum irrigation application rate for nonuniform water applications. The proposed root water uptake
model can be easily extended to three dimensions by including
an additional exponential term in Eq. [6]. Moreover, the model
can be adapted to account for root growth by allowing timedependent zm and r m values during a growing season.
To provide for root water uptake under water-stressed conditions, a soil water stress response function was included (vanGenuchten, 1987).

where S(h,r,z ) is the actual water uptake (T⫺1); Kc is the crop
coefficient (-); ET0 is the reference evapotranspiration (LT⫺1);
and Es denotes soil evaporation (LT⫺1). Hence, the actual
transpiration rate (Ta ) can be computed from Eq. [10].
zr2␲ m mTa ⫽冮 冮 rS(h,r,z)drdz
␲R 2 0 0[10]

Field Description and Measurements
The experimental plot is located in an almond orchard, andcovers about one quarter of the wetted area of a microsprinkler irrigating a single almond tree (Koumanov et al., 1997).
In the 2.0 m by 2.0 m instrumented area, 25 polyvinyl chloride
neutron probe access tubes were installed in a square grid of
0.50 m spacing to a depth of 90 cm (Fig. 3). The neutron probe
was calibrated 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- to 15-cm surface soil and the 30to 90-cm soil depth interval. Standard errors of estimate ofvolumetric water content curves were approximately 0.01 (15cm depth interval) and 0.02 m3 m⫺3 (all other measurements

Fig. 3. Schematic view of the experimental plot (after Koumanov et al., 1997).

1031

VRUGT ET AL.: CALIBRATION OF A TWO-DIMENSIONAL ROOT WATER UPTAKE MODEL

depths). The field is slightly undulating, and the soil is a shallow gravely loam (Andreu et al., 1997), overlaying a sloping
high density restricting clay layer at about the 90- to 120-cm
soil depth. The study by Andreu et al. (1997) indicated thatroot water uptake during the growing season was mainly limited to the top 60 cm, and that drainage of the soil primarily
occurred by lateral flow along the sloping restricting clay layer.
The measurements were carried out 13 September through29 September in the summer of 1995. First, the sprinkler system was used to moisten the whole soil profile. Neutron probe
measurements were taken on 13 September, immediately after
the irrigation at 1300, 1500, and 1800 h, during the period of
14 Septenber through 17 September, every 4 h at 600, 1000,
1400, and 1800 h, and during the period of 18 September
through 29 September, daily at 苲1000 h.
To simplify testing of the root water uptake model, thethree dimensional grid measurements of water content needed
to be reduced to two dimensions (r and z ). For this we assumed
that (i) the root water uptake around the tree was axisymmetrical and (ii) that the measurement volume of the neutron probe
water content was a sphere with a constant radius of 苲0.25 m.
Although the first assumption was rather arbitrary, there wasno reason to expect the contrary.
First, for each depth interval the rectangular measurement
grid of Fig. 3 was partitioned into five concentric adjacent
0.6-m wide circular strips with the origin of the circles defining
these soil strips. The soil strips were determined by the neutron
access pipe location closest to the tree trunk (see Fig. 3).
Second, a radial-average water content value was computedfor 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) 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 compared tothe original three-dimensional grid of water content measurements. Since the averaging procedure was applied to depth
intervals of 0 to 0.15, 0.15 to 0.3, 0.3 to 0.45, and 0.45 to
0.6, the final two-dimensional map included 20 average watercontent values at each measurement time (four depth intervals
and five radial distance increments) during the 13 September
to 29 September calibration period.

Water Flow Simulation
While assuming axial symmetry for an isotropic soil, theRichards equation for a rigid porous media can be written as(Inoue et al., 1998; Šimunek et al., 1999)

C

冢

冣

and

冦 冤

冧

m

K(Se) ⫽ KsSle 1 ⫺ 1 ⫺ Se(1⁄m)冥

2

[13]

where ␪s is the saturated water content (L3 L⫺3); ␪r is theresidual water content (L3 L⫺3); ␣(L⫺1), n, and l are curve
shape parameters (-); and Ks is the saturated hydraulic conductivity (L T⫺1). The simulated flow domain was 3.0 m long(radial direction) by 0.6 m deep, using a grid spacing of 0.05 m
in the radial and 0.05 m in the vertical direction.

Boundary and Initial Conditions
For the soil surface boundary conditions, HYDRUS–2D
requires estimates of the potential transpiration (Tpot ) and soil
evaporation (Es ). Daily reference evapotranspiration (ET0 )data was provided by a nearby weather station of the California Irrigation Management Information System (CIMIS).
Almond potential ETalm was calculated from ET0 and the appropriate crop coefficients (Kc ). Snyder et al. (1988) recommended a value for Kc of 0.91, corresponding to conditionsof 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 (Zea mays
L.) (Snyder et al., 1985), while we used an upper limit of Stage
1 cumulative evaporation of 6 mm and a partitioning factorof 0.4 between Stage 1 and Stage 2 evaporation (Ritchie, 1972).
The potential transpiration of almond trees (Talm) was obtainedby subtracting Es from ET (Eq. [9b]). Figure 4 presents thedaily boundary conditions as functions of time during themonitoring period that were used for the HYDRUS–2D
model simulations.
Due to the lack of flux information, we assumed a unit
hydraulic gradient as the lower boundary condition (gravityflow). The calculated two-dimensional soil moisture profile
immediately after irrigation, 13 September, was used as initialcondition for the numerical simulations.

Parameter Optimization
In addition to the root water uptake parameters zm, r m, z*,r *, pz and pr (Eq. [6]), and the stress response parameter h50(Eq. [8]), the soil hydraulic parameters n and Ks were also
optimized. Although the soil hydraulic properties of a nearby
location in the same almond orchard were reported by Andreu
et al. (1997), the large soil heterogeneity within the orchard
made it necessary to also optimize the soil hydraulic parameters simultaneously with the root water uptake model parame-冢 冣⳵h 1 ⳵
⳵h
⳵
⳵h
⳵K
rK
⫹
K
⫺
⫽
⫺ S(h,r,z)
⳵t
r ⳵r
⳵r
⳵z ⳵z
⳵z
[11]

where C is the water capacity (L⫺1); h is the soil water pressure
head (L); r is the radial coordinate; z is the vertical coordinate(positive downwards); t is time (T); and S(h,r,z) denotes root
water uptake (T⫺1). Equation [11] was solved with the HYDRUS–2D model (Šimunek et al., 1999) using the Galerkinfinite element method based on the mass conservative iterative
scheme proposed by Celia et al. (1990).
The unsaturated hydraulic properties in HYDRUS–2D are
defined by (van Genuchten, 1980; Mualem, 1976)␪ ⫺ ␪rn⫺1⫽ Se ⫽ [1 ⫹ (⫺␣h)n]⫺m; m ⫽␪s ⫺ ␪rn[12]Fig. 4. Soil surface boundary conditions during simulation period(Time 0 corresponds with 13 September).

1032

SOIL SCI. SOC. AM. J., VOL. 65, JULY–AUGUST 2001

ters. While fixing the parameters ␪s, ␪r, ␣, and l to reported
values of 0.28, 0.0, 9.4 m⫺1, and ⫺0.850, respectively (Andreu
et al., 1997), spatial variability in soil hydraulic properties was
assumed to be characterized by the fitting parameters n andKs (Eq. [12] and [13]). We fixed the parameters ␪s, ␪r, ␣,and l to their reported values to partly avoid problems withnonuniqueness of the parameter estimates, especially because
a relatively large number of root parameters are already involved in the optimization. Hence, in the calibration stage ofthis study, a combined total of nine root water uptake andsoil hydraulic parameters were optimized simultaneously.
Since optimization algorithms such as Levenberg–Marquardt or Simplex method are only generally applicable toidentify a limited number of unique parameters, an alternative
was needed to optimize this many parameters. Recently it
has been shown that GAs are a powerful tool for parameter
identification, when the number of fitted parameters is large(Bäck, 1996; Holland, 1975). Genetic algorithms were developed in evolution theory, based on the concepts of natural
selection and genetics. In this approach, variables are represented as genes on a chromosome. Genetic Algorithms featurea group of candidate solutions (population) on a response
surface. Through natural selection and using genetic operators, such as mutation and crossover, the objective of a GA
is to search for chromosomes with improved fitness, that is,a parameter set, which is closer to the global optimum in theobjective function. Natural selection guarantees that chromosomes with the best fitness will propagate in future populations. Using the crossover operator, the GA combines genes
from two parent chromosomes to form new chromosomes(children) that have a high probability of having better fitness
than their parents. Mutation allows near areas of the response
surface to be explored. Hence, GAs offer an improvement in
the fitness of the chromosomes through application of the GA
in reproduction, so that many generations will create chromosomes containing the optimized variable settings (Wang,
1991).
We applied the GA method presented in Penny and Lindfield (1995), with the small adaptation that the best performing
parameter combination is not mutated in the next generation.
We used a relatively high crossover value of 0.85 to ensure a
relatively fast convergence to the global optimum, where asa mutation factor of 0.15 was used to avoid optimized solutions
in local minima. The population size or number of possible
first generation parameter combinations was set to 150,
whereas the final selected optimized parameter combination
was determined after 100 generations. The fitness of a chromosome was calculated by:

Although GAs are an effective means of determining the
global minimum region, their use is not necessarily the mostefficient way of finding the exact optimum location. Therefore,
the results of the GA were used as initial values for the Simplex
optimization to determine the local minimum of Eq. [15]within the global minimum region. Using a sensitivity analysis
in which each parameter was varied by 10% around its true
value, while keeping the additional parameters fixed at theirvalue found by the GA, we selected those six parameters in
the parameter vector that were most sensitive to the model
output. Both the GAs and the simplex optimization were carried out using MATLAB, version 5.3 (The Mathworks, 1999).
The estimated standard deviation of each parameter bj of b
was determined from the diagonal elements of the parameter
covariance matrix C (Kool and Parker, 1988),

sj ⫽ √Cjj[15]

whereas final fitting results were expressed as RMSE values,
computed fromRMSE ⫽冪n兺 [␪*(ti) ⫺ ␪(ti,b)]2

i⫽1

[16]

n⫺m

where m is the total number of parameters.

RESULTS AND DISCUSSION
Final parameter values after performing the GA andthe Simplex optimization are presented in Table 4. This
table also includes final derived parameter values, suchas zmax, r max, and the 95% confidence intervals using theresiduals and Jacobian matrix at the final solution. The
fitting results are expressed by RMSE values. As theparameters zm, r m, z*, r*, n, and Ks showed the highest
sensitivity to the model output in the parameter solution
obtained with the GA, these six parameters were further
optimized using the SA. Because of the simultaneous
fitting of this many parameters, problems may occur
with the nonuniqueness of the parameter estimates
Table 4. Optimized parameter values after genetic algorithm(GA) and Simplex algorithm (SA). Also included are the 95%
parameter confidence intervals of the final solution and thecoefficient of variation of the final parameter estimates.

nOF(b) ⫽ 兺 [␪*(ti) ⫺ ␪(ti,b)]2

95% Confidence
interval

[14]

i⫽1

where n is the number of measurements; and ␪*(ti ) and ␪(ti,b )denote the measured and predicted water content values, respectively, at time ti. The parameter vector, b, characterizes
the chromosome with the genes representing the fitting parameters. The lower the value of objective function, OF(b ), themore fit is the chromosome. The allowable ranges of the parameters included in b are presented in Table 3 and are determined by physical constraints of possible parameter values.
Table 3. Range parameters values used in global optimization
with genetic algorithms.
Parameter
zm

rm

z*

r*

pz

m
Min
Max

0.00
0.60

0.00
5.00

pr

h50

n

Ks

0.10
15.0

m
⫺0.30
⫺50.0

1.2
4.0

m d ⫺1
1e-4
0.64

(-)
0.00
0.60

0.00
5.00

0.10
15.0

Parameter
zzm
(m)
rm
(m)
(-)
pz
pr
(-)
z*
(m)
r*
(m)
n
(-)
(cm d⫺1 )
Ks
h50
(m)
RMSE
(m3 m⫺3 )R2Derived Parameters
zmax
rmax

(m)
(m)

GA

SA

Lower

Upper

CV

0.426
4.174
3.214
2.918
0.300
2.052
1.673
0.408
⫺0.5330.0157
0.91

0.403
4.144
†
†
0.330
2.075
1.674
0.460
†
0.0154
0.92

0.326
3.463
2.319
2.412
0.189
1.862
1.545
0.350
⫺0.733

0.481
4.825
4.109
3.424
0.469
2.288
1.803
0.569
⫺0.333

%
9.6
8.3
13.9
8.7
21.2
5.1
3.9
11.9
18.7

0.294
2.053

0.278
2.075

† Open space indicates that parameter was held constant at value estimatedby GA.

VRUGT ET AL.: CALIBRATION OF A TWO-DIMENSIONAL ROOT WATER UPTAKE MODEL

caused by the presence of multiple local minima in theobjective function.
Since parameter estimation involves a variety of possible errors, including measurement, model and numerical errors, an uncertainty analysis of the optimized parameters constitutes an important part of parameter
estimation. Table 4 shows that the standard deviationsof the various parameters are typically small, revealingthat the information of the water content measurements
for most of the parameters is robust (e.g., Vrugt et al.,2001). Although a variety of errors determine the final
parameter standard deviation, the root parameters pz
and z*, have a relatively higher uncertainty as comparedto the other root parameters. Their larger values of
the coefficient of variation (CV), may be caused by thesmaller number of nodal points defining a strip in thedepth direction as compared to the radial direction. Furthermore, the 95% confidence interval of the n-parameter in the soil hydraulic functions lies within the rangeof n-values between 1.44 and 1.99 as reported by Andreu
et al. (1997) at a nearby location between the 0- and60-cm depth. The relatively high uncertainty in h50 asindicated by its CV value, is partly because of its extreme
large value, corresponding with a low water holding
capacity of this soil. Although not presented, inspectionof the parameter correlation matrix for the final solution
showed that correlations between the parameter estimates are typically low (R ⬍ 0.50), except between
parameter r m and z* (R ⫽ ⫺0.88), r m and r* (R ⫽ 0.92),
z* and r* (R ⫽ 0.88), and n with Ks (0.88).
The optimized soil water retention and unsaturated
soil hydraulic conductivity curves from the final SA
solution are presented in Figure 5. Also included are
the measured (␪,h) data using the multistep outflow
method from soil cores taken at a 30-cm depth at a
nearby location, and (K,␪) points as obtained using the
instantaneous profile method at a nearby location at
the 30-cm soil depth (Andreu et al., 1997). Both themeasured (␪,h) points and the optimized retention curve
clearly show the small water holding capacity of thisshallow gravely soil. Additionally, the optimized unsaturated hydraulic conductivity and measured (K,␪) points
show the rapid decrease of the hydraulic conductivity
with decreasing water content. Moreover, the large soil
water stress value (h50 ⫽ ⫺0.53 m) is in agreement withthe low water holding capacity of this soil. The optimized saturated conductivity of 0.46 cm d⫺1 of the SA
optimization is much lower than the reported range of
34.1 to 62.4 cm d⫺1 between 0- and 60-cm depth byAndreu et al., 1997. However, one should realize thatthe saturated hydraulic conductivity in this study, ismuch more a water balance parameter controlling themagnitude of the lower boundary flux than it is a soil
physical parameter affecting soil water flow in the soil
domain.
The range of the maximum rooting depth (zm ) as
found by the SA (0.33 m ⬍ zm ⬍ 0.48 m) is in excellent
agreement with the results obtained by Koumanov etal. (1997) for the same experimental plot, suggestingthat active root water uptake was limited to the top
40 cm only using seasonal root and soil water content

1033

Fig. 5. Optimized soil water retention and unsaturated hydraulic conductivity curves as obtained by the Simplex Algorithm optimizationand measured (␪,h ) using the multistep outflow method from soilcores taken at 30-cm depth at a nearby location and (K,␪) points,
as obtained using the instantaneous profile method at a nearby
location for the 30-cm soil depth.

observations.Whereas our optimization results indicate
the depth of maximum root water uptake is about
0.28 m, the study by Andreu et al. (1997) concluded
that maximum uptake occurred at the soil surface (0–15cm) and decreased further down the soil profile. However, their study did not include soil evaporation as apossible mechanism of soil water depletion near the soil
surface. The optimized radial position of maximum root
water uptake (r max ) of the almond tree (2.07 m) agreed
well with the region of highest irrigation application
amounts of the microsprinkler (Koumanov et al., 1997).
This was so, despite the location of the microsprinkler
at the far corner along the tree row (see Fig. 3), andwas caused by systematic nonuniform water applications
during the irrigation period. This finding is consistent
with the experimental data of Coelho and Or (1996),
who concluded that the applied irrigation strategy candetermine root development in both space and time.
Figure 6 presents the optimized spatially distributed
root water uptake model, ␤(r,z), as determined overthe 16-d monitoring period using the final optimized
root water uptake parameters. Clearly, the zone of maximum root water uptake is concentrated in a thin soil
layer between 0.1 and 0.35 m. In this area, the roots act
like a sink and provide for soil water potential differ-

1034

SOIL SCI. SOC. AM. J., VOL. 65, JULY–AUGUST 2001

Fig. 6. Optimized spatial distribution of potential transpiration ␤(r,z ).

ences in the soil that regulate soil water flow towards
this active soil region from above and below.
For calibration purposes, the two-dimensional flow
region of the rooting zone of the single tree was divided
into 20 blocks of 0.6 by 0.15 m. The residuals in volumetric water content of only 12 of these blocks (C1–12)
were included in the OF (Eq. 14), and served as calibration values. The eight remaining blocks are used forvalidation (V1–8). Measured and simulated average water content values for all 20 blocks as a function oftime are presented in Fig. 7. The position of the graphs
indicates the position of each block in the two-dimensional plane. Although simulated water contents using
the optimized soil hydraulic and root water uptake parameters agree remarkably well with measured watercontent values, there are some discrepancies. For instance, simulated and measured water content values
deviate considerably for blocks C7 and V6. Clearly, thepresented procedure includes restrictive assumptions,making a perfect fit unlikely. For example, we haveassumed that root water uptake is axisymmetrical. The
original three-dimensional measured water content values, however, showed that this can only be approximately true. Moreover, the root water uptake model

assumes that only a single region of maximum water
uptake exists, whereas multiple regions within a rooting
system may show maximum uptake as caused by nonuniform water application and soil environmental factors
affecting root growth.
The explained variances and time-averaged RMSE
between measured and optimized water content values
for all calibration and verification blocks of the monitoring period are presented in Table 5. The overall fit, as
computed from the RMSE values is excellent, considering that the standard error of the water content measurements with the neutron probe is between 0.01 m3 m⫺3(15-cm depth) and 0.02 m3 m⫺3 (all larger depths).
The final two-dimensional maps of simulated watercontent (m3 m⫺3 ) with corresponding root water uptake
intensity (m3 m⫺3 h⫺1) at three times during the calibration period are shown in Fig. 8, and were obtained frominterpolation of nodal values using SURFER (Golden
Software, 1996). At the beginning of this period, maximum water uptake rates approached 7.0 ⫻ 10⫺4 m3 m⫺3h⫺1. As these soil volumes become depleted in water,
regions of maximum root water uptake shift to otherlocations within the rooting zone where soil water ismore readily available. For example, close observation

Table 5. Final optimization results for both calibration (C) and validation (V) blocks.
Subregion

RMSE
R2

RMSE
R2

C1

C2

C3

C4

C5

C6

C7

C8

C9

C10

C11

C12

0.012
0.98

0.015
0.95

0.006
0.99

0.011
0.99

0.015
0.99
Subregion

0.008
0.99

0.036
0.54

0.014
0.98

0.008
0.99

0.011
0.99

0.017
0.90

0.011
1.00

V1

V2

V3

V4

V5

V6

V7

V8

Overall

0.016
0.95

0.010
0.94

0.004
0.99

0.007
0.99

0.008
0.99

0.034
0.05

0.004
0.97

0.011
1.00

0.01540.92

Fig. 7. Measured (open circles) and simulated (solid lines) water content values as a function of time across the measured spatial domain of the almond tree (C denotes calibration and V is validation).

1036

SOIL SCI. SOC. AM. J., VOL. 65, JULY–AUGUST 2001

Fig. 8. Two-dimensional maps of simulated water content values (m3 m⫺3) and root water uptake (d⫺1) at three different times during theexperimental period.

of Fig. 8 shows that root water uptake patterns changewith time towards regions further away from the centerof the initial highest water uptake
Finally, in Fig. 9a, we present the different components of the soil water balance as calculated by the
HYDRUS–2D model. The total simulated change insoil water storage of the axial symmetric two-dimensional grid during the period 13 September through29 September was 1692 L, which corresponds to a soil
water depletion of 52.8 mm for the 0.6 by 3.0 m soil
domain. Soil evaporation (13.2%) and root water uptake (69.7%) account for 82.9% of the total soil water
storage change, whereas the remaining depletion (17.1%)
is caused by drainage. The difference between cumulative root water uptake and cumulative potential transpiration is caused by soil water stress, which mainly occurs

at the end of the study period, as the rooting zone
becomes further depleted of water. Over the whole domain during the experimental period, the reduction intranspiration for Tpot is 23.5%, with the difference between potential and actual evapotranspiration being
about 11 mm. Although soil water storage is not explicitly included in the objective function, there was excellent agreement between measured and simulated soil
water depletion, as shown in Fig. 9b.

CONCLUSIONSIn this paper we have presented a two-dimensionalroot water uptake model which was based on the modelby Raats (1974). The developed root water uptake
model is extremely flexible and allows spatial and tem-

Fig. 9. (a) Components of the water balance, (b) Measured versus simulated soil water storage during the monitoring period.

VRUGT ET AL.: CALIBRATION OF A TWO-DIMENSIONAL ROOT WATER UPTAKE MODEL

poral variations of water uptake as influenced by uniform and nonuniform water application patterns. Using
the optimized flow and uptake parameters, the simulated and measured water contents values during the
16-d period were in excellent agreement, with R2 generally ranging between 0.94 and 0.99, and an average deviation of 0.015 m3 m⫺3. These results are excellent bearing
in mind that the standard error of the water content
measurements was between 0.01 m3 m⫺3 and 0.02 m⫺3m⫺3. Some of the optimized root water uptake parameters such as maximum rooting depth agreed well withexperimental observations in the same experimental
plot. The relatively high value of the optimized soil
water stress parameter agrees with the low water holding capacity of the gravely soil. Maximum root water
uptake rates approached 7.0 ⫻ 10⫺4 m3 m⫺3 h⫺1, withsoil regions of maximum root water uptake shifting toother locations within the rooting zone where soil water
is more readily available.
ACKNOWLEDGMENTS
The authors acknowledge A. Tuli and G. Schoups for stimulating discussions. The work period of the first author in the
Department of Land, Air, and Water Resources of the University of California in Davis was made possible through financial
support by the Dr Hendrik Muller’s Vaderlandsch Fonds andthe Schuurman Schimmel-van Outeren Stichting. We acknowledge Dr. K. Koumanov for providing his water content data.

REFERENCES
Andreu, L., J.W. Hopmans, and L.J. Schwankl. 1997. Spatial andtemporal distribution of soil water balance for a drip-irrigatedalmond tree. Agric. Water Manage. 35:123–146.
Bäck, T. 1996. Evolutionary algorithms in theory and practice: Evolution strategies, evolutionary programming, genetic algorithms. Oxford University Press, New York, NY.
Canadell, J., R.B. Jackson, J.R. Ehleringer, H.A. Mooney, O.E. Sala,
and E.D. Schulze. 1996. Maximum rooting depth of vegetation
types at the global scale. Oecologia 108:583–595.
Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general massconservative numerical solution for the unsaturated flow equation.
Water Resour. Res. 26:1483–1496.
Clothier, B.E., and S.R. Green. 1994. Rootzone processes and theefficient use of irrigation water. Agric. Water Manage. 25:1–12.
Coelho, F.E., and D. Or. 1996. A parametric model for two-dimensional water uptake intensity by corn roots under drip irrigation.
Soil Sci. Soc. Am. J. 60:1039–1049.
Clausnitzer, V., and J.W. Hopmans. 1994. Simultaneous modeling oftransient three-dimensional root growth and soil water flow. Plant
Soil 164:299–314.
Feddes, R.A., E. Bresler, and S.P. Neuman. 1974. Field test of amodified numerical model for water uptake by root systems. Water
Resour. Res. 10:1199–1205.
Golden Software. 1996. Surfer, version 6.04c. Golden Software,
Golden, CO.
Holland, J.H. 1975. Adaptation in natural and artificial systems. The
University of Michigan Press, Ann Arbor, MI.
Hoffman, G.J., and M.Th. van Genuchten. 1983. Soil properties andefficient water use: Water management for salinity control. p. 73–
85. In H.M. Taylor et al. (eds) Limitation to efficient water use incrop production. Am. Soc. Agron., Madison, WI.
Homaee, M. 1999. Root water uptake under non-uniform transient
salinity and water stress. Ph.D. thesis. Agricultural University Wageningen, the Netherlands.
Inoue, M., J. Šimunek, J.W. Hopmans, and V. Clausnitzer. 1998. In

1037

situ estimation of soil hydraulic functions using a multistep soilwater extraction technique. Water Resourc. Res. 34:1035–1050.
Jarvis, N.J. 1989. A simple empirical model of root water uptake.
J. Hydrol. 107:57–72.
Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem fortransient unsaturated flow. Water Resour. Res. 24:817–830.
Koumanov, K.S., J.W. Hopmans, L.J. Schwankl, L. Andreu, and A.
Tuli. 1997. Application efficiency of micro-sprinkler irrigation ofalmond trees. Agric. Water Manage. 34:247–263.
The Mathworks. 1994. MATLAB, version 4.2.c. The Mathworks, Natick, MA.
Molz, F. 1981. Models of water transport in the soil-plant system: A
review. Water Resour. Res. 17:1245–1260.
Mualem, Y.A. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
Musters, P.A.D. 1998. Temporal and spatial patterns of root water
uptake in an Austrian pine stand on sandy soil. Ph.D. thesis. University of Amsterdam, the Netherlands.
Musters, P.A.D., and W. Bouten. 1999. Assessing rooting depths ofan Austrian Pine stand by inverse modeling soil water content
maps. Water Resour. Res. 35:3041–3048.
Neuman, S.P., R.A. Feddes, and E. Bresler. 1975. Finite element
analysis of two-dimensional flow in soil considering water uptake
by roots: I. Theory. Soil Sci. Soc. Am. J. 35:224–230.
Penny, J., and G. Lindfield. 1995. Numerical methods using MATLAB. Ellsi Horwood, New York.
Prasad, R. 1988. A linear root water uptake model. J. Hydrol. 99:297–
306.
Raats, P.A.C. 1974. Steady flows of water and salt in uniform soil
profiles with plant roots. Soil Sci. Soc. Am. Proc. 38:717–722.
Ritchie, J.T. 1972. Model for predicting evaporation from a row cropwith incomplete cover. Water Resour. Res. 8:1204–1212.
Šimunek, J., M. Šejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D software package for simulating two-dimensonal movement of water, heat, and multiple solutes in variable saturated
media. Version 2.0, IGWMC-TPS-53, International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO.
Snyder, R.L., and W.O. Pruitt. 1989. Crop coefficients. p. 67. In Goldhamer and Snyder (ed.) Irrigation scheduling—Aguide for efficient
on-farm water management. Univer. of Calf. Div. Agric. Nat. Res.
Publ. No. 21454.
Snyder, R.L., W.O. Pruitt, D.W. Henderson, and A. Dong. 1985.
California irrigation management information system final report.
Vol. I. Dept. Land, Air and Water Resourc., Univer. of Calif.,
Davis, CA.
Somma, F., J.W. Hopmans, and V. Clausnitzer. 1998. Transient threedimensional modeling of soil water and solute transport with simultaneous root growth, root water and nutrient uptake. Plant Soil
202:281–293.
van Genuchten, M.Th. 1987. A numerical model for water and solute
movement in and below the root zone. Res. Rep. 121. U.S. Salinity
Lab, ARS USDA, Riverside, CA.
van Genuchten, M.Th., and S.K. Gupta. 1993. A reassessment of thecrop tolerance reponse function. Indian Soc. Soil Sci. 4:730–737.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
Vrugt, J.A., W. Bouten, and A.H. Weerts. 2001. Information content
of data for identifying soil hydraulic parameters from Outflow
experiments. Soil Sci. Soc. Am. J. 65:19–27.
Walton, B.T., and T.A. Anderson. 1990. Microbial degradation oftrichloroethylene in the rhizosphere: Potential application to biological remediation of waste sites. Appl. and Environ. Biol. 56:1012–
1016.
Wang, Q.J. 1991. The genetic Algorithm and its application to calibrating conceptual rainfall runoff models. Water Resour. Res.
27:2467–2471.
Warrick, A.W., D.O. Lomen, and A.A. Fard. 1980. Linearized moisture flow with root extraction for three dimensional, steady conditions. Soil Sci. Soc. Am. J. 44:911–914.
Whisler, F.D., A. Klute, and R.J. Millington. 1968. Analysis of steadystate evapotranspiration from a soil column. Soil Sci. Soc. Am.
Proc. 32:167–174.