Parameter Estimation Analysis of the Evaporation Method
for Determining Soil Hydraulic Properties
Jifi Simunek,* Ole Wendroth, and Martinus Th. van Genuchten
ABSTRACT

oped (Wind, 1968; Becher, 1971; Boels et al., 1978;
Schindler, 1980; Tamari et al., 1993; Wendroth et al.,
1993; Halbertsma and Veerman, 1994).
An important modification was suggested by Wind
(1968), who introduced an iterative graphical procedure.
He first estimated the water retention characteristicfrom average water content and pressure head readings
at several locations in a homogeneous soil sample, andsubsequently determined hydraulic conductivities fromthe measured pressure head profile and changes in thewater content distribution. The water content profile
was obtained from results of the first step. The Wind
iterative method was later automated by several researchers (e.g., Boels et al., 1978; Halbertsma and Veerman, 1994). Wendroth et al. (1993) developed a methodthat required measurement of the pressure head at only
two locations. In the traditional setup, only unsaturated
hydraulic conductivities in the pressure head range fromapproximately -50 to -700 cm could be determined.
Good estimates of the conductivity close to saturation
could not be obtained because the hydraulic gradientswere too small. To overcome these problems near saturation, Wendroth et al. (1993) imposed two different
evaporation rates: initially a relatively high evaporation
rate to obtain large pressure head gradients near saturation and, after reestablishing hydraulic equilibrium,
much lower evaporation rates as controlled by the prevailing laboratory conditions. Using numerically generated data, several researchers (Tamari et al., 1993; Wendroth et al., 1993; Mohrath et al., 1997) evaluated theassumptions of the evaporation method or the effects
of experimental errors.
An alternative method of analyzing transient flow
during an evaporation experiment is to use parameter
estimation techniques (Kool et al., 1987). Methods ofthis type typically involve the coupling of a numericalmodel for variably saturated water flow with a parameter optimization algorithm. The Levenberg-Marquardt
method (Marquardt, 1963) has been especially popular
for this purpose. Starting with the studies of Zachmann
et al. (1981) and Dane and Hruska (1983), the parameter
estimation method is increasingly being used for estimating unsaturated soil hydraulic functions. Computermodels applicable to laboratory column outflow measurements have been given by Kool et al. (1985a,b) and
Parker et al. (1985) for one-step outflow procedures,and by van Dam et al. (1992, 1994) and Eching andHopmans (1993) for multistep approaches. A more general parameter estimation model applicable to transientwater flow subject to less restrictive initial and boundary
conditions was developed by Kool and Parker (1987).
Applications to evaporation experiments have been presented by Feddes et al. (1988), Ciollaro and Romano(1995), and Santini et al. (1995). In a review, Feddes etSoil hydraulic properties are important parameters affecting waterflow in variably saturated soils. We estimated the hydraulic properties
from a laboratory evaporation experiment using both a parameter
estimation method and the modified Wind method. The parameter
estimation method combined a one-dimensional numerical solutionof the Richards equation with the Marquardt-Levenberg optimization
scheme. In our study we used both numerically generated data anddata measured in the laboratory. Two experiments were carried out
on 10-cm-high soil cores containing two different soils. Pressure heads
inside the cores were measured with five tensiometers, while evaporative water loss from the top was determined by weighing the soil
samples. The objective function for the parameter estimation analysis
was defined in terms of the final total water volume in the core andpressure head readings by one or several tensiometers. An analysis
of numerically generated data showed that the optimization method
was most sensitive to the shape factor (n) and the saturated watercontent (6S) and least to the residual water content (6,). Pressure
heads measured close to the soil surface were found to be morevaluable for the parameter estimation technique than those measured
at lower locations. The optimized hydraulic parameters corresponded
closely with those obtained using Wind's analysis. All optimizations
gave similar results for the soil hydraulic properties within the rangeof measured pressure heads (0 to -700 cm). Extrapolation beyond thisrange involved a high level of uncertainty because of high correlation
between parameters 0, and n.

M

ANY LABORATORY AND FIELD METHODS exist to determine soil hydraulic properties, especially forthe unsaturated hydraulic conductivity (Klute and Dirksen, 1986; Green et al., 1986). Most methods remain
relatively time consuming and costly, and are often limited to relatively narrow ranges of water content. One
fairly simple laboratory method for simultaneous estimation of both retention and unsaturated hydraulic conductivity data for the past 30 yr has been the evaporation
method. This method was first introduced by Gardnerand Miklich (1962), who imposed a series of constantfluxes on one side of an initially equilibrated sample,
and measured the pressure head response of two tensiometers. The flux needed to be sufficiently small toassume a constant hydraulic conductivity and diffusivity
in the sample (Halbertsma and Veerman, 1994). Becher(1971) simplified the evaporation method by using acontinuous evaporation rate. Several other modifications of the evaporation method with simultaneous measurements of evaporation rate and pressure heads at
different heights in the sample have since been develJ. Simunek and M.Th. van Genuchten, U.S. Salinity Lab., USDAARS, 450 West Big Springs Road, Riverside, CA 92507-4617; O.
Wendroth, Inst. for Soil Landscape Research, ZALF, EberswalderStr. 84, D-15374 Muncheberg, Germany. Received 15 Aug. 1997.
*Corresponding author (JSimunek@ussl.ars.usda.gov).
Published in Soil Sci. Soc. Am. J. 62:894-905 (1998).

894

SIMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION

