WATER RESOURCES RESEARCH, VOL. 32, NO. 9, PAGES 2683–2696, SEPTEMBER 1996

Estimating unsaturated soil hydraulic properties from tension
disc infiltrometer data by numerical inversion
J. Šimůnek and M. T. van Genuchten
U.S. Salinity Laboratory, U.S. Department of Agriculture, Agricultural Research Service, Riverside, California

Abstract. Tension disc infiltrometers are becoming increasingly popular devices for insitu measurement of the unsaturated hydraulic properties of soil. Tension infiltration data
are generally used to evaluate the parameters K s and a in Gardner’s exponential model ofthe unsaturated hydraulic conductivity. Either two measurements using different disc
diameters or measurements with a single disc but using multiple pressure heads are thenused. In this paper we describe a parameter estimation procedure which combines theLevenberg-Marquardt nonlinear parameter optimization method involving weighted leastsquares, with a quasi-three-dimensional numerical model which solves the variably
saturated flow equation. By numerical inversion of Richards’ equation the unknown
parameters in van Genuchten’s model of the unsaturated soil-hydraulic properties are
estimated from observed cumulative infiltration data during transient water flow.
Additional measurements of the pressure head or water content, as well as a penalty
function for constraining the unknown parameters to remain in some feasible region(Bayesian estimation), can be optionally included into the parameter estimation
procedure. The problem of optimal sampling design, that is, selecting the best points inspace and time for making measurements, is addressed by studying the sensitivity of theobjective function to changes in the optimized hydraulic parameters. We calculate
objective functions based on available cumulative infiltration, pressure head, and watercontent measurements and also on several combinations of these data. The behavior of
the objective function in three-dimensional parameter space is evaluated by means of aseries of two-dimensional response surfaces. The utility of the parameter estimation
procedure is demonstrated using numerically generated data. The sensitivity of theprocedure to different initial estimates of the model parameters is also discussed.

Introduction
Reliable application of computer models to field-scale flow
and transport problems demands a commensurate effort inquantifying a large number of model parameters. As increasingly more complicated flow and transport models are being
developed, the accuracy of numerical simulation depends upon
the accuracy with which various model parameters are estimated. Knowledge of the unsaturated soil hydraulic properties
is especially important when numerical models are used tosimulate variably saturated water flow and contaminant transport. Such simulations are generally based on the numericalsolutions of Richards’ equation, which requires information
about the soil water retention, u (h), and unsaturated hydraulic conductivity, K(h), functions involving the water content u,
the hydraulic conductivity K, and the soil-water pressure headh. Accurate measurement of these hydraulic properties is confounded by the extreme spatial heterogeneity of the subsurface
environment. The hydraulic properties frequently also show
significant variations in time because of cultivation or otheragricultural activities, shrink-swell phenomena of fine-textured
soil, the effect of particle dispersion and soil crusting, andchanges in the concentration and ionic composition of the soil
solution [van Genuchten and Šimůnek, 1996].
A variety of laboratory and field methods are available forThis paper is not subject to U.S. copyright. Published in 1996 by theAmerican Geophysical Union.
Paper number 96WR01525.

direct measurement of the hydraulic conductivity, K, or thesoil water diffusivity, D, as a function of the pressure head
and/or the water content [Klute and Dirksen, 1986; Green et al.,1986]. Most laboratory methods are steady state procedures
based on direct inversion of Darcy’s law. Transient methods
generally involve some type of approximation or simplificationof the Richards’ equation. Popular transient methods include
the Bruce and Klute [1956] horizontal infiltration method andvarious modifications thereof, such as the hot air method andthe sorptivity method. Popular field methods include the instantaneous profile method, various unit-gradient-type approaches, sorptivity methods following ponded infiltration, andthe crust method based on steady water flow. While relatively
simple in concept, these direct measurement methods have anumber of limitations that restrict their use in practice. For
example, most methods are very time consuming because of
the need to adhere to relatively strict initial and boundary
conditions. This is especially true for field gravity-drainage
experiments involving medium- and fine-textured soils or layered profiles. Methods requiring repeated steady state flow
situations or other equilibrium conditions are also tedious,
while linearizations and other approximations or interpolations to allow analytic or semianalytic inversions of the flow
equation introduce additional errors. Finally, information
about uncertainty in the estimated hydraulic parameters is notreadily obtained using direct inversion methods.
A more flexible approach for solving the inverse problem isthe use of parameter optimization methods. Optimization pro-

2683

2684

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

cedures also make it possible to simultaneously estimate theretention and hydraulic conductivity functions from transient
flow data [Kool et al., 1987]. Early parameter optimization
studies focused primarily on solute transport [e.g., van Genuchten, 1981; Parker and van Genuchten, 1984]. Having started
with the studies of Zachmann et al. [1981] and Dane andHruska [1983], the method is now increasingly being used also
for estimating the unsaturated soil hydraulic functions. Computer models 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 and Hopmans [1993] for multistep approaches. Considerable attention has also been given to
the estimation of soil-hydraulic properties from ponded infiltration experiments [Russo et al., 1991; Bohne et al., 1992].
While initially applied to laboratory-type experiments, inverse
methods are equally well applicable to field data [Kool andParker, 1988] or some appropriate combination of field andlaboratory data. An important advantage of inverse procedures, if formulated within the context of a parameter optimization problem, is that a detailed error analysis of the estimated parameter is more easily considered. While parameter
optimization methods provide several advantages, a number ofproblems related to computational efficiency, convergence,and parameter uniqueness remain to be solved, especiallywhen many hydraulic parameters must be estimated simultaneously [van Genuchten and Leij, 1992].
Tension disc infiltrometers are increasingly being used for insitu measurement of the unsaturated soil hydraulic properties[e.g., Perroux and White, 1988; Ankeny et al., 1991; Reynolds andElrick, 1991; Logsdon and Jaynes, 1993]. Tension infiltrometers
are especially useful for quantifying the effects of macroporesand preferential flow paths on infiltration in the field. Thus far,tension infiltration data have been used primarily for evaluating the saturated hydraulic conductivity K s and the sorptivity
parameter a in Gardner’s exponential model [Gardner, 1958]of the unsaturated hydraulic conductivity,
K~h! 5 K s exp ~ a h!

(1)

Analyses of this type require either two infiltration measurements using two different disc diameters [Smettem and Clothier, 1989] or measurements using a single disc diameter but
then imposing multiple tensions [Ankeny et al., 1991]. The
parameters K s and a are generally estimated using Wooding’s[1968] equation for unconfined steady state infiltration from a
disc,
Q~h 0! 5 p r 20 K~h 0! 1 4r 0f ~h 0!

(2)where Q is the steady state infiltration rate, r 0 is the radius of
the disc, h 0 is the wetting pressure head, and f is the matrix
flux potential or, equivalently, the linearized Kirchhoff transform given by [Gardner, 1958]f ~h 0, h i! 5

E

h0

K~h! dh(3)hi

