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. Simfinek 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 KS and (Y 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, 0 (h), and unsaturated hydraulic conductivity, K(h), functions involving the water content 0,
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 ,b&zek, 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

SIMUNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

hydraulic conductivity at the wetting pressure head and therefore can be neglected. Since (2) contains two unknown parameters, K and 4, or K, and (Y 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, and (Y 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- 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
K(h) = K, exp (cub)(1) which combines Levenberg-Marquardt nonlinear parameter
optimization involving weighted least squares with a quasiAnalyses of this type require either two infiltration measurethree-dimensional
numerical model which solves the variably
ments using two different disc diameters [Smettem and Clothsaturated
flow
equation.
The unknown parameters in vanier, 1989] or measurements using a single disc diameter but
then imposing multiple tensions [Ankeny et al., 1991]. The Genuchten’s model for the unsaturated soil properties areparameters K, and (Y are generally estimated using Wooding’s estimated from observed cumulative infiltration data during[ 1968] equation for unconfined steady state infiltration from a transient water flow by numerical inversion of Richards’ equation. Additional data, such as measured pressure heads or
disc,water contents, as well as a penalty function for constraining
QW = ~r%hJ + 4r04hJ(2) the optimized parameters to remain in some feasible region
using Bayesian estimation, can be optionally included in thewhere Q is the steady state infiltration rate, rO is the radius of parameter estimation procedure. Objective functions are forthe disc, h, is the wetting pressure head, and phi is the matrix mulated in terms of cumulative infiltration, pressure head, andflux potential or, equivalently, the linearized Kirchhoff trans- water content data, and on a combination of these measureform given by [Gardner, 1958]
ments. The behavior of the objective functions in threedimensional parameter space is studied by means of a series ofhu
response surfaces in two-dimensional parameter planes. TheK(h)dh
4th h,) =(3)
parameter estimation procedure is demonstrated using numerIh,
ically generated data; the sensitivity of the method to different
in which hi is the initial pressure head. Equation (2) considers initial estimates of the model parameters will be discussed also.
the two main forces driving the infiltration process: the first
One important issue addressed in this paper is the question
term represents one-dimensional gravity-driven flow, while the of whether infiltration data usually collected during tension
second term accounts for capillary-induced flow [Wooding, disc permeameter experiments are sufficient for the parameter
1968]. Explicit in (2) is the assumption that the hydraulic con- estimation procedure or whether additional data are needed toductivity at the initial pressure head is much lower than the guarantee a unique inverse solution. We will also examine
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, and the sorptivity
parameter a! in Gardner’s exponential model [Gardner, 1958]of the unsaturated hydraulic conductivity,2686SIMUNEK 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:@(P) = [cl* - s(P)W’[cl* - s(P)1(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, pO( 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*(P) = CL(P)P”(P)

where Qi + , and Qj are the objective functions at the previous
and current iteration level, R, is a positive definite matrix, andpi 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 v, and/or the step size pi [Bard, 1974]. Steepest descentmethod selects pi = 1 and R, = I, where I is an identity matrix,
whereas Newton’s method uses pi = 1 and R, = Hz:‘, where
H is the Hessian matrix of Q(B):

(15)
x

II

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 @:Q(P) = ]4* - s(P)l’v-‘[cl* - s(P)1+ (p* - pyv,yp* - p,(16)

where B* is the parameter vector containing the prior information (e.g., 0,, 8,, CX, IZ, and K,), fi is the associated predicted parameter vector, and V, 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 V,, 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 V, 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
Q(P) = i w,[qT- c7l(P)l’
r=l

(17)

where w, 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 pi of the unknown parameters to be estimated andthen studying how the objective function Q,(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,

Hi,(P) = apap
1

(20)

I

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 = J’J(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:

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 + A), where h is apositive scalar. When A is large, then the matrix is diagonallydominant resulting in the steepest descend method. On theother hand, when h 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 AB against
rounding errors in differencing which decrease with AB. A
common practice is to change the optimized parameters by one
percent to obtain estimates of the sensitivity coefficients,

84,

qL(P + Ape,) - qL(P)

aP i

APj

-zz

(23)

where ej is the jth unit vector, and AB = 0.01 B.

Data Generation
Pi+, = I% + P~V, = Pi - /~iR,p,(18)

decreases the value of the objective function@,+1 < Q,(19)

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

SIMtJNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

2687

2
,

- 1 .

Node(lO.-10).... - -2. Node (20,.10)

0

20000

40000

60000

80000

20000

Time [s]

40000

60000

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: 0, = 0.000, 0, = 0.480, (Y = 0.0150 cm-‘,
n = 1.592, and KS = 0.0007 cm s-‘. 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 -500 cm, whereas the pressure head in the discwas set at -3 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-Figure 3. Calculated pressure head versus time at two observation nodes, (r, z) = ( 1 0 , - 10) and (20, -l0), 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 19 (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 0,
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:i(r, z, t) = b

O(r, z, t) dV”

1h(r, z, t) = ~V-1. Node (10,-10). ..---Z.Node(ZO.-10)

--t
20000

40000

80000

Time [s]

60000

80000

Time [s]Figure 2. Calculated water content versus time at two observation nodes, (r, z) = (10, -10) and (20, -l0), for the
hypothetical disc permeameter infiltration experiment.

(24)h(r, z, t) dVv

where i!l and & 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-

SIMUNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

2688
0.3

0.25

Experimental Design

T

iII’I

I:,
F

---n

’
i

0.2

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:-a______~-.-.K]-,,15 ,:.\I

Time [s]Figure 4. Sensitivity of the instantaneous infiltration rate,u(t), to a 1% change in hydraulic parameters (Y, II, tlq, and KS.

= lq;(P + Ape,) - qAP)I

ditional 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 Qi 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 [hrhaek et al., 1996]. The cumulative disc permeameter infiltration rate, Q(t), was calculated asfollows:(25)
where t, is the starting time of the experiment, u(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.

(26)

where e,, is the change in the auxiliary variable qi corresponding to a 1% change in parameter pi.
Figure 4 shows the sensitivity of the instantaneous infiltration rate, v(t), to 1% changes in the parameters (Y, n, 19,, and
KS. Since the residual water content was assumed to be zero,
we did not use 0, in our sensitivity study. Because of physicalconsiderations, we expected that there would be a certain timeinterval (0, tgrav) at the beginning of the infiltration event
during which the infiltration rate would be very sensitive to theunknown parameters. The characteristic time tgrav 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, u 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

0

10000

20000

30000

40000

50000

60000

70000

80000

90000

Time [s]Figure 5. Sensitivity of the cumulative infiltration rate, Q(t),to a 1% change in hydraulic parameters (Y, n, O,, and KS.

2689
parameters. The cumulative infiltration rate is most sensitive tothe parameter 12 and least to the saturated water content 0,.
The sensitivity to IZ is about three times higher than the sensitivities to a! and the saturated hydraulic conductivity KS,
which are about the same. If the absolute values in (26) wereto ‘be removed, the sensitivity to CY would get a sign opposite tothat of KS. 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 KS, 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 (Y, IZ, 0,, and KS at observation nodes(r, z) located at (10, -10) (Figure 6a) and (20, -10) (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

0

1000

2000

3000

4000

5000

6000

7000

8000

Time [s]
0.005

Tb

0

10000

5000

15000

Time [s]Figure 7. Sensitivity of the water content at two observation
nodes (a) (r, z) = (10, -10) and (b) (r, z) = (20, - 10) to
a 1% change in hydraulic parameters 01, n, 0,, and K,.

0

1000

2000

3000

4000
Time [s]

i
5000

6000

7000

-a
___”

______#-.-.j(

Figure 6. Sensitivity of the pressure head at two observation
nodes (a) (r, z) = (10, - 10) and (b) (Y, z) = (20, -10) to
a 1% change in hydraulic parameters (Y, n, tl,, and K,y.

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 01, n, O,, and K, at thesame two observation nodes, (r, z) = (10, -10) and (20,-10), are shown in Figure 7. The initial sensitivities of 0 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 (CC,n, 0,) but not to the saturated hydraulic conductivity, K,v. The
water content during passage of the moisture front is very
sensitive to 01 and KS but not very sensitive to n and 8,. After
passage of the moisture front the sensitivity to o( and KS decreases gradually, although more slowly than for the pressure
head (Figure 6). By contrast, the sensitivity to the saturated
water content 0, 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 0,. Since thewater content remains sensitive to at least one hydraulic parameter during the entire infiltration event, water content mea-

SIMUNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

2690

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

Upper

Parameter

Parameter
Value

Parameter
Step Value

ParameterValuea, cm-’n [-]KS, cm SC’

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@(p, qm) = i
j=li

Wji Wt[q,Yti) i=l

4jCtt,

PII’

1

(27)

where m represents the different sets of measurements, that is,
those involving cumulative infiltration, pressure head, or watercontent data; nj is the number of measurements in a particular
set; qT(ti) are specific measurements at time ti for the jthmeasurement set; qj(ti, p) are the corresponding model predictions for the parameter vector 6; and wj and wi are weights
associated with a particular measurement set or point, respectively. We assumed that the weighting coefficients wi 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,, was also set
equal to 1, while those for the cumulative infiltration, w,, andthe pressure head, wh, sets were defined as

W

e=

1
1

-1

(28)
-1

(29)

respectively, where IZ 8, n/2, and ng represent the total numberof measurements obtained for the water content, pressure
head, and cumulative infiltration, respectively; and q *,, q i, andqk are the actual measurements of 0, h, and Q, respectively.
We calculated the objective functions for three parameter
planes (i.e., a-n, a-K,, and n-K,), 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 + h + 0, Q + h, Q + 0, and h + 0). 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 (a(Q) for the cumulative infiltration measurements (seealso Figure 1) for the three different parameter planes, a-n,a-K,, and n-K,. The a-n response surface (Figure 8a) showsa well-defined valley which starts at low 01 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 a(Q). We obtained
visually almost identical response surface in the parameter
plane ol-KS (Fig. 8b), although now several local minima appeared for relatively small values of the parameters (Y and KS.
Figures 8a and 8b indicate that a higher value of (Y and corresponding increases in n or K,? 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 KS interms of their effect on the objective function. An increase inn and a corresponding decrease in KS will produce the same
response in the objective function. This hyperbolic behavior of@(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 KS is

SIMdNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

a)

2691

b)

Alpha [l/cm]

Alpha [l/cm]

n [-lFigure 8. Contours of the objective function Q(Q) for the cumulative infiltration rate in the (a) 01-n plane,(b) a-K, plane, and (c) n-K, 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 II 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 (Y,n, and KS. His conclusion holds also in three dimensional
space for tension disc permeameters.
Figure 9 represents response surfaces of@(h) for the pressure head measured at ( I, z) = (20, -10) (Figure 3), again
in the three parameter planes o-n, a-K,, and n-K,. The
response surfaces for the h-based objective function, Q(h),are visually quite similar to those for the cumulative infiltration
rate, Q(Q). An increase in (Y and a simultaneous increases inK, or n again lead to approximately the same values for theobjective function. The contours of O(h) and the valleyaround its minimum in Figure 9a are almost parallel with thea)n axis, thus indicating that the parameter n is very difficult to
identify from available information. And similarly, as @(Q) in
Figure 8c, the pressure head based objective function @(h)reacts inversely to an increase in the parameters IZ and KS(Figure 9c). The contours of a(h) at relatively low values ofKS, as well as the long valley around the minimum, are almost
parallel with the II axis (Figure 9c), thus showing that 12 can not
be identified precisely when all three parameters 01, n, and KS
are to be estimated simultaneously. Still, the minimum of a(h)is somewhat better identifiable than the minimum of Q(Q).
Figure 10 shows similar response surfaces for the objectivefunction Q(0) for the water content measured at (r, z) =(20, - 10) (Figure 2). The objective function O(0) 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 0.

b)

Alpha [1/cm]c)Alpha [l/cm]

n [-]Figure 9. Contours of the objective function @(h) for the pressure head measured at (r, z) = (20, - 10)
in the (a) cr-n plane, (b) CY-KS plane, and (c) n-K, parameter plane.

2692

SIMdNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIESb)

Alpha [i/cm]

Alpha [l/cm]

Figure 10. Contours of objective function O( 0) for the water content measured at (r, z) = (20, - 10) inthe (a) (~-n plane, (b) a-K, plane, and (c) n-K, 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 @(Q, h)(= Q(Q) + a(h)) defined by both the cumulative infiltrationand the pressure head measured at (r, z) = (20, - 10). 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 Cp (Q) and m(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 (Y-n and a-K, (Figures lla and llb).
Better results are obtained when the objective functions
based on the cumulative infiltration, Q(Q), and the watercontent measured at (r, z) = (20, -10), a(0), are combined. Figure 12 shows that the response surfaces for thisobjective function @(Q, 19) are now well defined in all threeplanes (a-n, a-KS, and n-K,). An exception is perhaps thea-K, plane, which shows a valley that extends over a relativewide range of (Y and KS values. Finally, and for completeness,
we also present in Figures 13 and 14 contours of objectivefunctions @(h , 0 ) and @(Q, h , 19 ) , 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-a)

Alpha [l/cm]

Alpha [1/cm]

n I-l

Contours of the objective function @(Q, h) for the cumulative infiltration rate and the pressure
head measured at (r, z) = (20, -10) in the (a) (~-n plane, (b) a-K, plane, and (c) n-K, parameter plane.
Figure 11.

2694SIMUNEK AND VAN GENUCHTEN: ESTIMATING UNSATURATED SOIL HYDRAULIC PROPERTIES

Table 2.

Results Obtained With the Different Parameter Estimation Computer Runs

Case

FittedParameter

Initial Estimates

la*
b*
C*

2a
b
3:*
b*
C*

4a
b*
sea*
b
6:
b
7:*
b*
C*

8a
b
9:
b*
C

10a
b
11:*
b
C*

12a
b
C

13a
b
14:
b
C

15a
b
16:
b
17:
b
C

18a
b
a, n
Real parameters
C

Measurement
Setsv

Final Estimates

(Y

n

K,

(Y

n

K,

@

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
...

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
...

0.00010.00010.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.00035...

0.1718e-02
0.1719e-02
0.1853e-02
0.9209e-04
0.4022e-03
0.9260e-04
0.3482e+05
0.3181e+05
0.2113e+05
0.8023e+04
0.2283e+05
0.5087e+04
0.5309e+05
0.4485e-00
0.1326e-02
0.1386e-02
0.7290e-03
0.9520e-04
0.2598e-04
0.2232e-04
0.3984e-04
0.2845e-04
0.2850e-04
0.2891e-04
0.9381e+01
0.5243e+05
0.9209e+01
0.8755e+01
0.8751e+01
0.9996e+0l
0.1285e-03
0.4281e-04
0.1100e-03
0.3714e-04
0.3714e-04
0.3716e-04
0.2870e-04
0.2870e-04
0.2875e-04
0.2870e-04
0.2870e-04
0.2871e-04
0.3744e-04
0.3739e-04
0.3734e-04
0.3717e-04
0.3825e-04
0.3716e-04
0.5297e -04
0.7173e-04
0.3682e-04
0.3557e-04
0.3565e-04
0.3557e-04

0.0001...
...
...
0.00010.00010.0001...
...
...
0.00010.00010.0001...
...

...
0.00010.00010.0001...
...

...
0.00010.00010.0001...
...
...
0.00010.00010.0001

...
...
...

0.00010.00010.0001...
...
...
0.00010.00010.0001
...
...

...

0.00010.00010.0001

...
...
...
...

...
...

0.000430.000370.00041...
...
...
0.00057
0.000710.00070

...
...
...

0.000420.000380.00107

...
...
...

0.000700.000250.00070...
...
...
0.000410.000730.00101

...
...
...

0.000700.000700.00070

...
...
...

0.000700.000700.00070
...

...
...
0.000680.000670.00070

...
...
...

0.00070...

*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 01, n, and KS (cases 1,
3;**, 17) or only the shape parameters cx and 12 (cases 2, 4,
6;.., 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 + h for cases 9 and 1 0 ;Q + 0 for cases 11 and 12; h + 0 for cases 13 and 14; Q +h + 0 for cases 15 and 16 with one observation node; and Q +h + 0 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
01 and n.

SIM~NEK A ND VAN GENUCHTEN:E S T I M A T IN G U N S A T U R A T E D SOIL

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 (Y, n, and KS 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 (Y, n,
and KS. However, we were somewhat surprised to have several
unsuccessful runs when parameters were optimized using theobjective functions @(19) or @(Q, 0) 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 a(0) in Figure 10b
as well as the valley around the minimum are almost parallel tothe KS axis, thus suggesting that the solution will not be very
sensitive to the values of the saturated hydraulic conductivity.
Second, the (y-n, a-K,, and n-K, 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 @ 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 (~-n parameter plane but now assuming different values
for KS (plots not further shown here). Several local minima of
the type shown in Figure 8a and also visible in Figures lla andllb 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 (Y-~-K,~. 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,, and n-K,), additional local minima in three-dimensional parameter space still caused nonuniqueness in the op-

HYDRAULICHYDRUALICHYDRAULICPROPERTIES
C

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 01 and IZ wassuccessful 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 12 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 diffusivitv. Soil Sci. Soc. Am. Proc., 20. 458-462. 1956.
Celia, M. A., E. T. Bouloutas, anh 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, 19.58.
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

SIM~NEK A ND VAN GENUCHTEN: ES T IM A T I N G

oratory methods, in Methods of Soil Analysis, Agron. Monogr. 9, part1. 2nd ed.. edited bv A. Klute. p p . 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 Resow. Res., 23(2), 253272, 1987.
Kool, J . B., and J. C. Parker, Analysis of the inverse problem fortransient unsaturated flow, Water Resow. 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., l55,316324, 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.
SimBnek, J., M. Sejna, 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-U N S AT U R A T E D

soIL H Y DR AU L I C P RO P E R T I E S

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 Resow.
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. l-14, Univ. of Calif., Riverside,
1992.
van Genuchten, M. T., and J. Simunek, Evaluation of pollutant transport in the unsaturated zone, in Proceedings of Regional Approaches
to Water Pollution in the Environment, NATO ASZ 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 Resow. Res.,
28(5), 1319-1327, 1992.
White, I., and M. J. Sully, Macroscopic and microscopic capillary
length and timescales from field infiltration, Water Resow. Res., 23,
1514-1522, 1987.
Wooding, R. A., Steady infiltration from large shallow circular pond,
Water Resow. 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. Simunek 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.)