al. (1988) compared hydraulic conductivities determined from an evaporation experiment using both parameter estimation and Wind's method. While quite
successful, they applied their method to only one data
set and did not discuss problems of identifiability, stability, and uniqueness. As part of a study of the spatial
variability of soil hydraulic properties, Ciollaro and Romano (1995) used parameter estimation to determine
the soil hydraulic parameters of a large number of samples. Santini et al. (1995) also used parameter estimationin connection with evaporation experiments and compared the results with independently measured retention data and saturated hydraulic conductivities, Ks. No
attempts were made in these two studies to also compare
estimated unsaturated hydraulic conductivities with independent data. Both Ciollaro and Romano (1995) andSantini et al. (1995) assumed that 0r was equal to zero,and independently measured 0S. They defined the objective function only in terms of pressure head measurements, which prevented them from estimating 6r and 6S
directly from the evaporation experiment. In spite ofthese relatively promising results, Halbertsma (1996)
stated that evaporation experiments are not suitable forinverse modeling because of uniqueness problems.
The objective of this study was to evaluate in moredetail the potential of parameter optimization techniques for simultaneously estimating the water retention
and hydraulic conductivity relationships from an evaporation experiment. We addressed the questions of identifiability, uniqueness, and stability. The sensitivity ofan evaporation experiment to particular soil hydraulicparameters was investigated so as to evaluate the minimum amount of information needed to guarantee aunique solution. We also addressed the question whether
or not the two-rate evaporation method of Wendrothet al. (1993) provides better information for parameter
estimation than the traditional one-rate method. The
first part of our study was carried out using numericallygenerated error-free data. Numerically generated data
sets are preferred for this analysis since the true valuesof soil hydraulic parameters are known and since nomeasurement errors are present, thus causing less uncertainty in the analysis (Toorman et al., 1992). Different
sources of errors could bias some of the conclusions
that can be drawn from numerically generated data. We
next applied the parameter estimation technique to twolaboratory data sets. These data sets contained sourcesof errors typically present in experimental data sets,such as measurement errors, small inhomogeneities inthe samples, and other sources of noise, which would
make their analysis more complicated. The soil hydraulic characteristics obtained were compared with thoseusing the simplified Wind method (Wendroth et al.,1993).
METHODS
Governing Flow Equations
The governing flow equation for one-dimensional isothermal Darcian flow in a variably saturated rigid porous medium
is given by the following modified form of the Richards equa-

895

tion:
dt

[1]

dz \ dz