in which h i is the initial pressure head. Equation (2) considers
the two main forces driving the infiltration process: the firstterm represents one-dimensional gravity-driven flow, while thesecond term accounts for capillary-induced flow [Wooding,1968]. Explicit in (2) is the assumption that the hydraulic conductivity at the initial pressure head is much lower than thehydraulic conductivity at the wetting pressure head and therefore can be neglected. Since (2) contains two unknown parameters, K and f, or K s and a after using (1) and (3), the equation
has to be written twice for two different experiments as explained above, using two different disc diameters or two different tensions, unless one of the parameter can be estimatedindependently. For this purpose White and Sully [1987] defined
the matrix flux potential in terms of the sorptivity which they
estimated from early-time infiltration data. However, extreme
care is needed in certain cases to properly define the timeinterval during which the relation between the cumulative infiltration and the square root of time is linear, with sorptivity
being the proportionality constant [Quadri et al., 1994]. In adifferent approach, Logsdon and Jaynes [1993] proposed a regression method by which the parameters K s and a were fitted
simultaneously to more than two measurements involving different wetting pressure heads at the same location.
As compared to the approximate analytical approaches
above, relatively little work has been done in simulating unsaturated flow underneath a disc permeameter using more complete numerical solutions of Richards’ equation and even lessso in attempts to estimate the unsaturated soil hydraulic properties (including the soil-water retention curve) from tension
disc infiltration experiments by means of inverse solutions ofRichards’ equation. Quadri et al. [1994] developed a finite
difference numerical model for two-dimensional, axisymmetricwater flow and solute transport underneath a disc permeameter. They compared numerical results with laboratory dataobtained with a 1/4-sector disc permeameter. Good predictions
were obtained for the observed water content and Br 2 profiles
in the sand box as well as for the measured infiltration ratefrom the disc [Quadri et al., 1994]. In an earlier study, Warrick[1992] developed a finite element program, “Disc,” which heused for comparing steady state and transient axisymetric flow
rates underneath a disc permeameter with alternative solutionsof Richards’ equation.
In this paper we describe a parameter estimation procedure
which combines Levenberg-Marquardt nonlinear parameter
optimization involving weighted least squares with a quasithree-dimensional numerical model which solves the variably
saturated flow equation. The unknown parameters in van
Genuchten’s model for the unsaturated soil properties are
estimated from observed cumulative infiltration data duringtransient water flow by numerical inversion of Richards’ equation. Additional data, such as measured pressure heads orwater contents, as well as a penalty function for constraining
the optimized parameters to remain in some feasible region
using Bayesian estimation, can be optionally included in theparameter estimation procedure. Objective functions are formulated in terms of cumulative infiltration, pressure head, andwater content data, and on a combination of these measurements. The behavior of the objective functions in threedimensional parameter space is studied by means of a series of
response surfaces in two-dimensional parameter planes. The
parameter estimation procedure is demonstrated using numerically generated data; the sensitivity of the method to different
initial estimates of the model parameters will be discussed also.
One important issue addressed in this paper is the questionof whether infiltration data usually collected during tension
disc permeameter experiments are sufficient for the parameter
estimation procedure or whether additional data are needed toguarantee a unique inverse solution. We will also examine

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

which additional data should be collected so as to improve theparameter inversion process.

TheoryThe governing flow equation for radially symmetric isothermal Darcian flow in a variably saturated isotropic rigid porous
medium, under the assumption that the air phase plays an
insignificant role in the liquid flow process, is given by thefollowing modified form of Richards’ equation [Warrick, 1992]:­u 1 ­
5
­t
r ­r

S

rK

­h
­r

D

1

­
­z

S D
K

­h
­z

1

­K
­z

(4)

where r is a radial coordinate, z is vertical coordinate positive
upward, and t is time. Equation (4) was solved numerically forthe following initial and boundary equations applicable to adisc tension infiltrometer experiment [Warrick, 1992]:
h~r, z, t! 5 h ih~r, z, t! 5 h 0
2

t50

(5)