where 6 is the volumetric water content (L3 L~3), h is thesoil-water pressure head (L), K is the hydraulic conductivity(L T"1), z is a vertical coordinate (L) positive upward, and t
is time (T). Initial and boundary conditions applicable to an
evaporation experiment are as follows:

[2]
+ l

= <7evap(L, 0

[3]

= 0

[4]

where h, is the initial soil-water pressure head (L), qemf(t) is
the time-variable evaporation rate imposed at the soil surface(L T"1), and L is a coordinate of the soil surface. Equation[1], subject to the above initial and boundary conditions, wassolved numerically using the finite-element code HYDRUS1D as documented by Simunek et al. (1997).
The unsaturated soil hydraulic properties in this study wereassumed to be described by the following expressions (vanGenuchten, 1980):
1-e r- (1 - Sl"")m]2

[5]

[6]

where 5e is the effective water content, Ks is the saturated
hydraulic conductivity (L T"1), 6r and 9S denote the residual
and saturated water contents (L3 L~3), respectively, / is a poreconnectivity parameter, and a (L"1), n, and m (=1 — l/«) areempirical parameters. The predictive K(Q) model is based on
the capillary model of Mualem (1976) in conjunction with
Eq. [5]. The pore-connectivity parameter / in the hydraulicconductivity function was estimated by Mualem (1976) to be
0.5 as an average for many soils. The hydraulic characteristics
defined by Eq. [5] and [6] contain five unknown parameters:
6r, 6S, a, n, and Ks. The evaporation experiment in general isa drying process, which means that the hydraulic parametersin Eq. [5] and [6] represent drying branches of the unsaturated
hydraulic properties.

Wind Method
During the evaporation experiments we monitored thepressure head at selected times tt at five different depths z,(1, 3, 5, 7, and 9 cm) within the soil column (10-cm height),
and also repeatedly measured the total weight of the soil
sample. After conclusion of the experiment, the final watercontent was measured, from which the total amount of water,
Vmj (L), at each time, th could be calculated. A fourth-orderpolynomial function for log(/z)-0 was fitted to the pairs of
average soil water content and soil water pressure head dataobtained at different times according to Halbertsma and Veerman (1994). The polynomials were calculated by multiple regression. A polynomial function of fourth order was chosen
since a third-order polynomial function did not always provide
enough flexibility to follow the particular shape of the log(/z )-6
curve. The water content value 9y was subsequently calculated
for each pressure head reading for iteration number k = 1.
For a given time t/, this water content value was considered
representative for the depth compartment associated with a

896

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

particular tensiometric reading. Multiplying the thicknessesof the five compartments with the respective volumetric watercontents yielded water storage values, which, when summed,
gave the calculated total water storage, V^, at a given time.
A correction factor cj1 was determined from the ratio between
the calculated amount of water and the measured total amountof water, Vmj. This approach led to a set of; correction factors,
one for each time when the weight and tensiometer readingsof the sample were obtained. Assuming that the error invokedin the initial water retention curve was distributed equally tothe five depth compartments, the calculated water content
values 6^ were updated using[7]

to lead to new estimates of the water content, /z-6?,. The
resulting data pairs of h-Q^ were subsequently used to obtainan improved polynomial function for log(&)-6. Using thisfunction and the measured pressure heads, new water content
estimates Of/1 for iteration k + 1 were obtained, leading tocalculated water storage values that were closer to the measured values (i.e., the correction factors c/k became closer toone). With each iteration, the difference between the measured and calculated total water storage values decreased. The
solution was considered to converge when the maximum watercontent change between two iterations became insignificant(less than some small preset value). Once convergence in thewater retention curve was reached, temporal changes in thewater content of each of the five compartments were used to
calculated water flux profiles as follows
7

t,]

[8]

in which g,v (L T"1) denotes the flux from compartment i
upward into compartment / - 1 (the compartment number
increases with depth) during time interval Ar, = ti+1 - tt, where
tj and ty-H are the beginning and the end of the time interval.
The hydraulic conductivity for a compartment was subsequently derived according to Darcy:[9]+
where A/z,v is the hydraulic gradient between two consecutive
compartments (z>-\ and zi) as estimated from tensiometerreadings, and Az, = zi-\ — z,-. In Eq. [9] we used the average
gradient during a given time interval, Afy._The K value was
next plotted against the geometric mean, htj, of the pressure
head values within the two compartments at the beginningand the end of the time interval involved.

Inverse Solution
The objective function $, which is minimized during theparameter estimation process, is defined asp) = 2 v, 2 ivy \pf (t,) - Pi(th b)]2[10]

where m represents the different sets of measurements (observed pressure heads at different locations and the total water
volumes in the sample), n, is the number of measurements ina particular measurement set,p?(ti) are specific measurements
at time t, for the y'th measurement set, pfab ) are the corresponding model predictions for the vector of optimized parameters b (e.g., 6r , 6S , a, n, and Ks), and Vy and w-,j are weights
associated with a particular measurement set or point, respectively. For an evaporation experiment, the vector p* consistsof pressure heads measured at either one or several locations,and one value of the average water content in the soil sample,
usually measured at the end of the experiment. We assumed
in this study that the weighting coefficients vf,v in Eq. [10] areequal to one; that is, the variances of the errors inside ofa particular measurement set are the same. The weighting
coefficients vy are given by (Clausnitzer and Hopmans, 1995)
1
n,<rj

[11]

which shows that the objective function is defined as the average weighted squared deviation normalized by measurement
variances cr2. Since the final water volume is only one number
and the variance cannot be defined, the weight for this data
point is assumed to be one.
Minimization of the objective function <J> is accomplishedby using the effective Levenberg-Marquardt nonlinear minimization method (Marquardt, 1963), which has become a standard in nonlinear least-square fitting among soil scientists andhydrologists (van Genuchten, 1981; Kool et al., 1985a,b, 1987).
The method combines the Newton and steepest descent methods, and provides confidence intervals for the optimized parameters.

NUMERICAL EXPERIMENT
An evaporation experiment was first simulated numericallyusing average parameter values of the silty soil textural groupas estimated by Carsel and Parrish (1988) from an analysis ofa large number of soils. The soil hydraulic parameters for thishypothetical experiment are given in Table 1. The soil sample
was assumed to have a height of 10 cm, with five tensiometers
located in the soil core 1, 3, 5, 7, and 9 cm below the sample
surface. The soil profile was divided into 100 elements for our
numerical simulations using HYDRUS-1D. The thickness ofthe first element at the soil surface was 0.0205 cm. Subsequentelements increased linearly up to 0.1111 at the depth of the firsttensiometer, below which the grid spacing remained constantthroughout the sample.
As an initial condition, we assumed hydraulic equilibriumwith a zero pressure head at the bottom of the soil sample.
A zero-flux condition was imposed at the bottom boundary.
Two scenarios were used for the upper boundary condition.
For the first scenario we assumed a constant evaporation rateof 0.15 cm d~' for 14.4 d. Pressure head readings at intervalsof 0.1 d were used to define the objective function. In thesecond scenario a higher evaporation rate of 1.5 cm d"1 wasused during the first 0.5 d, followed by a zero surface flux
until 1 d, and a constant rate of 0.15 cm d"1 afterward until
the end of the simulation at 10.4 d. Shorter measurement
intervals of 0.05 d were employed for the higher evaporation
rate. The second scenario was used to evaluate whether or
not a two-rate evaporation experiment (Wendroth et al., 1993)
could benefit the parameter optimization procedure. Bothsimulations were terminated when the pressure head of theuppermost tensiometer dropped below —700 cm. Pressure
head readings for both scenarios are shown in Fig. 1. Notice
the much faster decrease in pressure head for the second (twostep evaporation) scenario at the beginning of the experiment.
Table 1. Parameter spacings used for the parameter planes of the
hypothetical evaporation experiment.
Parametere,e,a, cm *nK,, cm d-'True value

Lower
parameter
value

Parameter
step value

Upperparameter
value

0.034
0.46
0.016
1.37
6.0

0.00
0.307
0.002
1.02
0.00001

0.0035
0.007
0.002
0.02
1.0

0.1015
0.51
0.06
1.6
30.0

SIMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION

Evaporation losses during the first 5 d for the first scenario
were about equal to the evaporation loss during the first day
for the second scenario; pressure head readings afterward
were almost identical except for the time difference of 4 d.
Two-rate evaporation experiments have two important advantages compared with more traditional one-rate experiments. First, by initially increasing the evaporation rate, thetotal amount of time required for the experiment can be decreased (e.g., by 4 d in our example). Second, higher evaporation rates require higher fluxes inside the soil sample, andcorrespondingly higher pressure head gradients. By increasing
the pressure head gradients, the effects of experimental errors
associated with the pressure head readings can be reduced
substantially (Wendroth et al, 1993). This is especially important when hydraulic conductivities are calculated directly frompressure head gradients, as is the case with the Wind method
or its modifications. The question remains whether or not
adapting two-rate evaporation experiments will also benefit
numerical inversions, which use only pressure head readings
and do not require the estimation of gradients from measured
data. If such a benefit exists, then the sensitivity of pressure
head readings to the unknown soil hydraulic parameters would
have to increase significantly during the high-evaporation
stage of a two-rate evaporation experiment. We first studied
this problem by means of a sensitivity analysis.

Sensitivity Analysis
In general, an experiment should be designed such thatmeasurements are made that yield the most information about
the unknown parameters to be optimized, i.e., measurementsthat are most sensitive to changes in the unknown parameters.
One could expect that different locations of the tensiometersand measurements taken at different times will provide different degrees of information for the parameter estimation
method. The optimal positions of the tensiometers are best
evaluated by means of a sensitivity analysis.
Sensitivity coefficients, s(z,t,bj), for the hypothetical experiment were calculated from (Simunek and van Genuchten,
1996)
o
-100
-200
~ -300

J. -400^ -500
-600
-700

-800
2.5

5

7.5

10

12.5

15

Time [d]
0

b)

-100
-200
~ -300A) = 0.016,

897- 0.016= h(b + Abe,) - A(b)[12]where s(z,t,bj) is the change in the variable h (the pressurehead) corresponding to a 1% change in parameter b/, e, isthe y'th unit vector, and Ab = O.Olb. Equation [12] allows acomparison of the sensitivities to different parameters, independent of the invoked units or their absolute values. We
stress here that Eq. [12] calculates sensitivity coefficients,
which characterize the behavior of the objective function ata particular location in parameter space, presumably in thevicinity of the true parameter values. A high sensitivity in thisrespect means that the minimum is well defined, and that one
should be able to estimate the parameters with relatively highprecision once the global minimum is identified. The approach,however, does not give any information about possible local
minima in the objective function elsewhere in parameter
space.
Figure 2 shows sensitivity coefficients for all five hydraulicparameters as calculated with Eq. [12] for five different tensiometer locations in the soil core. Results are for the constantevaporation case; similar results were obtained also for thetwo-rate evaporation experiment (data not shown here). Results for both scenarios, as in Fig. 1, have similar shapes interms of the time translation of 4 d. The sensitivity coefficientsfor both scenarios increased substantially as the experiments
progressed in time, suggesting that long-duration experiments
will provide much more information for the parameter optimization process. The second scenario, however, shows only arelatively insignificant bump during the first day when theevaporation rate is high. These slightly higher sensitivity coefficients are not expected to have an important effect on theoverall parameter estimation procedure, except for short-duration experiments.
The absolute sensitivity of the pressure head readings to
all soil hydraulic parameters increased with time as the soil
profile became drier. With the exception of a, Tensiometer1 yielded much higher sensitivity coefficients to changes in
the hydraulic parameters than the other tensiometers. Figure
2 indicates that the pressure head is most sensitive to theparameter n, followed by 9S. Sensitivities to 6r, a, and Ks are,by comparison, much smaller. The results presented in Fig. 2
also suggest that the longer the measurements are taken, themore valuable the information becomes.
One can expect the sensitivity of tensiometer readings toincrease with time. As the soil becomes drier, the retention
curve becomes steeper, resulting in larger changes in pressure
head with small changes in parameter values and water contents. Figure 3 shows sensitivity coefficients expressed in termsof the water content for all five hydraulic parameters forfive different locations in the soil core. Notice that sensitivity
coefficients, again, increased with time for all parameters except a for the location closest to the soil surface. The increase
is, however, less steep due to the nonlinear nature of theretention curve. Sensitivity coefficients at other locations also
increased initially, but reached a maximum at about 10 to
12 d, after which they remained more or less constant.

J, -400
£

-500
-600
-700

-Tensiometer 1- — — Tensiometer 2------Tensiometer3- - - - Tensiometer 4- - - - Tensiometer 5

Response Surfaces

-800

2.5

7.5

12.5

Time [d]Fig. 1. Simulated pressure heads h for a silty soil for (a) one- and (b)two-step evaporation experiments.

We tested whether data measured conventionally during
an evaporation experiment (i.e., pressure head readings atseveral locations and average water contents) provide enoughinformation to enable the identification of a unique set of soil
hydraulic parameters from the inverse problem. We approached this question in a similar way as was done previously
by Toorman et al. (1992) and Simunek and van Genuchten

898

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

— - - Tensiometer 2
Tensiometer 3— - — - Tensiometer 4— • • — Tensiometer 5———Tensiometer 1— — — Tensiometer 2—-----Tensiometer 3— - — • Tensiometer 4— - - — Tensiometer 5

15

Tensiometer 1
Tensiometer 2
Tensiometer 3
Tensiometer 4
Tensiometer 5

c)
10

15

Ks

,0,

2

— — — Tensiometer 2

1— - — - Tensiometer 4— - - — Tensiometer 5

^^^

e)

""^'••^

»

/
/

0

{

5

10

1

Time [d]Fig. 2. Sensitivity of the pressure head h to parameters (a) residual soil water content, 6t, (b) saturated water content, 6(, (c) a, (d) n, and (e)saturated hydraulic conductivity, A,, as a function of time for a silty soil assuming a one-step evaporation experiment.

0.000008
0.000006"0.000004
ai

0.00000———Tens. 1— — — Tens. 2•Tens. 3— - — - Tens. 4— - - - Tens. 5b)-0.00002
r-0.00004• -0.00006

0.000002-0.00008a)0.000000-0.00010— — — Tens. 2-Tens. 3— - — - Tens. 4— • • — Tens. 5

100.00000

10

15

10

15

0.00000
Tens. 1-0.00005-0.00001,.,-0.00010* -0.00015-0.00002-0.00020-0.00003-0.00025

10— - - Tens. 2-Tens. 3— • — • Tens. 4— - - — Tens. 5

15
Ks0.000006
0.000004
0.000002; 0.000000
i -0.000002-0.000004-0.000006-0.000008--- Tens. 2•Tens. 3- - - - Tens. 4- • • — Tens. 5

510Time [d]15

Fig. 3. Sensitivity of the water content B to parameters (a) residual soil water content, 6t, (b) saturated water content, 08, (c) a, (d) n, and (e)hydraulic conductivity, K,, as a function of time for a silty soil assuming a one-step evaporation experiment.

SIMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION(1997) in studies of the one-dimensional one-step outflow
method and the axisymmetrical three-dimensional tension infiltrometer experiment, respectively.
The uniqueness of the inverse problem was evaluated interms of two-dimensional response surfaces of the objectivefunction as a function of pairs of soil hydraulic parameters.
Response surfaces were obtained by changing two selected soil
hydraulic parameters around their true values, while keeping
other parameters constant at their true values. The objectivefunction used for this purpose is given by Eq. [10]. We calculated the objective functions for 10 parameter planes (i.e.,a-n, a-Ks, n~Ks, a-6s, «-9s, A:s-9s, a-6r, «-9r, K-Qt, and 9r-9s)
for all data combined (i.e., five tensiometer readings and thefinal average water content) and for each tensiometer separately combined with the final average water content. The
response surfaces were calculated on a rectangular grid withparameter values given in Table 1. Each parameter domain
was discretized into 30 discrete points, resulting in 900 grid
points for each response surface.
We emphasize that the different parameter planes given
below represent only cross sections of the full five-dimensionalparameter space. The behavior of the objective function inthese parameter planes can only suggest how the objectivefunction might behave in the five-dimensional continuum. For
example, local minima of the objective function <S could exist
and not show up in the cross-sectional planes (Simunek andvan Genuchten, 1996). Nevertheless, the response surfaces
provide a useful approximate view of the behavior of theobjective function in the entire parameter space. The inverse
parameter estimation technique is expected to be unsuccessful
if response surfaces do not display a clearly defined global
minimum in the two-dimensional parameter planes.
Contours of the objective function <t> in 10 parameter planes
are presented in Fig. 4. The objective function <t> was definedin terms of all five tensiometer measurements and the final
total water volume in the sample. All response surfaces, withthe exception of the «-9r parameter plane (Fig. 4h), show
relatively well-defined global minima, and no additional local
minima. The structure of the response surfaces is very consistent with plots of the sensitivity coefficients presented above.
Because of the high sensitivity of the pressure head to thesaturated water content (Fig. 2), 9S yields well-defined minima
in all parameter planes (Fig. 4d, 4e, 4f, and 4j). By comparison,
the sensitivity to Ks was relatively low (Fig. 2), with the response surfaces showing also less certainty toward this parameter (Fig. 4b, 4c, and 4f). While the different minima for Ks arestill relatively well defined, many of the contours away fromthe minima are nearly parallel with the Ks axis, especially at
the higher K, values. This suggests that initial estimates of Ks
for parameter optimization should be relatively small. Pressure head measurements were also quite sensitive to n, withthis parameter again yielding relatively well-defined minima(Fig. 4a, 4c, and 4e). One exception is the n-6r parameter
plane (Fig. 4h), mainly because of uncertainty in 9r, toward
which the measurements showed the lowest sensitivity. The
parameters n and 0r were found to be highly correlated, as
will be discussed below.
Inverse Solutions
We next examined whether the global minimum of the
response surfaces could be identified numerically via inverse
simulation using the Levenberg-Marquardt optimization
method. We defined the objective function of Eq. [10] in termsof the final total water volume in the soil sample and either

899(i) all available pressure head data (at five locations and foreach time), (ii) only one tensiometer reading at a time, or (iii)
data from two tensiometers (1 and 5). For each case we chose
three different sets of initial estimates. For the one-rate evaporation experiment, all three optimizations converged when all
five or two tensiometers were used; for optimizations using
only one particular tensiometer data set, at least two runs
were always successful. The results were slightly worse for thetwo-rate evaporation experiment in that only two of the three
optimizations converged toward the true parameters for all
simulations with different definitions of the objective function.
The value of the objective function for unsuccessful optimizations was about four orders of magnitude higher than forthe successful runs, while the parameter values also yieldedphysically unrealistic values. Failure in the optimization wascaused primarily by uncertainty in the parameters 9, and n,
and partially in Ks. When 9r was fixed at its true value, or
when one (Q,h) pair for h = —150 m was included in theobjective function (Kool et al, 1985b), all optimizations forboth scenarios converged rapidly to the true parameter values.
Table 2 presents a correlation matrix for the optimized
parameters when all tensiometer readings were used for the
one-rate evaporation experiment. Note the high correlation
between n and 9r, which corresponds well with our findings
for the response surfaces (Fig. 4h). Because of this high correlation, these two parameters show a high level of uncertainty.
The uncertainty in 9r is due to the fact that within the rangeof measurement (0 to —700 cm) for the silty soil used here,
the effective water content (Se) dropped only to about 0.40,
still far from the residual water content. As will be shown
below for our laboratory data, several combinations of n and9r can yield similar retention curves within the range of measurements. Although the optimized parameters describe theretention curve well within the measurement range, extrapolation beyond this range will introduce a high level of uncertainty. A similar problem was also encountered by Parker etal. (1985) when analyzing one-step outflow data. They foundthat using only cumulative outflow data resulted in satisfactory
results within the range of water contents observed in theexperiment, but not for much lower 9. Including the watercontent at h = -150 m into the objective function extended
the range of validity of the predicted properties to the dry
end with only minor effects on the predictions at high 9 values(Parker et al., 1985). One alternative solution for this problem
is to estimate the parameter 9r independently, such as byindependent measurement or even using pedotransfer functions (Schaap et al., 1998), and then to treat 9r as a known
parameter. If estimated in this way, 9r can then also be included
into the objective function and the solution can be penalized
for any deviations from this value (Russo et al., 1991). This
technique yields Bayesian estimates of 9r (e.g., Bard, 1974;
Yeh, 1986; Simunek and van Genuchten, 1996). Inclusion of
the true value of 9r in the objective function always resulted
in response surfaces with clearly defined minima in terms of
9r (results not further shown), while also all optimizations
converged accurately and rapidly toward the true parameter values.
Not addressed in this study is the effect on the optimization
results of using initially different estimates of the parameters.
For our optimizations we selected three sets that were known
to be physically realistic for the soils used in our numerical experiments.
LABORATORY EXPERIMENT
Two undisturbed soil core samples with a height of 10 cm
and inside diameter of 10 cm were placed on a ceramic plateFig. 4. Contours of the objective function <I> in the (a) a-n, (b) a-Ks, (c) n-Kn (d) a-9s, (e) n-6s, (f) XS-6S, (g) a-6r, (h) n-0r, (i) *Ts-er, and (j) 8r-6s parameter planes. A dot represents
the true value of parameters.