0 , r , r 0, z 5 0(6)­h~r, z, t!
2150­zh~r, z, t! 5 h ir . r 0, z 5 0
r2 1 z2 3 `

(7)
(8)

Equation (4), subject to the above initial and boundary conditions, was solved using a quasi-three-dimensional (axisymmetric) finite element scheme as documented by Šimůnek et al.
[1996]. The numerical solution was based on the massconservative iterative scheme proposed by Celia et al. [1990].
A model of the unsaturated soil hydraulic properties must be
selected prior to application of the numerical solution of Richards’ equation. In this study we will limit ourselves to unsaturated soil hydraulic functions of the form [van Genuchten,
1980]

u e~h! 5

u ~h! 2 u r
1
5
us 2 ur~1 1 u a hu n! mu ~h! 5 u s

(9)

h$0

m 2
K~ u ! 5 K su 1/e 2@1 2 ~1 2 u 1/m
e ! #K~ u ! 5 K s

h,0

h$0

h,0

(10)

where u e is the effective water content; u r and u s denote theresidual and saturated water contents, respectively; and a , n,and m (51 2 1/n) are empirical parameters. Note that ahere is different from that used in (1). The hydraulic characteristics defined by (9) and (10) contain five unknown parameters: u r , u s , a , n, and K s . The saturated hydraulic conductivity, K s , is perceived here and in the remainder of the paper
as a fitted hydraulic conductivity at zero pressure head. In
reality, the saturated hydraulic conductivity might be different
from this value owing to the effects of macropores which saturate only after a zero or positive pressure head is applied[Logsdon et al., 1993]. Following Logsdon and Jaynes [1993], weview large macropores as not being predictive of hydraulicproperties of the bulk soil matrix.
In the analysis below we will assume that the residual andsaturated water contents have been measured independently.
The values of the shape parameters a and n will be sought bynumerical inversion of the flow problem. The saturated hydraulic conductivity K s (in the sense of being the hydraulic2685

conductivity at or extrapolated to zero pressure head) will beeither assumed known through independent measurement or
kept as an unknown parameter to be estimated also by numerical inversion of Richards’ equation. Since the tension infiltration experiments by definition do not reach the saturated watercontent, u s , and the saturated hydraulic conductivity, K s , thesetwo parameters are to be interpreted within this study as extrapolated, empirical parameters outside the range of the discexperiment. Also, tension disc infiltration in general is a wetting process (assuming that one can neglect internal drainage
at the initial pressure head); the hydraulic parameters in (9)and (10) should represent wetting branches of the unsaturated
hydraulic properties.
Formulation of the Inverse Problem
Several direct and indirect methods [Neuman, 1973] may be
used to carry out the inverse problem of parameter identification. Direct methods formally treat the model parameters as
dependent variables in an inverse boundary value problem[Yeh, 1986]. By contrast, indirect methods are based upon theminimization of a suitable objective function which expresses
the discrepancy between the observed values and the predicted
system response. Initial estimates of the parameters are theniteratively improved upon during this minimization processuntil a desired precision is obtained. The technique used in thispaper constitutes such an indirect approach.
When measurement errors follow a multivariate normal distribution with zero mean and covariance matrix V, the likelihood function can be written as [Bard, 1974]
L~b! 5 ~2 p ! 2n/ 2 det21/ 2 V? exp $~21/ 2!@q* 2 q~b!# TV 21@q* 2 q~b!#%(11)

where L(b) is the likelihood function, b 5 { b 1 , b 2 , z z z , b m }is the vector of optimized parameters (e.g., u r , u s , a , n, andK s ), m is the number of optimized parameters, q* 5 {q *1 ,q *2 , z z z , q *n } is a vector of observations (such as pressure
heads h, water contents u, and/or cumulative infiltration rates
Q), q(b) 5 {q 1 , q 2 , z z z , q n } is a corresponding vector ofmodel predictions as a function of the unknown parameters
being optimized, and n is the number of observations. The
likelihood function L(b) is defined as the joint probability
density function (pdf) of the observations and is considered afunction of the unknown parameters b [Bard, 1974]. The maximum likelihood estimate is that value of the unknown parameter vector b which maximizes the value of the same likelihood
function. Since the logarithm is a monotonically increasing
function of its argument, the value of b which maximizes thelikelihood function L(b) also maximizes ln L(b) [Bard, 1974].
This property of the logarithm is frequently used in parameter
identification studies since ln L is often a simpler function oreasier to use than L itself. Equation (1) is hence reformulated
asln L~b! 5 2

n1
ln ~2 p ! 2 2 ln (det V)21

2 2 @q* 2 q~b!# TV 21@q* 2 q~b!#(12)

The maximum of the likelihood function must satisfy the set ofb likelihood equations­ ln L~b!
50
­b

(13)

2686

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

If all elements of the covariance matrix V are known, then thevalue of the unknown parameter vector b which maximizes(12) must minimize the following equation:
F~b! 5 @q* 2 q~b!# TV 21@q* 2 q~b!#(14)

If something about the distribution of the fitted parameters
is known before the inversion, that information can be included into the parameter identification procedure by multiplying the likelihood function by the prior pdf, p 0 (b), whichsummarizes the prior information. Estimates which make useof prior information are known as Bayesian estimates and leadto the maximizing of a posterior pdf, p*(b), given by
p*~b! 5 cL~b! p 0~b!

(15)

in which c is a constant. The posterior density function is
proportional to the likelihood function when the prior distribution is uniform. Inclusion of the prior information leads tothe maximization of the following equation for F:
F~b! 5 @q* 2 q~b!# TV 21@q* 2 q~b!#
1 ~b* 2 b̂! TV 21
b ~b* 2 b̂!

(16)where b* is the parameter vector containing the prior information (e.g., u r , u s , a , n, and K s ), b̂ is the associated predicted parameter vector, and Vb is a covariance matrix for theparameter vector b. The second term of (16), sometimes called
the penalty function, insures that the obtained estimate is constrained to some feasible region around b* so as to remain
physically meaningful. Russo et al. [1991] showed that the useof a penalty function can significantly improve the uniquenessof the estimate.
The covariance matrices V and Vb, sometimes also called theweighting matrices, provide information about the measurement accuracy as well as about any possible correlation between measurement errors and between parameters [Kool etal., 1987]. When the covariance matrix V is diagonal and allelements of Vb are equal to zero, that is, the measurement
errors are uncorrelated and no prior information about the
optimized parameters exists, the problem simplifies to theweighted least squares problem

On

F~b! 5w i@q *i 2 q i~b!# 2

(17)

i51

where w i is the weight of a particular measured point.
Solution of the Inverse Problem
Many techniques have been developed in the past to solve
the nonlinear minimization/maximization problem [Bard,
1974; Beck and Arnold, 1977; Yeh, 1986; Kool et al., 1987]. Mostof these methods are iterative by first starting with an initial
estimate bi of the unknown parameters to be estimated andthen studying how the objective function F(b) behaves in thevicinity of the initial estimate. On the basis of this behavior,
one selects a direction vector vi such that the new values of theunknown parameter vector,
b i11 5 b i 1 r iv i 5 b i 2 r iR ip i

(18)

decreases the value of the objective function
F i11 , F i

where F i11 and F i are the objective functions at the previous
and current iteration level, Ri is a positive definite matrix, andr i is a scalar that insures that the iteration step is acceptable.
Methods based on (18) are called gradient methods. Differences among the various gradient methods presented in theliterature (e.g., steepest descent, Newton’s method, directional
discrimination, Marquardt’s method, the Gauss method, variable metric methods, and the interpolation-extrapolation
method) are a result of differences in choosing the step direction vi and/or the step size r i [Bard, 1974]. Steepest descentmethod selects r i 5 1 and Ri 5 I, where I is an identity matrix,
whereas Newton’s method uses r i 5 1 and Ri 5 H21, wherei
H is the Hessian matrix of F(b):(19)

H ij~b! 5­ 2F­ b i­ b j(20)

The Gauss-Newton method simply neglects the higher-orderderivatives in the definition of the Hessian matrix and assumes
that H can be approximated by a matrix N using only the
first-order derivatives. For nonlinear weighted least squares
this leads to
H < N 5 J TJ

(21)

in which J is the Jacobian matrix whose elements are given by
the sensitivity coefficients multiplied by a square root of theweight of a particular data point as follows:
J ij 5­@q *i 2 q i~b!#w 1/i 2­q i~b!
5 2w 1/i 2
­bj
­bj

(22)

Marquardt [1963] proposed a very effective method which
has become a standard in nonlinear least squares fitting. The
method represents a compromise between the inverse Hessianand steepest descend methods by using the steepest descentmethod when the objective function is far from its minimum
and switching to the inverse Hessian method as the minimum
is approached. This switch is accomplished by multiplying the
diagonal in the Hessian matrix (or its approximation N), sometimes called the curvature matrix, with (1 1 l), where l is apositive scalar. When l is large, then the matrix is diagonallydominant resulting in the steepest descend method. On theother hand, when l is zero, the inverse Hessian method will
result.
The sensitivity coefficients may be calculated using threedifferent methods [Yeh, 1986]: the influence coefficient method(finite differences), a sensitivity equation method, and a variational method. When using the influence method, one has tobalance truncation errors which increase with Db against
rounding errors in differencing which decrease with Db. A
common practice is to change the optimized parameters by one
percent to obtain estimates of the sensitivity coefficients,­q i q i~b 1 Dbe j! 2 q i~b!
5
­bj
Dbj

(23)

where e j is the jth unit vector, and Db 5 0.01 b.

Data Generation
Infiltration data describing flow from a tension disc werenumerically generated for three hypothetical soils (sand, loam,
and clay [Hillel and van Bavel, 1976]) and used to evaluate thefeasibility of the proposed method for estimating unsaturated

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Figure 1. Cumulative infiltration versus time for a hypothetical disc permeameter infiltration experiment.

soil hydraulic parameters. Only results for the loam will bereported here since the results and conclusions obtained forthis soil were qualitatively very much the same as those for the
hypothetical sand and clay soils.
The soil hydraulic parameters of the hypothetical loam weretaken as follows: u r 5 0.000, u s 5 0.480, a 5 0.0150 cm21,
n 5 1.592, and K s 5 0.0007 cm s21. The radius of the discpermeameter was assumed to be 10 cm. The initial pressure
head of the homogeneous and isotropic soil beneath the disc
was taken as 2500 cm, whereas the pressure head in the discwas set at 23 cm. Only measurements during the first 12 hoursof the experiment were included in the analysis below. The
ratio of the infiltration rate at 12 hours to those at 24 and 48
hours was 1.021 and 1.039, respectively. This result should
signal a caution to field workers since steady state infiltration
was apparently not reached during the first 12 hours. Still, therelatively small changes in infiltration rate in this example after
12 hours are probably negligible from a practical point of view.
We will follow an analysis somewhat different than that usedpreviously by Kool and Parker [1988], who considered a hypothetical one-dimensional ponded infiltration experiment fol-

2687

Figure 3. Calculated pressure head versus time at two observation nodes, (r, z) 5 (10, 210) and (20, 210), for the
hypothetical disc permeameter infiltration experiment.

lowed by redistribution and evapotranspiration. Our analysis
involves the following steps and calculations.
1. The known hypothetical soil-hydraulic properties u (h)and K(h) were first used in the direct problem to simulate thetension disc infiltration event for given initial and boundary
conditions. The finite element discretization was selected sothat the mass balance error for the direct solutions always
remained less than 0.05%.
2. The instantaneous and cumulative (Figure 1) infiltration
rates that were generated were discretized at given times,
whereas the calculated water content (Figure 2) and pressure
head (Figure 3) profiles were sampled at two given locations,again at discrete times. Both nodes were at 10 cm depth, one
at the edge of the disc, the other at twice the radius. In thisstudy we used local nodal values of the pressure head andwater content for the hypothetical observations. We emphasize, however, that experimental devices for measuring h and u,such as tensiometers and time domain reflectometry, havealways some finite length rather than yielding point measurements, that is, they average the measured variable over a finite
volume of soil:ū ~r, z, t! 5

1
V

E
E

u ~r, z, t! dVVh̄~r, z, t! 5

1
V

(24)
h~r, z, t! dVV

Figure 2. Calculated water content versus time at two observation nodes, (r, z) 5 (10, 210) and (20, 210), for the
hypothetical disc permeameter infiltration experiment.

where ū and h̄ are the spatial averages of the water content andpressure head, respectively, over a volume V representing the
averaging volume of the particular measurement device. Integrating pressure heads or water contents over a given volume
around a specified location for a particular measurement device should not pose any significant numerical problem as such.
However, interesting questions could be raised, such as in what
direction the measurement devices should be placed (horizontally or vertically) to obtain the greatest benefit in terms ofparameter identification [Vogeler et al., 1996], or whether thesimultaneous measurement of both horizontally and verticallyaveraged pressure heads or water contents could provide ad-

2688

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Experimental Design
In general, experiments should be designed such that measurements are made which yield the most information about
the unknown parameters to be optimized, that is, measurements which are most sensitive to changes in the unknown
parameters. The sensitivity of the objective function to the
optimized parameters is usually studied in terms of sensitivity
coefficients as given by (23). Unfortunately, sensitivity coefficients calculated in this manner are difficult to compare witheach other since they depend on the absolute values of thehydraulic parameters involved as well as on the invoked units.
Therefore we will evaluate the sensitivities in a different way so
as to remove their dependency on parameter units and absolute values:
Figure 4. Sensitivity of the instantaneous infiltration rate,v(t), to a 1% change in hydraulic parameters a , n, u s , and K s .

e ij 5 100 b j

U ­­qb U 5 100 b uq ~b 11.01Dbeb 2! 2b q ~b!uii

5 uq i~b 1 Dbe j! 2 q i~b!uditional information for improving parameter identifiability (B.
E. Clothier, personal communication, 1996). Since these questions go beyond the scope of this paper, they will not be further
addressed here. Infiltration rates were calculated as the sum of
the actual nodal fluxes Q i associated with nodes having prescribed Dirichlet boundary conditions representing the discpermeameter. The nodal fluxes could be calculated explicitly
and accurately from the original finite element equations associated with these nodes [Šimůnek et al., 1996]. The cumulative disc permeameter infiltration rate, Q(t), was calculated asfollows:
Q~t! 5

E

t

t0

EO
t

v~t! dt 5n

Q i~t! dt(25)

t0 i51

where t 0 is the starting time of the experiment, v(t) is theinstantaneous infiltration rate, and n is the number of nodes
representing the disc permeameter.
3. The data obtained in steps 1 and 2 were used to calculate sensitivities of the objective function to the optimized
parameters. A similar study on behavior of sensitivities in the
one-dimensional advection-dispersion equation was carriedout by Knopman and Voss [1987].
4. The data obtained in steps 1 and 2 were next used to
calculate response surfaces of the objective function as a function of a particular hydraulic parameter or a combination ofparameters so as to determine possible uniqueness problems of
the inverse procedure.
5. The data obtained in steps 1 and 2 were used as input
data for the inverse problem. In a first group of inversions,instantaneous infiltration rates, cumulative infiltration data, oradditional “data” (water contents and pressure heads measured at a given location) were used individually, without any
additional information. Combinations of these primary measurement sets were subsequently used in a second group of
inversions to evaluate whether additional information would
substantially improve the error estimates on the parameters aswell as on the goodness of fit. Sensitivities to the initial parameter estimates were also studied. Sensitivity coefficients in thisstudy were calculated using the Levenberg-Marquardt method,combined with the influence coefficient method.

j

i

j

j

j

j

(26)

where e ij is the change in the auxiliary variable q i corresponding to a 1% change in parameter b j .
Figure 4 shows the sensitivity of the instantaneous infiltration rate, v(t), to 1% changes in the parameters a , n, u s , andK s . Since the residual water content was assumed to be zero,
we did not use u r in our sensitivity study. Because of physicalconsiderations, we expected that there would be a certain timeinterval (0, t grav) at the beginning of the infiltration event
during which the infiltration rate would be very sensitive to theunknown parameters. The characteristic time t grav would define the point in time when gravity starts dominating the infiltration process [Warrick, 1992]. Figure 4 shows that the infiltration rate is very sensitive to changes in all four parameters
at only very short times, that is, less than 10 s. After that, thesensitivity of the infiltration rate becomes essentially constantand almost independent of time. Since it is practically impossible to obtain several measurements within such a short time
period immediately after initiating the infiltration event, v is
not a recommended variable for use in the inverse problem
and hence will not be further considered here for calculating
response surfaces.
The sensitivities of the cumulative infiltration, Q(t), to 1%changes in the hydraulic parameters is presented in Figure 5.
Notice that the sensitivities increase with time for all four

Figure 5. Sensitivity of the cumulative infiltration rate, Q(t),to a 1% change in hydraulic parameters a , n, u s , and K s .

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

2689

parameters. The cumulative infiltration rate is most sensitive tothe parameter n and least to the saturated water content u s .
The sensitivity to n is about three times higher than the sensitivities to a and the saturated hydraulic conductivity K s ,which are about the same. If the absolute values in (26) were
to be removed, the sensitivity to a would get a sign opposite tothat of K s . When comparing the sensitivities to the different
parameters, it is important to recognize that a 1% change inparameter n, which can change in order of units and therefore
1% change covers larger range of possible values, will be much
more important than a 1% change in K s , which can vary several orders of magnitude in the field. The results presented in
Figure 5 suggest that the further in time the measurements aretaken the more valuable the information becomes, with thelimiting situation occurring when steady state has beenreached.
Figure 6 presents the sensitivity of the pressure head, h, tothe hydraulic parameters a , n, u s , and K s at observation nodes(r, z) located at (10, 210) (Figure 6a) and (20, 210) (Figure6b), that is, one node slightly slanted below the disc and one
right below the edge of the permeameter. The highest sensitivities of the pressure head at both observation nodes occur
during the time period when the moisture front passes throughthese points. As with the cumulative infiltration rate, 1%changes in n induce the biggest changes in the pressure head,
although the sensitivities to the other three parameters are notmuch smaller. The sensitivities of h are about the same at thetwo observation nodes, although the second node located further to the side of the source maintains its high sensitivity for
Figure 7. Sensitivity of the water content at two observation
nodes (a) (r, z) 5 (10, 210) and (b) (r, z) 5 (20, 210) to
a 1% change in hydraulic parameters a , n, u s , and K s .

Figure 6. Sensitivity of the pressure head at two observation
nodes (a) (r, z) 5 (10, 210) and (b) (r, z) 5 (20, 210) to
a 1% change in hydraulic parameters a , n, u s , and K s .

a much longer time period because of a more diffuse moisture
front. Hence, as compared to measurement points below thesource, locations at the side of the disc should yield more
information for the inverse problem, but this will be at theexpense of a prolonged experimental effort and more computer time. In practice, a compromise between these two locations may be needed. Still, whatever the location of measurement, data of the pressure head should be taken when the
moisture front passes through that location.
Sensitivities of the water content to a , n, u s , and K s at thesame two observation nodes, (r, z) 5 (10, 210) and (20,210), are shown in Figure 7. The initial sensitivities of u before
arrival of the moisture front are caused by the initial condition
being prescribed in terms of the pressure head. This condition
causes the water content at early times to be sensitive to allparameters used in the description of the retention curve ( a ,n, u s ) but not to the saturated hydraulic conductivity, K s . The
water content during passage of the moisture front is very
sensitive to a and K s but not very sensitive to n and u s . After
passage of the moisture front the sensitivity to a and K s decreases gradually, although more slowly than for the pressure
head (Figure 6). By contrast, the sensitivity to the saturated
water content u s increases after passage of the moisture front;
this is to be expected since the soil water content is now slowly
approaching (and hence defining) the value of u s . Since thewater content remains sensitive to at least one hydraulic parameter during the entire infiltration event, water content mea-

2690

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Table 1. Grid Spacings Used for the Parameter Planes ofthe Hypothetical Disc Permeameter Infiltration ExperimentParameter

Lower
Parameter
Value

Parameter
Step Value

Upper
ParameterValuea, cm21n [2]K s , cm s21

0.002
1.03330.00001

0.002
0.03333
0.000050.06
2.0
0.0015surements should cover the entirety of the infiltration experiment.
The results presented in Figures 4 through 7 can be interpreted in two different ways. First, observed data should betaken at points in time and space which show the highest
sensitivities of those observations to the hydraulic parameters
and which are easily accessible to measure. The inverse procedure will not improve by including in the objective function
data which are insensitive to the parameters to be optimized.
For example, including pressure head data measured before
the arrival or after passage of the moisture front will not
improve the identifiability of the hydraulic parameters since
they will lead to similar predicted pressure heads during thosetimes. Second, higher weights should be placed on “measurements” with show the highest sensitivity to the optimized parameters, and relative small weights should be placed on other,more insensitive observations.

Response Surfaces
Russo et al. [1991] reported for a one-dimensional ponded
infiltration experiment that the cumulative infiltration curve
does not provide enough information to obtain a unique set of
optimized parameters. Similar results were also obtained forone-step outflow experiments when only the outflow curve wasused in the inverse problem [van Dam et al., 1992, 1994; Echingand Hopmans, 1993; Toorman et al., 1992] except when relativelarge pressures were imposed on the samples so as to provide
a high resolution in the average water content of the sample[Parker et al., 1985]. Additional information of the pressure
head inside of the sample provided enough information to
successfully apply the inverse problem by yielding a unique setof optimized parameters. In this section we test whether theconventionally measured data on cumulative infiltration ratefrom a disc permeameter provide enough information to enable the identification of a unique set of soil hydraulic parameters from the inverse problem and, if not, what additional
information is necessary to successfully carry out the inversion.
We will approach this question in a way similar to that ofToorman et al. [1992] in their study of the uniqueness of theinverse problem for one-dimensional one-step outflow experiments.
The uniqueness of the inverse problem will be evaluated interms of two-dimensional response surfaces of the objectivefunction as a function of pairs of soil-hydraulic parameters.
The objective function used for this purpose is given as

OS OmF~b, q m! 5

nj

wj

j51

i51

w i@q *j~t i! 2 q j~t i, b!# 2

D

(27)

where m represents the different sets of measurements, that is,
those involving cumulative infiltration, pressure head, or watercontent data; n j is the number of measurements in a particular
set; q *j (t i ) are specific measurements at time t i for the jthmeasurement set; q j (t i , b) are the corresponding model predictions for the parameter vector b; and w j and w i are weights
associated with a particular measurement set or point, respectively. We assumed that the weighting coefficients w i in (27)are equal to 1, that is, the variances of the errors inside of aparticular measurement set are the same. The weighting coefficient for the water content set as a whole, w u , was also set
equal to 1, while those for the cumulative infiltration, w Q , andthe pressure head, w h , sets were defined as

3 43 4O O
3 43 4
O
nu

wh 5

O
nh

q *uj

j51

q *hi

i51

nu

nu

nQ

j51

nu

(28)

nh

q *uj

wQ 5

21

21

q *Qi

i51

nQ

(29)

respectively, where n u , n h , and n Q represent the total numberof measurements obtained for the water content, pressure
head, and cumulative infiltration, respectively; and q *u , q *h , andq *Q are the actual measurements of u , h, and Q, respectively.
We calculated the objective functions for three parameter
planes (i.e., a -n, a -K s , and n-K s ), for three separate sets ofmeasurements (i.e., for cumulative infiltration, pressure head,
and water content data), and for four combined sets of data
(i.e., Q 1 h 1 u , Q 1 h, Q 1 u , and h 1 u ). The response
surfaces were calculated on a rectangular grid with parametervalues given in Table 1. Each parameter domain was discretized into 30 discrete points, resulting in 900 grid points foreach response surface. A total of more than 60 response surfaces was calculated in this manner (three soil types, although
only the loam results are presented here, and three twodimensional parameter planes for each of three types of separate or four combined data measurements).
Figure 8 represents response surfaces of the objective function F(Q) for the cumulative infiltration measurements (seealso Figure 1) for the three different parameter planes, a -n,
a -K s , and n-K s . The a -n response surface (Figure 8a) showsa well-defined valley which starts at low a and n values andextends linearly through nearly the entire parameter space.
Notice that we decreased the distance between the contour
levels toward the low values of the objective function so as tomore precisely locate the minimum of F(Q). We obtained
visually almost identical response surface in the parameter
plane a -K s (Fig. 8b), although now several local minima appeared for relatively small values of the parameters a and K s .
Figures 8a and 8b indicate that a higher value of a and corresponding increases in n or K s will lead to very similar values ofthe objective function, thus indicating possible difficulties infinding a unique inverse solution. Figure 8c also displays awell-defined narrow valley, now of hyperbolic shape, suggesting an inverse relationship between the parameters n and K s interms of their effect on the objective function. An increase inn and a corresponding decrease in K s will produce the same
response in the objective function. This hyperbolic behavior of
F(Q) indicates that the saturated hydraulic conductivity ismore identifiable when n is relatively large (e.g., greater thanabout 1.4), while n is similarly better identifiable when K s is

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

2691

Figure 8. Contours of the objective function F(Q) for the cumulative infiltration rate in the (a) a -n plane,(b) a -K s plane, and (c) n-K s parameter plane.

large (e.g., .0.008 cm/s). When n is larger than about 1.4, asmall change in the saturated hydraulic conductivity results ina relatively large change in the objective function. The same istrue for the parameter n when the saturated hydraulic conductivity is large. The nonexistence of a well-defined minimum inthe response surfaces in Figure 8 permits a broadening of theconclusion of Russo et al. [1991] that identical infiltration
curves for one-dimensional ponded infiltration can be generated by an infinite number of combinations of parameters a ,n, and K s . His conclusion holds also in three dimensional
space for tension disc permeameters.
Figure 9 represents response surfaces of F(h) for the pressure head measured at (r, z) 5 (20, 210) (Figure 3), again
in the three parameter planes a -n, a -K s , and n-K s . The
response surfaces for the h-based objective function, F(h),are visually quite similar to those for the cumulative infiltration
rate, F(Q). An increase in a and a simultaneous increases in
K s or n again lead to approximately the same values for theobjective function. The contours of F(h) and the valleyaround its minimum in Figure 9a are almost parallel with then axis, thus indicating that the parameter n is very difficult to
identify from available information. And similarly, as F(Q) in
Figure 8c, the pressure head based objective function F(h)reacts inversely to an increase in the parameters n and K s(Figure 9c). The contours of F(h) at relatively low values of
K s , as well as the long valley around the minimum, are almost
parallel with the n axis (Figure 9c), thus showing that n can not
be identified precisely when all three parameters a , n, and K s
are to be estimated simultaneously. Still, the minimum of F(h)is somewhat better identifiable than the minimum of F(Q).
Figure 10 shows similar response surfaces for the objectivefunction F(u) for the water content measured at (r, z) 5
(20, 210) (Figure 2). The objective function F(u) has very
well defined minima in all three parameter planes. Increased
identifiability of the parameters is partly caused by the fact thatthe initial condition was specified in terms of the pressure
head. Adding one water content measurement will then cause
one point of the retention curve to become known.
The response surfaces shown thus far are for objective functions associated with only one measured variable: Q, h, or u.

Figure 9. Contours of the objective function F(h) for the pressure head measured at (r, z) 5 (20, 210)
in the (a) a -n plane, (b) a -K s plane, and (c) n-K s parameter plane.

2692

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Figure 10. Contours of objective function F( u ) for the water content measured at (r, z) 5 (20, 210) inthe (a) a -n plane, (b) a -K s plane, and (c) n-K s parameter plane.

We now address the question whether the identifiability of theparameters will increase when the objective functions for two
variables are combined. Figure 11 presents in all three parameter planes the contours of the objective function F(Q, h)(5F(Q) 1 F(h)) defined by both the cumulative infiltrationand the pressure head measured at (r, z) 5 (20, 210). One
can readily see that the identifiability of the parameters does
not markedly increase in comparison with the response surfaces of the separate objective functions F(Q) and F(h). This
result is mainly due to the similarity of the response surfaces
when the two variables are considered separately. Figures 8
and 9 show that the valleys surrounding the minima in bothcases have almost an identical shape; hence simply adding thetwo objective functions will not improve the definition of theminima. A more optimal inverse scenario would be if thevalleys of the response surfaces surrounding the minima wereperpendicular for the two variables, thereby improving thedefinition of the minima when the objective functions are combined. This is not the case here where combination of theobjective functions for the pressure head and the cumulative

infiltration actually resulted in several local minima in theparameter planes a -n and a -K s (Figures 11a and 11b).
Better results are obtained when the objective functions
based on the cumulative infiltration, F(Q), and the watercontent measured at (r, z) 5 (20, 210), F( u ), are combined. Figure 12 shows that the response surfaces for thisobjective function F(Q, u ) are now well defined in all threeplanes ( a -n, a -K s , and n-K s ). An exception is perhaps thea -K s plane, which shows a valley that extends over a relativewide range of a and K s values. Finally, and for completeness,
we also present in Figures 13 and 14 contours of objectivefunctions F(h, u ) and F(Q, h, u ), respectively, in the threeparameter planes. The minima of the objective functions in all
the planes are again very well identifiable.

Inverse Solutions
To evaluate the robustness of the parameter optimization
procedure and the dependence of the final solution on theinitial estimates of the unknown hydraulic parameters, we car-Figure 11. Contours of the objective function F(Q, h) for the cumulative infiltration rate and the pressure
head measured at (r, z) 5 (20, 210) in the (a) a -n plane, (b) a -K s plane, and (c) n-K s parameter plane.

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Figure 12. Contours of the objective function F(Q, u ) for the cumulative infiltration rate and the watercontent measured at (r, z) 5 (20, 210) in the (a) a -n plane, (b) a -K s plane, and (c) n-K s parameter plane.

Figure 13. Contours of the objective function F(h, u ) for the water content and the pressure headmeasured at (r, z) 5 (20, 210) in the (a) a -n plane, (b) a -K s plane, and (c) n-K s parameter plane.

Figure 14. Contours of the objective functions F(Q, h, u ) for the cumulative infiltration, the water content
and the pressure head measured at (r, z) 5 (20, 210) in the (a) a -n plane, (b) a -K s plane, and (c) n-K s
parameter plane.

2693

2694

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Table 2.

Results Obtained With the Different Parameter Estimation Computer Runs

Case

FittedParameter

Initial Estimates

1a*a , n, K sb*a , n, K sc*a , n, K s
2a
a, n
b
a, n
c
a, n
3a*
a , n, K sb*a , n, K sc*a , n, K s
4a
a, n
b*
a, n
c
a, n
5a*
a , n, K sba , n, K sca , n, K s
6a
a, n
b
a, n
c
a, n
7a*
a , n, K sb*a , n, K sc*a , n, K s
8a
a, n
b
a, n
c
a, n
9a
a , n, K sb*a , n, K sca , n, K s
10a
a, n
b
a, n
c
a, n
11a*
a , n, K sba , n, K sc*a , n, K s
12a
a, n
b
a, n
c
a, n
13a
a , n, K sba , n, K sca , n, K s
14a
a, n
b
a, n
c
a, n
15a
a , n, K sba , n, K sca , n, K s
16a
a, n
b
a, n
c
a, n
17a
a , n, K sba , n, K sca , n, K s
18a
a, n
b
a, n
c
a, n
Real parameters

Final Estimates

Measurement
Sets

a

n

Ks

a

n

Ks

F

v
v
v
v
v
v
Q
Q
Q
Q
Q
Q
h1
h1
h1
h1
h1
h1
u1
u1
u1
u1
u1
u1
Q, h 1
Q, h 1
Q, h 1
Q, h 1
Q, h 1
Q, h 1Q, u 1Q, u 1Q, u 1Q, u 1Q, u 1Q, u 1h 1, u 1
h 1, u 1
h 1, u 1
h 1, u 1
h 1, u 1
h 1, u 1
Q, h 1 , u 1
Q, h 1 , u 1
Q, h 1 , u 1
Q, h 1 , u 1
Q, h 1 , u 1
Q, h 1 , u 1
Q, h 1 , h 2 , u 1 , u 2
Q, h 1 , h 2 , u 1 , u 2
Q, h 1 , h 2 , u 1 , u 2
Q, h 1 , h 2 , u 1 , u 2
Q, h 1 , h 2 , u 1 , u 2
Q, h 1 , h 2 , u 1 , u 2
zzz

0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
0.02
0.03
0.006
zzz

1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
1.25
1.66
1.12
zzz

0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001zzz
zzz
zzz
0.00010.00010.0001
zzz
zzz
zzz
zzz

0.009220.009160.009020.015080.015470.015070.01038
0.02187
0.009350.014930.02398
0.01501
0.00761
0.015100.015020.015020.015020.015020.01057
0.009910.02003
0.01503
0.01503
0.01504
0.015070.01475
0.01497
0.015020.01503
0.015030.01038
0.01552
0.01930
0.015020.015020.015020.01503
0.01503
0.01504
0.01503
0.015020.015030.01498
0.01498
0.01497
0.015020.015020.015020.014710.01451
0.01501
0.015020.015020.015020.01502

1.715
1.716
1.722
1.594
1.605
1.594
1.677
3.195
1.640
1.589
1.931
1.592
1.010
1.580
1.592
1.592
1.592
1.592
1.712
1.739
1.519
1.592
1.592
1.592
1.592
3.500
1.596
1.592
1.592
1.593
1.719
1.583
1.528
1.592
1.592
1.592
1.592
1.593
1.593
1.592
1.592
1.592
1.593
1.593
1.593
1.592
1.592
1.592
1.598
1.604
1.593
1.592
1.592
1.592
1.592

0.000360.000360.00035zzz
zzz
zzz
0.000430.000370.00041zzz
zzz
zzz
0.00057
0.000710.00070zzz
zzz
zzz
0.000420.000380.00107zzz
zzz
zzz
0.000700.000250.00070zzz
zzz
zzz
0.000410.000730.00101zzz
zzz
zzz
0.000700.000700.00070zzz
zzz
zzz
0.000700.000700.00070zzz
zzz
zzz
0.000680.000670.00070zzz
zzz
zzz
0.000700.1718e202
0.1719e202
0.1853e202
0.9209e204
0.4022e203
0.9260e204
0.3482e105
0.3181e105
0.2113e1050.8023e1040.2283e105
0.5087e104
0.5309e105
0.4485e200
0.1326e202
0.1386e202
0.7290e203
0.9520e204
0.2598e204
0.2232e204
0.3984e204
0.2845e204
0.2850e204
0.2891e204
0.9381e101
0.5243e105
0.9209e101
0.8755e101
0.8751e101
0.9996e101
0.1285e203
0.4281e204
0.1100e203
0.3714e204
0.3714e204
0.3716e204
0.2870e204
0.2870e204
0.2875e204
0.2870e204
0.2870e204
0.2871e204
0.3744e204
0.3739e204
0.3734e204
0.3717e204
0.3825e204
0.3716e204
0.5297e204
0.7173e204
0.3682e204
0.3557e204
0.3565e204
0.3557e204
zzz*Unsuccessful runs.

ried out 18 different optimization runs. The various optimization cases differed in terms of the number of parameters being
optimized and the information used to define the objectivefunction. Either all three parameters a , n, and K s (cases 1,
3, z z z , 17) or only the shape parameters a and n (cases 2, 4,
6, z z z , 18) were fitted. Similarly, as for the response surfaces,
we used for the inverse procedure the objective functions defined by the three different sets of measurements, that is, thecumulative infiltration (cases 3 and 4), the pressure head (cases

5 and 6), the water content (cases 7 and 8), or certain combinations of these measurements (Q 1 h for cases 9 and 10;
Q 1 u for cases 11 and 12; h 1 u for cases 13 and 14; Q 1
h 1 u for cases 15 and 16 with one observation node; and Q 1
h 1 u for cases 17 and 18 with two observation nodes). In
addition, we used the infiltration rates for the inverse solution(cases 1 and 2). Each optimization case was repeated threetimes using different initial estimates of the shape parameters
a and n.

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Results of the different parameter estimation runs are summarized in Table 2. Unsuccessful runs, whose cases are marked
by an asterisk, are defined as those where the optimized parameters deviated by more than 5% from the true value. Notethat with only one exception, all inversion attempts involving
only the two parameters a and n were successful. The exception was a case where only cumulative infiltration data wereused to define the objective function. Using different initial
estimates for the same case actually did produce the correct
hydraulic parameter values.
The results were far less successful when all three parameters a , n, and K s were optimized simultaneously, yielding
unsuccessful runs in 13 out of 27 cases. The low success rate
was fully expected when the parameters were estimated fromeither the instantaneous or cumulative infiltration curve (thesetwo curves should roughly hold the same information). As wediscussed earlier, the same infiltration curve can be calculatedusing an infinite number of combinations of values for a , n,
and K s . However, we were somewhat surprised to have several
unsuccessful runs when parameters were optimized using theobjective functions F(u) or F(Q, u ) since the response surfaces for these two objective functions (Figures 10 and 12)
showed relatively well-defined minima. This result may be explained in two ways. First, the contours of F(u) in Figure 10b
as well as the valley around the minimum are almost parallel tothe K s axis, thus suggesting that the solution will not be very
sensitive to the values of the saturated hydraulic conductivity.
Second, the a -n, a -K s , and n-K s planes represent only three
two-dimensional cross sections of the full three-dimensionalparameter space. The behavior of the objective function inthese three planes can only suggest how the objective function
might behave in the rest of the three-dimensional space. Other
local minima of the objective function F may exist which do
not show up with the three cross-sectional planes. To verify this
explanation, we calculated several additional response surfaces
for the a -n parameter plane but now assuming different values
for K s (plots not further shown here). Several local minima of
the type shown in Figure 8a and also visible in Figures 11a and11b were found on these additional plots. We note that someof the very small local minima on these and other plots may
have been caused also, or exacerbated, by minor oscillations in
the numerical results. Still, our study indicates that having
well-defined minima in the three two-dimensional parameter
planes does not automatically guarantee that no other local
minima are present and hence that the hydraulic parameters
obtained by the inverse solution are always unique.

Conclusions
We have developed a numerical code for identifying soilhydraulic parameters from unsaturated flow data typically observed during a three-dimensional disc permeameter infiltration experiment. The results presented in this study indicatethat measurements of the instantaneous or cumulative infiltration rates alone will not provide a unique solution in thethree-dimensional parameter space a -n-K s . An analysis of
inverse solutions shows that similar conclusion also hold forthe water content when measured at only one or two locationsin the field. Although the response surfaces for the objectivefunction based on the water content showed well-defined minima in each of the three selected parameter planes ( a -n,
a -K s , and n-K s ), additional local minima in three-dimensional parameter space still caused nonuniqueness in the op-

2695

timized parameters. Additional information of the pressure
head measured in the profile significantly increased theuniqueness of the solution. Actually, using only pressure headmeasurements already guaranteed convergence to the correct
parameters. Having a combination of the objective functions
for the pressure head and the cumulative infiltration rate did
not improve the finale fit. Best identifiability was obtained
when both the pressure head and the water content were measured simultaneously. In this case the simultaneously measured
pressure head and water content data should directly define
several points of the retention curve, in which case the saturated hydraulic conductivity becomes the main unknown parameter driving the inverse solution process.
Identification of the two shape parameters a and n was
successful using any combination of measurements. This finding suggests that it may be possible to obtain from disc permeameter experiments the saturated hydraulic conductivity (inthe sense of the unsaturated conductivity as extrapolated tozero pressure head) based on the traditional analysis of Wooding [1968] and the assumption of having a Gardner-type exponential unsaturated hydraulic conductivity function and to
identify only the shape parameters a and n in the retention
curve by inverse solution of Richards’ equation.
Finally, we emphasize that this study was conducted using
numerically generated synthetic data. No additional noise wasadded to the generated data to reflect measurement or calibration errors, which are unavoidable in experimental studies.
Additional studies of how measurement errors will affect theresults may need to be carried out before the inverse technique
should be used with real data.
Acknowledgment. Thanks are due to one anonymous and one notso-anonymous (B. E. Clothier) reviewer for their many useful suggestions and comments.

References
Ankeny, M. D., M. Ahmed, T. C. Kaspar, and R. Horton, Simple field
method for determining unsaturated hydraulic conductivity, Soil Sci.
Soc. Am. J., 55, 467– 470, 1991.
Bard, Y., Nonlinear Parameter Estimation, 341 pp., Academic, San
Diego, Calif., 1974.
Beck, J. V., and K. J. Arnold, Parameter Estimation in Engineering andScience, John Wiley, New York, 1977.
Bohne, K., C. Roth, F. J. Leij, and M. T. van Genuchten, Rapid
method for estimating the unsaturated hydraulic conductivity frominfiltration measurement, Soil Sci., 155(4), 237–244, 1992.
Bruce, R. R., and A. Klute, The measurement of soil moisture diffusivity, Soil Sci. Soc. Am. Proc., 20, 458 – 462, 1956.
Celia, M. A., E. T. Bouloutas, and R. L. Zarba, A general massconservative numerical solution for the unsaturated flow equation,
Water Resour. Res., 26(7), 1483–1496, 1990.
Dane, J. H., and S. Hruska, In-situ determination of soil hydraulicproperties during drainage, Soil Sci. Soc. Am. J., 47, 619 – 624, 1983.
Eching, S. O., and J. W. Hopmans, Optimization of hydraulic functions
from transient outflow and soil water pressure data, Soil Sci. Soc.
Am. J., 57, 1167–1175, 1993.
Gardner, W. R., Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a watertable, Soil Sci., 85, 228 –232, 1958.
Green, R. E., L. R. Ahuja, and S. K. Chong, Hydraulic conductivity,
diffusivity, and sorptivity of unsaturated soils: Field methods, in
Methods of Soil Analysis, Agron. Monogr. 9, part 1, 2nd ed., edited by
A. Klute, pp. 771–798, Soil Sci. Soc. Am., Madison, Wis., 1986.
Hillel, D., and C. H. M. van Bavel, Simulation of profile water storageas related to soil hydrologic properties, Soil Sci. Soc. Am. J., 40,
807– 815, 1976.
Klute, A., and C. Dirksen, Hydraulic conductivity and diffusivity: Lab-

2696

ŠIMŮNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

oratory methods, in Methods of Soil Analysis, Agron. Monogr. 9, part1, 2nd ed., edited by A. Klute, pp. 687–734, Soil Sci. Soc. Am.,
Madison, Wis., 1986.
Knopman, D. A., and C. I. Voss, Behavior of sensitivities in the onedimensional advection-dispersion equation: Implication for parameter estimation and sampling design, Water Resour. Res., 23(2), 253–
272, 1987.
Kool, J. B., and J. C. Parker, Analysis of the inverse problem fortransient unsaturated flow, Water Resour. Res., 24(6), 817– 830, 1988.
Kool, J. B., J. C. Parker, and M. T. van Genuchten, ONESTEP: A
nonlinear parameter estimation program for evaluating soil hydraulic properties from one-step outflow experiments, Bull. Va. Agric.
Exp. Stn. 85-3, 1985a.
Kool, J. B., J. C. Parker, and M. T. van Genuchten, Determining soil
hydraulic properties from one-step outflow experiments by parameter estimation, I, Theory and numerical studies, Soil Sci. Soc.
Am. J., 49, 1348 –1354, 1985b.
Kool, J. B., J. C. Parker, and M. T. van Genuchten, Parameter estimation for unsaturated flow and transport models—A review, J.
Hydrol., 91, 255–293, 1987.
Logsdon, S. D., and D. B. Jaynes, Methodology for determining hydraulic conductivity with tension infiltrometers, Soil Sci. Soc. Am. J.,
57, 1426 –1431, 1993.
Logsdon, S. D., E. L. McCoy, R. R. Allmaras, and D. R. Linden,
Macropore characterization by indirect methods, Soil Sci., 155, 316 –
324, 1993.
Marquardt, D. W., An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math., 11, 431– 441, 1963.
Neuman, S. P., Calibration of distributed parameter groundwater flow
models viewed as a multiple-objective decision process under uncertainty, Water Resour. Res., 9(4), 1006 –1021, 1973.
Parker, J. C., and M. T. van Genuchten, Determining transport parameters from laboratory and field tracer experiments, Bull. Va.
Agric. Exp. Stn. 84-3, 96 pp., 1984.
Parker, J. C., J. B. Kool, and M. T. van Genuchten, Determining soil
hydraulic properties from one-step outflow experiments by parameter estimation, II, Experimental studies, Soil Sci. Soc. Am. J., 49,
1354 –1359, 1985.
Perroux, K. M., and I. White, Design for disc permeameters, Soil Sci.
Soc. Am. J., 52, 1205–1215, 1988.
Reynolds, W. D., and D. E. Elrick, Determination of hydraulic conductivity using a tension infiltrometer, Soil Sci. Soc. Am. J., 55,
633– 639, 1991.
Russo, D., E. Bresler, U. Shani, and J. C. Parker, Analysis of infiltration events in relation to determining soil hydraulic properties by
inverse problem methodology, Water Resour. Res., 27(6), 1361–1373,
1991.
Quadri, M. B., B. E. Clothier, R. Angulo-Jaramillo, M. Vauclin, and S.
R. Green, Axisymetric transport of water and solute underneath a
disk permeameter: Experiments and numerical model, Soil Sci. Soc.
Am. J., 58, 696 –703, 1994.
Šimůnek, J., M. Šejna, and M. T. van Genuchten, The HYDRUS-2D
software package for simulating water flow and solute transport in
two-dimensional variably saturated media, version 1.0, IGWMC—

TPS-53, Int. Ground Water Model. Cent., Colo. School of Mines,
Golden, Colo., 1996.
Smettem, K. R. J., and B. E. Clothier, Measuring unsaturated sorptivity and hydraulic conductivity using multiple disk permeameters, J.
Soil Sci., 40, 563–568, 1989.
Toorman, A. F., P. J. Wierenga, and R. G. Hills, Parameter estimationof hydraulic properties from one-step outflow data, Water Resour.
Res., 28(11), 3021–3028, 1992.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method for
determining soil hydraulic functions from one-step outflow experiment, Soil Sci. Soc. Am. Proc., 56, 1042–1050, 1992.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method todetermine soil hydraulic functions from multistep outflow experiment, Soil Sci. Soc. Am. Proc., 58, 647– 652, 1994.
van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
van Genuchten, M. T., Non-equilibrium transport parameters frommiscible displacement experiments, Res. Rep. 119, U.S. Salinity Lab.,
U.S. Dep. of Agric., Agric. Res. Serv., Riverside, Calif., 1981.
van Genuchten, M. T., and F. J. Leij, On estimating the hydraulicproperties of unsaturated soils, in Indirect Methods for Estimating theHydraulic Properties of Unsaturated Soils, edited by M. T. van Genuchten, F. J. Leij, and L. J. Lund, pp. 1–14, Univ. of Calif., Riverside,
1992.
van Genuchten, M. T., and J. Šimůnek, Evaluation of pollutant transport in the unsaturated zone, in Proceedings of Regional Approaches
to Water Pollution in the Environment, NATO ASI Ser., Ser 2, Environment, edited by P. E. Rijtema and V. Elias, Kluwer, Norwell,Mass., in press, 1996.
Vogeler, I., B. E. Clothier, S. R. Green, D. R. Scotter, and R. W.
Tillman, Characterizing water and solute movement by time domainreflectometry and disk permeametry, Soil Sci. Soc. Am. J., 60, 5–12,
1996.
Warrick, A. W., Model for disc infiltrometers, Water Resour. Res.,
28(5), 1319 –1327, 1992.
White, I., and M. J. Sully, Macroscopic and microscopic capillary
length and timescales from field infiltration, Water Resour. Res., 23,1514 –1522, 1987.
Wooding, R. A., Steady infiltration from large shallow circular pond,
Water Resour. Res., 4, 1259 –1273, 1968.
Yeh, W. W.-G., Review of parameter identification procedures ingroundwater hydrology: The inverse problem, Water Resour. Res.,
22(2), 95–108, 1986.
Zachmann, D. W., P. C. Duchateau, and A. Klute, Simultaneous approximation of water capacity and soil hydraulic conductivity byparameter identification, Soil Sci., 134, 157–163, 1981.
J. Šimůnek and M. T. van Genuchten, U.S. Salinity Laboratory, U.S.
Department of Agriculture, ARS, 450 W. Big Springs Rd., Riverside,
CA 92507-4617. (e-mail: jsimunek@ussc.ars.usda.gov)(Received February 15, 1996; revised May 13, 1996;
accepted May 17, 1996.)