r

o

i

o

o
r

w

Q

V3

SIMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION

et al. (1993). We obtained 179 and 303 data pairs
and 94 and 119 data pairs K(h) for Experiments I andII, respectively. The lowest tension for which we could
calculate the hydraulic conductivity was 26 cm. Hence,
we also included in the data sets for each experiment
three K(h) data points measured with a tension infiltrometer at pressure heads of —1, -5, and -10 cm. The
hydraulic parameters for both soils were obtained bysimultaneous fitting Eq. [5] and [6] to the resulting 6(/z)
and K(h) data using the RETC code (van Genuchten
et al., 1991). We obtained an excellent fit with R2 =0.992 and 0.963 for Experiments I and II, respectively.
There was some uncertainty in the parameters Ks andn for Experiment I because of high correlation (-0.992)
between these two parameters. The final parameters forboth experiments are given in Tables 3 and 4. Dataobtained with Wind's method and the fitted hydraulicfunctions are shown in Fig. 5 and 7 for Experiments I
and II, respectively.
The soil hydraulic parameters were also estimatedfrom the evaporation experiment using parameter inversion. We used, for this purpose, the tensiometer readings as a function of time and the total water volume
at the end of the experiment. All other water volume
measurements were used to calculate evaporation rates
needed for the upper boundary condition in the simulations. Inclusion of these measurements in the objectivefunction would have resulted in duplicate information,
without any benefit for the numerical inversion process.
Only one value of the water volume (irrespective of

Table 2. Correlation matrix for the optimized parameters for the
hypothetical evaporation experiment.
Parameter

e.
e,
a
n

K,

0,

e,

a

n

K,

1.0
-0.7772-0.54250.9850-0.23181.0
0.4903-0.79640.25491.0-0.67450.88931.0-0.3871

1.0

901

and saturated with deionized water. Soil used in the first coresample (Experiment I) had a bulk density of 1.59 g cm~3, andsand, silt, and clay fractions of 7.4, 79.3, and 13.3%, respectively. The second soil core sample (Experiment II) had abulk density of 1.58 g cm"3 and sand, silt, and clay fractionsof 45.2, 43.6, and 11.2%, respectively. Five tensiometers withcups of 6-cm length and 0.6-cm outside diameter were horizontally inserted into drill holes in the soil cores 1,3,5,7, and 9 cm
from the sample surface. The ceramic cups were connected topressure transducers. Following previous experiments on thesame soil columns with infiltrometer disks, the samples weremoved onto impermeable plates for the evaporation experiment. Initial pressure heads of -15.4 cm for Experiment I
and -19.3 cm for Experiment II were measured in the middleof the soil sample.
Evaporation was subsequently allowed to start. After eachpressure reading, transducer wires were disconnected and thesoil samples with the tensiometers weighted to determine theevaporative water loss as a function of time. At the beginningof the experiments the evaporation rate was increased to approximately 1.2 cm d"1 by using of a fan to blow air awayfrom the soil surface at room temperature (Wendroth et al.,1993). Once the gradient between the tensiometers reacheda value between 1.5 and 2.5 m m"1, the top of the soil sample
was covered to prevent further evaporation. After reestablishing hydraulic equilibrium in the samples, evaporation wasallowed to continue, without the fan, at a rate of approximately
0.2 cm d"1. Measurements were taken every 30 min during
higher evaporation rate period, and every 4 h during the lower
rates. The evaporation experiment was terminated after theupper tensiometer recorded a pressure head value below —650
cm. This limit was originally set because of the close proximityof the tensiometer measurement range, and because Wind's
method requires simultaneous readings of all installed tensiometers. For the purpose of parameter optimization study,
the experiment could have continued as long as at least one
tensiometer was functioning. The water loss between particular measurements was used to calculate the average evaporation rate for a given time interval; this information was
subsequently used as the upper boundary condition in the
numerical simulations.

when the measurement was made) is needed in theobjective function in order to position the retention
curve along the 6 axis. Without this information, neither
6S nor 0r can be estimated because of their mutual correlation. The objective function therefore always involved
the total volume of water at the end of the experiment.
Three different scenarios were used for the data analysis(Tables 3 and 4). Pressure head measurements of all
five tensiometers were included in the objective function
for the first scenario. Only one tensiometer data set at
a time was used for the second scenario. Readings fromtensiometers located closest to the upper and lower
boundaries (Tensiometers 1 and 5) were included in theobjective function for the third scenario. In this way we
hoped to find the minimum information needed forsuccessful numerical inversion. The resulting optimized
soil hydraulic parameters are listed in Tables 3 and 4
for both experiments, together with the final values ofthe objective function <£ and the R2 for regression ofpredicted vs. measured values. Soil hydraulic character-

RESULTS AND DISCUSSION
The laboratory experiments were first analyzed using
the modified Wind method as described by Wendroth

Table 3. Results of parameter estimation and Wind's method for Experiment I.
Scenario

Run

<I>

e,

es

a

n
1

1
2

3

Wind's method

All tensiometers
Tensiometer 1
Tensiometer 2
Tensiometer 3
Tensiometer 4
Tensiometer 5
Tensiometers 1 and 5

R2

K,
1

0.321

cm"
0.02740.02840.03170.02100.02290.02790.0277

1.22
1.20
1.22
1.55
1.49
1.23
1.21

and'
93.1
181.
45.1
24.8
33.1
131.
105.

0.99870.99890.99950.99940.99940.99900.99880.321

0.0249

1.23

73.3

0.992

0.00318
0.002500.001270.001500.001540.00250
0.002880.00550.0001
0.02280.13650.13170.02100.0000

0.321
0.321
0.322
0.307
0.309
0.321

-

0.0045

902

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

Table 4. Results of parameter estimation and Wind's method for Experiment II.

e,

O

Run

Scenario

e,

n

a

and'
9.69
13.2
2.99
4.09

cm-'

0.004370.004120.003410.003430.00304
0.003130.00445-

All tensiometers
Tensiometer 1
Tensiometer 2
Tensiometer 3
Tensiometer 4
Tensiometer 5
Tensiometers 1 and 5

1
2

3Wind's method

0.136

0.134
0.142
0.150
0.154
0.157

0.133
0.0796

istics obtained by numerical inversion and by Wind's
method are shown in Fig. 5 and 7 for both experimental runs.
Notice a very good correspondence between retention curves obtained by all optimization scenarios andthe 0(/z) data points determined using Wind's methodor their analytical fit for Experiment I (Fig. 5a). The
soil water retention parameters obtained with Wind's
method and by optimization of the pressure head readings at all tensiometers, Tensiometer 1 or 5, and both
Tensiometers 1 and 5 are almost identical and indistinguishable from each other in Fig. 5a. A close match of
all retention curves for Experiment I was obtained eventhough different hydraulic parameters were obtained
for numerical inversions using pressure head readings
from Tensiometers 3 or 4. Optimizations using these twotensiometers resulted in substantially higher predictedvalues of both n and 6r, and slightly lower values of 6S.
This result illustrates that a close match or predictionof the soil hydraulic properties within the experimental
range (-20 to -700 cm) does not guarantee accurate
estimation outside the range of measurements, espe-

0.004290.004320.004350.004470.00451
0.004530.004290.00380

0.399
0.399
0.393
0.391
0.389
0.388
0.399
0.395

R1

K,

2.90

2.84
3.07

3.15
3.21
3.27
2.85
2.52

1

0.99930.99930.99950.99950.99950.99950.9993
0.963

5.00
5.55

8.76
2.03

cially toward the higher tensions. This is due to a low
sensitivity of the parameter optimization technique to
0r. If, however, an independent estimate of 6r is available, the estimated retention curve can be extrapolatedwith much more certainty beyond the experimental
range (Parker et al., 1985). As an estimate for 6r one
could either use the true value of the residual watercontent directly (in which case 6r is no longer an optimized parameter) or could include 6r into the objectivefunction in a Bayesian sense as explained above. Although this approach can only worsen the model fit,
incorporation of prior information into the inverse problem will reduce parameter uncertainty (Yeh, 1986). The
final fit cannot be improved by using prior information
since this knowledge represents an additional constraint
on the inverse solution. This lack of improvement is,however, balanced by more stable and reliable parameter estimates. By comparison, the value of 6S can bepredicted better due to knowledge of the final average
o
-100

0.35

a)
r-

-200
?-300

0.3

• Measured 1- Fitted 1* Measured 2- - - Fitted 2
x Measured 3......Fitted 3+ Measured 4- - - - Fitted 4- Measured 5- - - - Fitted 5u•= -400

00.25

All tensiometers
Tensiometer 1-Tensiometer2- Tensiometer 3
Tensiometer 4
Tensiometer 5
Fitted Wind's Points
Wind Method

0.2

-500
-600
-700

2
0.5

1.0

1.5

2.0

2.5

2.5

3.5

4.5

Time [d]3.0log h [cm]-900
-1 = 0.043d- - - t = 0.334 d......t= 1.674 d- - - - 1 = 2.667 d--•-(- 3.507 d

-800
-700

-1 = 4.167 d

-600

& -1

0.0

1.0

1.5

+

Measured-500———All Tensiometers— — — Tensiometer 1—•---•Tensiometer2— • - - Tensiometer 3— - - — Tensiometer 4——— Tensiometer 5———Fitted Wind's Points+ Wind Method• Tension Disc Infiltration-400 ——'———'——r.._.••*•""

-300
-200
-100
2.0

2.5

3.0

log h [cm]Fig. 5. (a) Water retention curves and (b) hydraulic conductivity functions determined with inverse parameter estimation and Wind's
method for Experiment I.

0

-10

-8

-

6

-

4

-

2

0

z[cm]

Fig. 6. Measured and fitted tensiometer readings for Experiment I asa function of (a) time and (b) depth.

903

SlMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION

water content and the evaporation loss during the entire experiment.
We obtained relatively large differences in terms of
the unsaturated hydraulic conductivity (Fig. 5b). The
estimated unsaturated hydraulic conductivity functions
are, again, almost identical for optimizations using all
tensiometer readings and for Wind's method. Also, thehydraulic conductivity function resulting from fitting thepressure heads from Tensiometers 1, 3, 4, and 5 follows
the same pattern of data points obtained by Wind's
method. However, there is now more scatter in the K(h)data points determined by Wind's method; differences
are about a half order of magnitude. Unsaturated hydraulic conductivities obtained with the tension infiltrometer disk on the same sample for pressure heads of—1, -5, and -10 cm are displayed as black squares inFig. 5b. These values were not used in the parameter
optimization, although they could have been easily included in the objective function, thus forcing the unsaturated hydraulic conductivity function to be located close
to these points.
Measured and fitted tensiometer readings for Experiment I are shown in Fig. 6. The highest deviations wereabout 5 and 20 cm for the first and second evaporation
rates, respectively, with the majority of values being
much lower. The model was unable to accurately predict
redistribution during the interruption between the firstand second evaporation rates (Fig. 6a). This could havebeen caused in part by hysteresis effects since, during
redistribution, the upper part of the soil profile is getting
wetter. Hysteresis was not considered in our direct simulations.
The resulting soil hydraulic characteristics for Experi-

ment II are shown in Fig. 7, and their hydraulic parameters are given in Table 4. As for Experiment I, a close
match was obtained for all estimated retention functions
and data points derived with Wind's method. All inversion scenarios resulted in very similar soil hydraulicparameters. Although the retention curves did have similar shapes within the experimental range, the soil hydraulic parameters differed slightly from those obtainedby fitting 0(/z) data points determined by Wind's
method. The values of 6r and n were overestimated byabout 0.05 and 0.4, respectively. These differences leadto different retention values at the higher tensions beyond the experimental pressure head range (greaterthan -700 cm). The differences in unsaturated hydraulicconductivities as determined by different optimization
scenarios or methods differed again by about half anorder of magnitude, which appear acceptable for mostpractical applications. Figure 7b displays unsaturated
hydraulic conductivities obtained with a tension disk
infiltrometer as black squares. These three values lie
well within the range of values covered by the optimized
hydraulic conductivity functions and the fitted Wind
conductivity data. The three infiltrometer data points
were not considered in inverse solution, but used together with the Wind data points to calculate the fitted
unsaturated hydraulic conductivity function. The threevalues are actually somewhat lower than the Wind datafor pressure heads up to -75 cm. Hence, excluding thetension infiltrometer data would have resulted in an
unsaturated hydraulic function that would have corre-

a)
-100

I

-200

Vi

0.45

a)
0.4

\mxfTCtaf-300' 0.35

0.25
0.2

1.5

• Measured 1- Fitted 1* Measured 2-400-All Tensiometers- — — Tensiometer 1--••-•Tensiometer2- - - - Tensiometer 3- • • - Tensiometer 4- — Tensiometer 5-Fitted Wind Points+ Wind Method

0.3

x Measured 3------FittedS+ Measured 4- - - - Fitted 4- Measured 5- - - - Fitted 5

-500

-600

10
Time [d]2.0log h [cm]-900——— = 0.042 d•-- =0.25d- - - - - = 1.257 d• - - • = 5.502 d• - - - = 8.008 d——— = 9.494 d+ Measured

-800
-700
-600
-500

———All tensiometers— — — Tensiometer 1— - - - - - Tensiometer 2— - - - Tensiometer 3— • - - Tensiometer 4——— Tensiometer 5———Fitted Wind's Points+ Wind Method• Tension Disc Infiltration

0.0

0.5

1.0

1.5

-400
-300
-200
-100

2.0

• 2.5

3.0

log h [cm]Fig. 7. (a) Water retention curves and (b) hydraulic conductivity functions determined with inverse parameter estimation and Wind's
method for Experiment II.

0
-10

-8

-

6

-

4

-2

z [cm]
Fig. 8. Measured and fitted tensiometer readings for Experiment II
as a function of (a) time and (b) depth.

904

SOIL SCI. SOC. AM. J., VOL. 62, JULY-AUGUST 1998

sponded more closely with the inverse solution. Finally,in Fig. 8 we show an excellent match of the inverse
solution with the tensiometer readings as a function oftime and depth.
Analysis of both experiments showed that the retention curves within the measurement range can be obtained with far less information required than thatneeded for Wind's method. Moreover, retention curves
obtained using pressure head readings from only one
tensiometer were always very similar to those determined when all information was used, or when Wind's
method was employed. Although the unsaturated hydraulic conductivity estimates showed much more uncertainty, this uncertainty was of the same order of magnitude as the scatter in K(h) data points determinedwith Wind's method, which in turn is probably causedby small-scale soil heterogeneities between the different
compartments, and perhaps some variability in the tensiometer readings.

SUMMARY AND CONCLUSIONSThe evaporation experiment is a widely used experimental technique for simultaneous estimation of soil
hydraulic functions, especially in European laboratories(Halbertsma and Veerman, 1994; Halbertsma, 1996;
Wendroth et al., 1993; Ciollaro and Romano, 1995; Santini et al., 1995). The experiment is usually analyzed bymeans of a modified Wind (1968) method; the perception exists that parameter estimation techniques, in spiteof their successful application in several studies (Ciollaro and Romano, 1995; Santini et al., 1995), are unsuitable for its analysis (Halbertsma, 1996).
In this study we showed that the parameter estimation
technique can be successfully used to analyze an evaporation experiment. Although our study was carried out
on three specific examples (one numerical and two experimental), we discussed a relatively wide range of soils
with n values ranging from 1.23 to 2.52 and a values
from 0.0038 to 0.0249 cm"1. Although some uncertainty
exists in the soil hydraulic parameters caused by highcorrelation between 0r and n, both the retention curve
and the unsaturated hydraulic conductivity function
were predicted well within the range of the Wind-type
measurements. Extrapolation beyond the range of measurement is associated with a high level of uncertainty.
Inclusion of independently measured information beyond the measurement range, i.e., water contents atsome high suction, or directly the residual water content,
could greatly decrease this uncertainty.
Evaporation experiments are usually carried out withfour to five tensiometers (Halbertsma and Veerman,
1994). Methods have been proposed to decrease thenumber of tensiometers to two (Becher, 1971; Schindler,
1980; Wendroth et al., 1993). Using our parameter estimation technique on both numerically generated andlaboratory data, we showed that for the existing experimental setup (five tensiometers in a 10-cm-long sample)
one tensiometer reading is already sufficient to guarantee precise estimation of the soil hydraulic characteristics within the range of measurements. We also showedthat measurements taken closer to the soil surface exhibit a higher sensitivity to the optimized parametersthan measurements from deeper locations, at least for
the relatively short measurement time scale of our study.
While a two-rate evaporation experiment has important advantages over a one-rate experiment in termsof determining the unsaturated hydraulic conductivity
using the modified Wind method, its use in combination
with the parameter estimation technique did not show
any benefits compared with the traditional one-rate approach, except for accelerating the experiment by several days.

SIMUNEK ET AL.: PARAMETER ESTIMATION ANALYSIS FOR SOIL HYDRAULIC PROPERTY DETERMINATION

905