WATER RESOURCES RESEARCH, VOL. 33, NO. 8, PAGES 1879 –1892, AUGUST 1997

A sequential uncertainty domain inverse procedure for
estimating subsurface flow and transport parameters
K. C. Abbaspour,1 M. T. van Genuchten,2 R. Schulin,1 and E. Schläppi3
Abstract. A parameter estimation procedure, sequential uncertainty domain parameter
fitting (SUFI), is presented and has the following characteristics. The procedure issequential in nature, meaning that one more iteration can always be made before choosing
the final estimates. The procedure has a Bayesian framework, indicating that the method
operates within uncertainty domains (prior, posterior) associated with each parameter.
The procedure is a fitting procedure, conditioning the unknown parameter estimates on anarray of observed values. Finally, the procedure is iterative, requiring a stopping rule
which is provided by a critical value of a goal function. Performance of the SUFI
parameter estimation procedure is demonstrated using three examples of increasing
complexity: (1) analysis of a solute breakthrough curve measured in the laboratory during
steady state water flow, (2) estimation of the unsaturated soil hydraulic parameters from atransient drainage experiment carried out in a 6-m deep lysimeter, and (3) estimation of
selected flow and transport parameters from a hypothetical ring infiltrometer experiment.
The procedure was found to be general, stable, and always convergent.

Introduction
During the past several decades, considerable advances havebeen made in the conceptual understanding and mathematicaldescription of water flow and solute transport processes invariably saturated soils and groundwater systems. A large number of models of different degrees of complexity and dimensionality are now available for predicting subsurface flow andtransport [e.g., National Research Council (NRC), 1990; Ségol,1994]. Still, effective application of these models to practicalfield situations suffers from two general problems. One issue
concerns the inability of models, however complex and sophisticated numerically, to consider all of the important simultaneous physicochemical and hydrologic processes operative inthe subsurface, such as the presence of physical and chemical
nonequilibrium transport processes, preferential flow, various
boundary processes, and spatial and temporal variability in soil
hydraulic properties. The second issue is that even if a model
were to be developed that would account for all pertinent
processes, model users would still have difficulty collecting
enough meaningful field data to effectively run the model. The
past 30 years or so have probably seen more advances innumerical modeling per se than in the methodology of fielddata collection and how best to deal with parameter uncertainty in predictive modeling.
Various schemes have been proposed to somehow deal withthe problem of uncertainty in model input parameters. The
problems of parameter estimation and uncertainty should be
addressed in a comprehensive manner using all available information. No experimental work is ever performed in a com1

Department of Soil Protection, Swiss Federal Institute of Technology, Schlieren, Switzerland.
2
U.S. Salinity Laboratory, Agricultural Research Service, U.S. Department of Agriculture, Riverside, California.
3
Colombi Schmutz Dorthe, Bern, Switzerland.
Copyright 1997 by the American Geophysical Union.
Paper number 97WR01230.
0043-1397/97/97WR-01230$09.00

plete vacuum, and prior information is nearly always available
in some manner; such information must be used to advantage
by providing direction to further sampling. Nearly every project
involves the collection of data of different types and qualities;
these data should all be used to improve the understanding of the
problem. For example, most environmental problems are spatialin nature, characterized by autocorrelated subsurface flow andtransport parameters; autocorrelation in parameters hence must
be revealed and taken advantage of. Moreover, measured valuesof the primary outputs of hydrologic simulation models often
reveal the nature of the input fields and hence should also be
utilized. While the collection of new or additional data is generally the most direct way of reducing input uncertainty, limiting
factors such as cost, time, site destruction, and accuracy may notalways permit further sampling. Comprehensive algorithms areneeded to link the above seemingly different issues and considerations. An algorithm developed by Abbaspour et al. [1996], theBayesian uncertainty development algorithm (BUDA), can utilize all of the above techniques to achieve a higher reduction in
uncertainty in environmental projects. The inverse procedure inthis work is part of BUDA.
Of the different options mentioned above the use of measured primary outputs of hydrologic simulation models to estimate model input parameters has received special attention.
Much literature is devoted to this process alternatively referred
to as inverse modeling, model calibration, parameter fitting,parameter estimation, and history matching. Reviews of thesubject are given by Yeh [1986] and Kool et al. [1987]. The
procedure generally involves minimization of the square difference function of some measured and simulated flow ortransport variable. Yeh and Yoon [1981] studied the parameter
identifiability of an aquifer with optimum dimension in parameterization. Cooley [1982, 1983] incorporated prior information
about the parameters into a nonlinear regression groundwaterflow model. Parker and van Genuchten [1984] applied an inverse technique to laboratory and field tracer experiments,
while Parker et al. [1985] used one-step laboratory outflow
experiments to estimate the unsaturated soil hydraulic param-

1879

1880

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Figure 1. Division of a distributed parameter into equally
sized strata for exhaustive stratified sampling. The dashed lines
locate the first moments of the strata.

eters. Carrera and Neuman [1986] employed a maximum likelihood method and incorporated prior information in theiraquifer parameter estimation procedure assuming transientand steady state conditions. Van Dam et al. [1990] modified theone-step outflow technique of Parker et al. [1985] into a multistep procedure for determining soil hydraulic functions. Eching and Hopmans [1993], similarly, studied optimization of thehydraulic functions using both transient outflow and soil waterpressure head data. Dane and Hruska [1983], Kool et al. [1987],
and Sisson and van Genuchten [1991] determined soil hydraulicproperties from in situ water content data during verticaldrainage, whereas Kool and Parker [1988] applied an inverse
procedure to a flow process consisting of ponded infiltration
followed by gravity drainage with evaporation at the soil surface. Some studies combined parameter estimation with hydraulic scaling capabilities [Shouse et al., 1992; Eching et al.,1994], while others coupled numerical inverse problems suchas heat or mass transport with unsaturated or saturated flow
[Carrera, 1987; Mishra and Parker, 1989; Sun and Yeh, 1990].
Yet others introduced geostatistics considerations in the parameter estimation process [Clifton and Neuman, 1982; Kitanidis and Vomvoris, 1983; Dagan, 1985; Hoeksema and Kitanidis,
1984; Kitanidis, 1995; Yeh et al., 1996].
A common problem with most of the inverse methods isstability and convergence [Yeh, 1986], and more robust procedures are desirable. The objective of this work is to describe
and demonstrate the use of a very different approach, sequential uncertainty domain parameter fitting (SUFI), for parameter estimation. The procedure is general, forward, sequential,
iterative, and Bayesian in nature. The method begins with prioruncertainty domains on the input parameters, usually invoking
relatively large uncertainties, and, subsequently, conditions themodel parameters on the measured data through an objective
or goal function. Ultimately, posterior uncertainty domains are
obtained with much reduced uncertainties. The stopping (convergence) rule is dictated by a tolerance imposed on the goal
function. On the basis of our experience in working with SUFIwith different problems, three of which are discussed later inthis work, the proposed procedure was found to be stable,
always converging, and also well suited for global optimization.

Sequential Uncertainty Domain
Parameter Fitting
The first step in running SUFI is to identify the domain of
uncertainty for each parameter. This is achieved by invoking aprobabilistic description of the input data. Unfortunately, expressing data in probabilistic forms is not a common practice;still, given the experimental difficulties in determining most
inputs, all parameters are subject to some uncertainty andhence projects should benefit from probabilistic rather thanabsolute statements regarding the input data. The form of theprobability function for a given parameter will depend on theinformation available for that parameter. For example, in depicting hydraulic conductivity, knowledge of a lognormal distribution will speed up the process of convergence. In a worst-

Figure 2. A schematic representation of uncertainty propagation in an exhaustive stratified sampling. The
250 possible values of the goal function are used to build a cumulative frequency distribution of the goal.

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

1881

Figure 3. Frequency distribution of the number of hits in each stratum. The strata with the largest numberof hits are most likely to contain the desired estimates.

case scenario the uncertain domain of a parameter can bedepicted as a finite interval, i.e., a uniform probability distribution between physically realistic lower and upper bounds ofthat parameter. In addition, nominal parameters can be considered in the analysis. For example, if the type of boundary
condition to be imposed is uncertain, a variable A can consistof A 5 {constant head, constant flux}. Similarly, a variable H
could depict uncertainty in the type of chemical equilibrium or
nonequilibrium model to be adopted, i.e., H 5 {local equilibrium, first-order kinetics, two-site sorption kinetics}.
During the second step, uncertainty in each parameter must
be propagated through a simulation model to a goal function.
Some of the commonly used goal functions in parameter estimation procedures are of the form

not yet been given due attention in the literature. This issue isfurther discussed at the end of this paper.
If measured values ( x im , i 5 1, z z z , N) contain error, thenwe treat each measured point as a random variable and obtain
the Bayesian distribution [Benjamin and Cornell, 1970] of thegoal function as
fbG ~ g! 5

E E
···

1

fG ~ gux1m , · · · , xNm ! f~x1m ! · · · f~xNm ! dx1m · · · dxNm

N

(4)

Absolute error

O
N

g5

ux m 2 x su i

(1)

i51

Root-mean-square error (RMSE)
g5

Î

1
N

O
N

~ x m 2 x s! 2i

(2)

i51

Logarithmic form of RMSE (LRMSE)
g5

Î

1
N

O
N

~log x m 2 log x s! 2i

(3)

i51

where x m is a measured value, x s is a simulated value, and N is
the number of measurements. We emphasize here that different goal functions can lead to quite different values for theestimated parameters. While a suitable goal function must
always be formulated, its form may depend also on the objectives of the problem being investigated. The potentially very
important effects of different goal functions on the results have

Figure 4. Comparison of measured (circles) effluent andsimulated (solid line) concentrations using parameters obtained by both van Genuchten and Wierenga [1986] and thesequential uncertainty domain parameter fitting (SUFI) estimation approach.

1882

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Table 1. Parameters of Chromium Column Displacement
Experiment in Example 1

Parameters

Values Obtained by
van Genuchten and
Wierenga [1986]

Values Obtained
by SUFI*

19.65
1.349
0.01406

19.66†
1.348‡
0.01398

Peclet number P
Retardation factor R
Value of goal function
expressed as (2)*SUFI is defined as the sequential uncertainty domain parameter
fitting.
†The range of uncertainty is (19.659 –19.665).
‡The range of uncertainty is (1.34766 –1.34772).

SUFI employs a Latin hypercube procedure [McKay, 1979] toobtain a solution of (4).
To propagate uncertainty in the parameters, SUFI can invoke either an exhaustive stratified sampling scheme or a random stratified sampling scheme. The stratified sampling procedure is based on a division of the cumulative probability
scale range into equal probability strata for each uncertain
input parameter. The first moment of each stratum on theparameter scale range is taken to represent that stratum. The
procedure is illustrated in Figure 1 for a lognormally distributed parameter such as the saturated hydraulic conductivity.
For parameters having uniform probability distributions theinterval depicting the range of the parameter is simply divided
into equal strata and the midpoint of each stratum is chosen to
represent that stratum. Categorical or nominal parameters are
already discretized and need no further treatment. In the exhaustive stratified sampling a hydrologic simulation program issubsequently run for all possible combinations of strata for all
input parameters, and for each run the goal function is calculated (Figure 2). For the case of the random stratified sampling, only a random subset of the exhaustive case is used for
simulation. This exercise allows uncertainties in the input parameters to be propagated to the goal function. The choice ofthe exhaustive stratified sampling or the random stratified sampling depends on the number of parameters and the speed of
the simulation program. If there are too many parameters or
the simulation program is too slow, then the random stratified
sampling should be invoked.
Each parameter stratum in the SUFI analysis carries a score
which is initially set to zero. After each simulation the goal
function is calculated, and if it meets a certain tolerance cri-

terion (e.g., when g , 0.2), a hit is obtained and a positive
point is added to the score of each of the strata of that run. For
the case where the goal function possesses a distribution as theresult of measurement errors, if the tolerance is greater thanthe lower limit on the 95% confidence interval of the goal, thenwe obtain a hit. When all runs are completed, a frequency
distribution of hits is constructed for each stratum as illustrated
in Figure 3, where the numbers on the horizontal axes indicate
the range of the parameters and the numbers on each column
refer to the number of hits. Next, strata having zero or a small
number of hits at both ends of each interval are eliminated,thus providing an updated uncertainty domain for the next
iteration. If the initial estimate of the uncertainty range for aparameter was too small and the actual value was to the right of
the interval, then all hits would fall on the rightmost stratum asshown in Figure 3 for the porosity. For the next iteration therefore we would select the interval to be [0.28, 0.32], for example.
With each iteration the uncertainty domain for each parametershould get smaller, while the goal function will also decrease.
Choosing the number of strata for each parameter should be
problem dependent. While having more strata may promote
faster convergence because of possibly fewer iterations, too
many strata could result in an unacceptable number of simulations for each iteration, especially when the simulation modelruns slowly and many parameters are to be estimated. For
example, if 10 parameters are present and each parameter
range is divided into 3 strata, then the total number of simulations required would be 59,049 for the exhaustive stratified
sampling case. For this scenario, as mentioned above, SUFI
allows random stratified sampling which selects randomly asubset of simulations to be performed. For example, if we limit
the number of simulations to 1% of the above exhaustive case,then each of the 59,049 simulations would be performed witha probability of 1%. This would reduce the number if simulations to about 590. Another option in SUFI for obtaining afaster result is to run the program in parallel. For example, if
we would like to perform 4000 runs in an iteration, we could
submit simultaneously 10 jobs of 400 runs each.

Test Cases and Parameter Estimation Results
The performance of the sequential uncertainty domain parameter fitting (SUFI) procedure is illustrated below by means of
three examples of increasing complexity: (1) analysis of a solute
breakthrough curve measured in the laboratory during steadystate water flow, (2) estimation of unsaturated soil hydraulic pa-

Table 2. Progression to Convergence for Example 1
Numberof
Iterations
1
2
3
4
5
6

Parameter Value
Parameters

Minimum

Maximum

Number ofStrataNumber of Hits in Each Stratum

Tolerance

P
R
P
R
P
R
P
R
P
R
P
R

5.0
1.0
9.5
1.2
15.8
1.32
18.32
1.338
19.2651.345219.492
1.347

50.0
2.0
41.0
1.5
22.1
1.38
20.21
1.356
20.0211.351519.73
1.3492

10
10
10
10
10
10
20
20
20
20
40
40

0, 1, 2, 3, 2, 1, 1, 1, 0, 0
0, 0, 2, 7, 2, 0, 0, 0, 0, 0
0, 0, 2, 2, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 2, 2, 0, 0, 0, 0
0, 0, 0, 0, 0, 2, 2, 0, 0, 0
0, 0, 0, 0, 2, 2, 0, 0, 0, 0
10 3 0, 4, 6, 7, 7, 6, 6, 5, 3, 0, 0
8 3 0, 5, 7, 8, 8, 7, 6, 3, 0, 0, 0, 0, 0
7 3 0, 1, 2, 1, 10 3 0
9 3 0, 2, 0, 2, 8 3 0
28 3 0, 1, 11 3 0
12 3 0, 1, 27 3 0

0.06
0.02
0.01440.01410.01400.01398

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Table 3.
Numberof
Iterations

1883

Distribution of the Number of Hits for Iteration 1 in Table 2 as a Function of Tolerance
Parameter Value
Parameters

Minimum

Maximum

Numberof Strata

P
R
P
R
P
R
P
R
P
R
P
R

5.0
1.0
5.0
1.0
5.0
1.0
5.0
1.0
5.0
1.0
5.0
1.0

50.0
2.0
50.0
2.0
50.0
2.0
50.0
2.0
50.0
2.0
50.0
2.0

10
10
10
10
10
10
10
10
10
10
10
10

1

rameters from a transient drainage experiment carried out usinga large 6-m deep lysimeter, and (3) estimation of selected flow
and transport parameters from a hypothetical ring infiltrometer
experiment. We also used example 1 to show the effect of a 50%
error in measured data, and we performed example 3 with both
exhaustive and random stratified sampling schemes.

c e~T! 5

Number of Hits in Each Stratum
2, 3, 3, 3, 3, 3, 3, 3, 3, 20, 0, 9, 10, 9, 0, 0, 0, 0, 0
0, 3, 3, 3, 3, 3, 2, 2, 2, 1
0, 0, 8, 9, 5, 0, 0, 0, 0, 0
0, 1, 2, 3, 2, 1, 1, 1, 0, 0
0, 0, 2, 7, 2, 0, 0, 0, 0, 0
0, 0, 1, 1, 1, 1, 0, 0, 0, 0
0, 0, 0, 4, 0, 0, 0, 0, 0, 0
0, 0, 0, 1, 0, 0, 0, 0, 0, 0
0, 0, 0, 1, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 0, 0, 0, 0

1
erfc
2
1

Example 1: Analysis of Observed Solute BreakthroughCurve
The first example deals with a column experiment treatedpreviously by van Genuchten and Wierenga [1986]. In this experiment, transport of chromium was studied through a 5-cm
long column of sand. Figure 4 shows the observed effluent
curve. The governing transport model for this example is the
one-dimensional convection dispersion equation
R

­C
­C
­ 2C
2v
5D
­t
­ x2
­x

(5)

subject to the initial and boundary conditions
C~ x, 0! 5 0

S

2D

(6)

D

­C
1 vC u x50 5 vC 0
­x

F ÎÎ

P

2 RT

1 P
e erfc2~R 2 T!

F ÎÎ

0.08
0.06
0.04
0.02
,0.02

G

(9)
(10)

P 5 vL/D(11)

where T is dimensionless time (pore volume) and P is thecolumn Peclet number. Measured c e (T) data are plotted in
Figure 4.
Several methods for estimating the unknown parameters P
and R in (9), i.e., trial and error, analysis of the slope of theeffluent curve, analysis of a lognormal plot of the data, andusing a nonlinear least squares analysis based on Marquardt’s[1963] maximum neighborhood method, are discussed by van
Genuchten and Wierenga [1986]. Of these four methods theleast squares approach was shown to be the most accurate andobjective method. The least squares approach minimized thefollowing goal function

O
N

g5
(8)

where C is the volume-averaged concentration, R is the retardation factor accounting for linear sorption, D is the dispersion
coefficient, x is the distance, t is the time, v is the pore watervelocity, and C 0 is the inlet concentration. The analytical solution of this transport problem for the flux-averaged relative
effluent concentration c e at the end of the column x 5 L isgiven by~R 1 T!

0.10

T 5 vt/L(7)~­C/­ x!~`, t! 5 0

P

2 RT

G

Tolerance~c e 2 c! 2

(12)

i51

where c e is the calculated concentration from (9), c is theobserved effluent concentration, and N is the number of observed data points. The residual sum of squares, g, was minimized using the least squares program of van Genuchten[1980], leading to the values for P and R listed in Table 1. The
results obtained with SUFI are also listed in Table 1 and shownin Figure 4 as the solid curve which is the same as the curve

Table 4. Progression to Convergence for Example 1 for the Case With Measurement Errors
Parameter Value
Number of
Iterations
1
2
3

Parameters

Minimum

Maximum

Number ofStrata

P
R
P
R
P
R

5.0
1.0
5.0
1.3
5.0
1.3

50.0
2.0
27.5
1.7
20.75
1.58

10
10
10
10
10
10

Number of Hits in Each Stratum
3, 2, 3, 1, 1, 0, 0, 0, 0, 0
0, 0, 0, 0, 1, 5, 3, 1, 0, 0
1, 3, 3, 3, 3, 3, 0, 0, 0, 0
1, 0, 3, 4, 4, 3, 1, 0, 0, 0
1, 1, 3, 5, 5, 6, 3, 3, 2, 1
1, 0, 0, 3, 5, 5, 4, 5, 4, 3

Tolerance
0.215
0.208
0.208

1884

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Table 5. Soil Hydraulic Parameters of Crushed Bandelier Tuff Used in the Caisson Drainage Experiments (Example 2)
Values Obtained by van Genuchten et al.
[1987]

Values Obtained by SUFI

HydraulicParameters

Method 1

Method 2

Scenario 1

Scenario 2

Scenario 3

urusa, cm21n
K s , cm d21
lValue of goal
function

0.0*0.331*
0.01433
1.506
25.0
0.5*
0.011623(equation (16))

0.02550.33200.015451.474
33.71
0.49460.015655(equation (2))

0.0*0.331*
0.005025
1.79
9.85
0.5*
0.0096(equation (2))

0.0*0.331*
0.011915
1.579
12.45
0.5*
0.004945(equation (16))

0.043083
0.36975
0.00595
1.867528.1670.5*
0.007238(equation (2))

Parameters were estimated using method 1, estimating three parameters using water content u and pressure heads h; method 2, estimating
six parameters using u and h; scenario 1, estimating three parameters using u only; scenario 2, estimating three parameters using u and h; andscenario 3, estimating five parameters using u only.
*Parameter values assumed to be known.

obtained by van Genuchten and Wierenga [1986]. The value ofthe goal function based on the RMSE formulation (equation(2)) using estimated parameters obtained by SUFI was essentially the same as that obtained using the least squares approach of van Genuchten and Wierenga [1986].
Table 2 shows the progression of SUFI to convergence using
as initial parameter estimates uniform distributions in the intervals P 5 [5.0, 50.0] and R 5 [1.0, 2.0]. The minimum
and maximum parameter values in Table 2 depict the lower
and upper limits of the parameter uncertainty domains; noticethat these intervals become smaller and smaller as the iterative
search proceeds. As mentioned earlier, no exact rules exist forchoosing the number of strata. A large number of strata are
possible in this example because of the small number of parameters and the simple analytical model involved (equation(9).
We emphasize that the value of the goal function during
each iteration can be adjusted by a tolerance, g , tolerance,

allowing us to choose prior uncertainty domains for the next
iteration. Table 3 demonstrates this point for iteration 1 inTable 2 and for different tolerance levels. As the magnitude of
the tolerance is decreased, the number of hits decreases until
there are no hits in any of the strata for each parameter. Our
experience indicates that it is better to be relatively conservative. For example, in the present case it is better to use for thesecond iteration the domain having a tolerance of 0.06 rather
than domains having tolerances of 0.04 or 0.02. Having ahigher tolerance will avoid the possibility of falling into local
minima as the iterative process continues. The data in Table 3
also reveal that for the given uncertinty domains, the goal
function is far more sensitive to the retardation factor R than to
the Peclet number P since the strata for R are much better discriminated (leading to narrower intervals with hits) than for P.
In a second run we assumed that each measured datum
could be in error by 50%. In other words, the actual value of

Figure 5. A comparison of measured water contents (datapoints) with simulated values using parameters obtained withthe three-parameter method 1 approach of van Genuchten etal. [1987] (dashed line) and the three-parameter (u only) SUFI
scenario 1 approach (solid line).

Figure 6. A comparison of measured water contents (datapoints) with simulated values using parameters obtained withthe three-parameter method 1 approach of van Genuchten etal. [1987] (dashed line) and the three-parameter (u and h data)
SUFI scenario 2 approach (solid line).

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

1885

Example 2: Estimating the Unsaturated
Soil Hydraulic Properties
The second example involves transient drainage of waterfrom a large caisson 6 m deep and 3 m in diameter. The
experimental study was conducted by Abeele [1984] at Los
Alamos National Laboratory for the purpose of measuring theunsaturated soil hydraulic properties of Bandelier Tuff in connection with studies of radionuclide migration from waste disposal facilities. The initially dry tuff in the caisson, instrumented with tensiometers and neutron probe access tubes, wasfirst saturated with water and then allowed to drain undergravity for 100 days. The experimental data were previously
analyzed by van Genuchten et al. [1987] and Kool et al. [1987]
using nonlinear least squares optimization methods based on
the Levenberg-Marquardt method [Marquardt, 1963]. The data
were used to estimate the hydraulic parameters in the following model for the unsaturated soil hydraulic parameters [vanGenuchten, 1980]:

u 5 ur 1
Figure 7. A comparison of measured water contents (datapoints) with simulated values using parameters estimated withthe six-parameter method 2 approach of van Genuchten et al.
[1987] (dashed line) and the five-parameter (u only) SUFI
scenario 3 approach (solid line).

iieach x m
was uniformly distributed in the interval [0.5x m,
1.5x im ], where i 5 1, z z z , N. The results of this example areshown in Table 4. We notice the larger values of the tolerances
as compared with the previous case and the inability to reduce
the tolerance after the third run. In general, measurement
errors have the effect of producing more hits, and this limits
the resolution with which we can choose the posterior interval.
The overall effect of the measurement errors, based on our
definition of a hit, is to give a larger uncertainty interval foreach parameter. Given the 50% error in the measurements, wecannot get more precise values for the parameters than theintervals (5.0, 20.75) for P and (1.3, 1.58) for R. This constitutes the uncertainty in these parameters after fitting. Compare these values with (19.659, 19.665) for P and (1.34766,
1.34772) for R obtained in the case without measurement
errors.

Table 6.
Numberof
Iterations
1
2
3
4
5

us 2 ur~1 1 uahu n! m

(13)

m 2
K 5 K sS le@1 2 ~1 2 S 1/m
e ! #

(14)

where u is volumetric water content, u r is residual water content, u s is saturated water content, h is pressure head, a and n
are shape factors, m 5 1 2 1/n, K s is saturated hydraulicconductivity, l is a pore connectivity parameter, and S e iseffective fluid saturation:
S e~ u ! 5

u 2 urus 2 ur(15)

The drainage data were analyzed by van Genuchten et al.
[1987] and Kool et al. [1987] using two different methods withdifferent goal functions and input data. In the first method(method 1) three parameters, a , n, and K s , were estimatedusing the following goal function

OO
M

g~B! 5

P

i51 j51

OP@ u *ij 2 u ij~B!# 2 1 W 2@h *j 2 h j~B!# 2

(16)

j51

where u *ij represents measured water contents at M 5 5depths and P 5 6 different times, h *j is the measured pressure
head at x 5 0.4 m, u ij (B) and h j (B) are model-predicted u andProgression to Convergence for Example 2, Scenario 2
Parameter Value
Parameters
21

a, cm
n
K s , cm d21
a, cm21n
K s , cm d21
a, cm21n
K s , cm d21
a, cm21n
K s , cm d21
a, cm21n
K s , cm d21

Minimum

Maximum

Number ofStrata

0.005
1.0
1.0
0.00951.3
10.9
0.01081.39
10.9
0.01081.474
10.9
0.0112
1.54
11.2

0.05
2.0
100.0
0.023
1.6
60.4
0.01761.6
30.7
0.01421.6
18.8
0.0125
1.6
14.2

10
10
10
10
10
10
10
10
20
10
10
30
10
10
30

Number of Hits in Each Stratum

Tolerance

0, 1, 1, 2, 0, 0, 0, 0, 0, 0
0, 0, 0, 2, 1, 1, 0, 0, 0, 0
0, 1, 1, 0, 1, 1, 0, 0, 0, 0
0, 3, 4, 4, 2, 1, 0, 0, 0, 0
0, 0, 0, 1, 1, 2, 2, 3, 2, 2
6, 4, 2, 1, 0, 0, 0, 0, 0, 0
1, 6, 5, 3, 1, 0, 0, 0, 0, 0
0, 0, 0, 0, 1, 1, 4, 4, 3, 3
2, 3, 4, 2, 2, 2, 0, 1, 12 3 0
0, 0, 1, 1, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 0, 1, 1, 0
0, 0, 0, 0, 1, 0, 1, 23 3 0
0, 0, 0, 0, 0, 1, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 1, 0, 0, 0
12 3 0, 1, 17 3 0

0.01
0.007
0.00540.005
0.004945

1886

Table 7.
Numberof
Iterations
1

2

5

6

7

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Progression to Convergence for Example 2, Scenario 3
Parameter Value
Parameters

Minimum

Maximum

Number ofStrataurus
a , cm21n
K s , cm d21
urus
a , cm21n
K s , cm d21
urus
a , cm21n
K s , cm d21
urus
a , cm21n
K s , cm d21
urus
a , cm21n
K s , cm d21

0.001
0.3
0.005
1.0
1.0
0.001
0.36
0.005
1.4
18.0
0.00920.36750.005421.75
28.0
0.03670.369
0.0056
1.82
28.0
0.039
0.369
0.0057
1.86
28.0

0.1
0.4
0.05
2.0
100.0
0.1
0.4
0.02752.0
83.0
0.05050.375
0.006251.9
34.0
0.05050.373
0.00625
1.9
34.0
0.046
0.371
0.0061
1.89
30.0

4
4
6
6
6
4
4
6
6
6
6
4
4
6
6
6
4
4
6
6
6
4
4
6
6

h values corresponding to parameter vector B 5 ( a , n, K s ),and W 5 0.0016 is a weighting coefficient chosen such thatthe two terms of (16) attain roughly the same value [Kool etal., 1987]. The estimation method involved an inverse solution of the one-dimensional unsaturated flow (Richards)
equation for an initially saturated soil profile, i.e., h( x) 5 0,
subject to a no flow boundary condition at x 5 0 and a

Number of Hits in Each Stratum

Tolerance

1, 1, 2, 1
0, 0, 3, 2
3, 1, 1, 0, 0, 0
0, 0, 0, 2, 2, 1
0, 1, 2, 1, 1, 0
3, 2, 1, 0
1, 2, 3, 0
5, 1, 0, 0, 0, 0
0, 1, 2, 0, 2, 1
0, 3, 2, 0, 1, 0
0, 0, 0, 0, 4, 4
0, 4, 4, 0
0, 1, 4, 30, 0, 0, 1, 4, 3
4, 3, 1, 0, 0, 0
0, 1, 0, 2, 1, 0
1, 2, 1, 0
0, 1, 3, 0
0, 0, 0, 2, 2, 0
4, 0, 0, 0, 0, 0
0, 0, 0, 1, 0, 0
0, 1, 0, 0
0, 0, 1, 0
0, 1, 0, 0, 0, 0
1, 0, 0, 0, 0, 0

0.00850.008

0.007290.007258

0.007238

first-type lower boundary condition at x 5 4.23 m by linearly
extrapolating between observed pressure head data at that
depth.
In the second method (method 2), six parameters, B 5 {K s ,a , n, u r , u s , l}, were estimated directly from reported laboratory and caisson u (h) and K(h) data [Abeele, 1979, 1984]
using the following goal function

Figure 8. A comparison of measured pressure heads (data points) with simulated values using parameters
estimated by method 1 of van Genuchten et al. [1987] (dashed line) and scenarios 1, 2, and 3 of SUFI (solidlines).

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Figure 9. Hydraulic properties of Bandelier Tuff. (a) Water
retention u (h) and (b) hydraulic conductivity K( u ) curves areshown using different methods. Data points are measured values from Abeele [1984].

OM

g~B! 5i51

ON@ u *i 2 u i~B!# 2 1 V$log ~K *j ! 2 log @K j~B!#% 2

j51

(17)
where u *i and u i (B) are observed and predicted water contents
at M values of the pressure head, K *j and K j (B) are observed
and predicted conductivities at N values of h, and V is aweighting factor that makes the two terms of (17) of roughly
equal value. The results obtained by van Genuchten et al.
[1987] for the two runs are shown in Table 5.
We ran three different scenarios with SUFI. First, the goal
function (2) with observed water contents only was used toperform a three-parameter ( a , n, and K s ) estimation. Second,
the goal function (16) was used with measured water contents
and pressure heads to again estimate a , n, and K s . Third, a1887

five-parameter ( u r , u s , a , n, and K s ) estimation was carriedout, again with goal function (2) and again using only watercontents as the conditioning data. The latter scenario was used
to examine the observation by Kool et al. [1987] that simultaneous estimation of three or more parameters required more
information than just water content data. The HYDRUS [Vogel et al., 1996] code was used in each case to obtain simulated
water content and pressure head data. Initial and boundary
conditions were similar to those used by Kool et al. [1987].
Hydraulic parameters estimated by SUFI using the three
scenarios above are also listed in Table 5 and compared withthe results obtained by van Genuchten et al. [1987] (see also
Kool et al., 1987). Figures 5, 6, and 7 show the measured watercontents as a function of depth after 1, 4, 20, and 100 days for
scenarios 1, 2, and 3, respectively. For scenarios 1 and 2 thecurves based on the three-parameter (method 1) estimation
results obtained by van Genuchten et al. [1987] are also shown
for reference, while for scenario 3 (Figure 7) curves based on
the six-parameter method 2 estimation are shown. In all cases
the parameters estimated by SUFI lead to a closer fit to thedata. Tables 6 and 7 show the progression of SUFI toward
convergence for scenarios 2 and 3, respectively. We note herethat the measured value of K s as obtained by Abeele [1984] was
12.4 cm d21.
Two important observations can be made for this test case.
First, contrary to the observation of Kool et al. [1987] withrespect to their nonlinear least squares analysis, SUFI wascapable of a five-parameter inversion using only observed water content values. Actually, SUFI was found to have no limitation in this regard, and any number of parameters can be
fitted using any pertinent observed data. As previously mentioned, the only limitation of SUFI may be the total number ofsimulations that need to be carried out as a function of thestratification of the parameters and the speed of the simulation
code. Second, Figures 5, 6, and 7 seem to suggest that the best
fit was obtained with the five-parameter case (u data only)
followed by the three-parameter case without h data, and finally the three-parameter case using both u and h data. However, this conclusion is correct only if we consider water content as the output. The five-parameter estimation is highly
conditioned on the 30 observed water content data and provides the best fit because of more degrees of freedom in thefitting process. This conditioning on only one variable, in thiscase water content, may not necessarily provide better estimates for other output variables such as pressure heads, drainage rates, or concentrations if measured. This point is further
illustrated in Figures 8, 9a, and 9b.
Figure 8 compares the observed pressure head data withpredicted values when the parameter sets of scenarios 1, 2, and3, as well as those of method 1 of van Genuchten et al. [1987](Table 5), are used in the forward problem. Clearly, goal functions which included observed pressure head values provided
the best estimates. The same is evident when the calculated
hydraulic properties (equations (13) and (14)) using SUFIderived parameters (scenarios 1, 2, and 3) are compared withthe measured hydraulic data (data points) and the method 2
curves [van Genuchten et al., 1987] fitted to those data (Figures
9a and 9b). The data points in Figures 9a and 9b are a composite of data derived from the caisson drainage experimentusing a standard instantaneous profile (gravity drainage) analysis in the wet range (u . 0.164) and independently measured
retention and conductivity data in the dry range (u , 0.164)
obtained in the laboratory using small cores [Abeele, 1979,

1888

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Figure 10. Schematic of the system considered in test example 3. Points indicate locations of the 14
concentration sampling points.

1984]. Notice that the simultaneous use of observed pressure
heads and water contents in the caisson (scenario 2) now provides a much better fit than when only water content data areused (scenarios 1 and 3). These results indicate that in general,
it is best to somehow include data of different types in the goal
function for an all-purpose simulation. Such an approach requires selection of proper coefficients, such as W in (17), forthe different data types used in the goal function. As discussed
in more detail later, this task is not a trivial one.
Example 3: Hypothetical Water and Solute Infiltration Test
The third example considered here is a hypothetical caseinvolving the radially symmetric three-dimensional infiltrationof water and a dissolved solute from a single-ring infiltrometer
into a variably saturated, layered soil profile (Figure 10). This
problem was previously used to mathematically verify part of
the SWMS 2D code version 1.12 [Simunek et al., 1994]. The
layered soil profile was assumed to be that of the Hupselse
Beek watershed in the Netherlands, consisting of a 40-cm thickA horizon on top of a 300-cm B/C horizon. The hydraulicparameters describing the observed mean-scaled hydraulicfunctions and the assumed transport parameters of the two soil
layers [Cislerova, 1987; Hopmans and Stricker, 1989] are given
in Table 8. For our purposes these parameters are treated here
as true parameters.
Calculations were carried out for a period of 5 days using the
SWMS 2D code with the true parameters as listed in Table 8.
The soil profile was assumed to be initially free of any solute.
All sides of the flow region were considered to be impervious
except for a ponded region (h 5 0) inside the ring infiltrometer around the origin at the surface (Figure 10). Figure 11

presents calculated pressure head and concentration profiles at
two different times. Calculated concentrations after 2 days at
14 arbitrary locations (Figure 10) were chosen to represent themeasured values.
We assumed that six of the hydraulic and solute parametersin Table 8 were unknown: the saturated hydraulic conductivities of the two soil layers (K s1 and K s2 ), the longitudinal andtransverse dispersivity coefficients (D L and D T ), the adsorption coefficient (k), and a first-order rate constant for solutedegradation in the liquid phase ( m w ). The adsorption coefficient appears in the retardation factor R as follows
Table 8. Soil Hydraulic Parameters of the Hupselse Beek
Area Used in the Verification Example of the
SWMS 2D CodeHydraulic andTransport
Parameters

True Values Used in
SWMS 2D Forward
Simulation

Values Obtained by
SUFILayer 1

Layer 2

Layer 1

Layer 2usura , cm21n
K s , cm d21
D L , cm
D T , cmk cm3 g21m w , min21

0.399
0.00010.0174
1.37570.02070.5
0.1
0.1
23.472E-5

0.339
0.00010.01391.60240.03150.5
0.1
0.1
23.472E-5

0.399*
0.0001*0.0174*
1.3757*0.02090.500025
0.101
0.103223.625E-5

0.339*
0.0001*0.0139*1.6024*.0316660.500025
0.101
0.103223.625E-5

Read 23.472E-5 as 23.472 3 1025.
*Parameters assumed to be known.

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Figure 11. Simulated (left) pressure head profiles and (right)
concentration profiles at (top) t 5 0.25 and (bottom) t 5 5days for test example 3.

R511

rk
u

(18)

where r is the soil bulk density. The dispersion, adsorption,and degradation parameters were assumed to be the same forthe two soil layers. We subsequently used SUFI with (2) as thegoal function to estimate the 6 unknowns using the 14 “measured” concentration values as the conditioning information.
The initial parameter values used to start SUFI were assumed to contain large uncertainties and were depicted with

Table 9.

Progression to Convergence for Test Example 3

Numberof
Iterations

Parameter Value

1

3

4

8

uniform distributions within the following intervals: K s1 5
[0.001, 0.1], Ks2 5 [0.001, 0.1], DL 5 [0.15, 0.75], DT 5 [0.01, 0.4],
k 5 [0.05, 0.5], and m w 5 [21.0 3 10 24 , 21.0 3 10 26 ].
The final estimates of these parameters as obtained with SUFI
using the exhaustive stratified sampling option are given in
Table 8. Table 9 shows the progression to convergence. The
predicted concentration distributions using the final estimates
were essentially identical to those shown in Figure 11.
Since the true parameters for this hypothetical example were
known, the simulation model should be able to fully explain themeasured concentrations; that is, we should be able to achievea value of zero for the goal function. Such a scenario is highly
unlikely in reality since measurements are always subject tosome error. Also, the selected simulation model with its limited
number of invoked processes and associated parameters will
hardly ever correctly describe the true physical processes operating in the subsurface, thus leading to additional approximations or inaccuracies for the estimated parameters in thatsimulation model. Numerical approximations in the simulation
model itself may also further limit the final accuracy of theparameter estimation process. To obtain a goal function ofzero in this example, SUFI requires that the true values of thesix parameters all fall in the middle of an interval in theirrespective strata. The data shown in Table 9 suggest that a goal
function of zero eventually would be obtained if the iterative
process were continued. Since having a goal function of zero is
not the objective, the iterative process was terminated after theeighth iteration when the tolerance reached a very small value.
An interesting observation in Table 9 is that during the firstand second iterations, SUFI could not discriminate among thedifferent strata of D L and D T despite the relatively large uncertainty domains used for these two parameters. This result
indicates that the goal function was insensitive to these parameters for the given level of uncertainty in the other parameters.
When the other parameters became more defined during theParameters

Minimum

Maximum

Number ofStrataD L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21

0.15
0.01
0.05
21.0E-4
0.001
0.001
0.25
0.01
0.05
21.0E-4
0.017
0.025
0.45
0.01
0.087
25.5E-5
0.02
0.025
0.45
0.01
0.087
24.0E-5
0.02
0.025

0.75
0.4
0.5
21.0E-6
0.1
0.1
0.65
0.4
0.16
21.0E-5
0.026
0.05
0.55
0.205
0.123
21.0E-5
0.02150.037
0.51670.14
0.105
21.0E-5
0.021
0.033

2
2
2
2
2
2
4
2
3
2
2
2
3
3
4
3
3
3
2
5
5
4
5
3

Read 21.0E-4 as 21.0 3 1024.

1889

Number of Hits in
Each Stratum
2, 2
2, 2
4, 0
4, 0
4, 0
4, 0
0, 0, 1, 0
1, 00, 1, 00, 10, 0, 1, 0, 0, 01, 01, 1, 01, 1, 01, 1, 00, 2, 01, 1, 01, 1, 00, 10, 0, 0, 1, 0
0, 0, 0, 0, 1
1, 0, 0, 0
0, 0, 0, 0, 1
0, 0, 1

Tolerance
0.04

0.005

0.003

0.00021890

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

Table 10. Progression to Convergence for Test Example 3 Using Random Stratified Sampling Procedure
Numberof
Iterations
1

3

4

10

Parameter Value
Parameters

Minimum

Maximum

Number ofStrataD L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21
D L , cm
D T , cmk , cm3 g21
m W , d21
K s1 , cm d21
K s2 , cm d21

0.15
0.01
0.05
21.0E-4
0.001
0.001
0.32
0.075
0.05
21.0E-4
0.01750.0285
0.32
0.075
0.05
21.0E-4
0.01750.02850.498
0.09850.097
23.5E-50.020540.0313

0.75
0.4
0.5
21.0E-6
0.1
0.1
0.65
0.4
0.5
21.0E-6
0.045
0.056
0.65
0.18330.2
21.17E-5
0.02670.0376
0.515
0.106
0.103
23.2E-50.02090.0317

6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6

third iteration, D L and D T started to exert more effect on thegoal function. This example hence shows that the sensitivity ofan objective function to a certain parameter may depend on
the invoked or calculated uncertainty domains of the otherparameters. The example also indicates that SUFI is a very
appropriate tool for carrying out sensitivity analyses. Most
traditional methods of carrying out sensitivity analyses invoke
systematic changes in a certain parameter of interest while
keeping all other parameters constant. Such an approach as-

Number of Hits in Each Stratum
0, 6, 9, 4, 8, 0
5, 0, 3, 8, 0, 10
0, 9, 8, 0, 0, 7
4, 4, 0, 8, 0, 8
4, 21, 12, 0, 0, 0
3, 5, 6, 9, 3, 0
9, 0, 6, 0, 0, 5
1, 9, 8, 0, 0, 0
1, 8, 0, 0, 0, 0
12, 7, 0, 0, 3, 0
1, 8, 0, 0, 0, 0
7, 4, 0, 0, 0, 0
0, 9, 3, 4, 0, 0
4, 7, 0, 0, 0, 0
4, 0, 0, 11, 0, 0
0, 0, 0, 4, 11, 4
4, 0, 0, 13, 0, 0
0, 3, 7, 3, 0, 3
2, 2, 2, 1, 2, 5
2, 1, 3, 1, 5, 1
1, 0, 2, 10, 0, 0
0, 0, 1, 1, 5, 5
0, 1, 0, 2, 5, 4
1, 3, 3, 3, 3, 0

Tolerance
0.08

0.012

0.00780.0005sumes exact (true) knowledge of the other parameters andhence is not appropriate for cases where the true values ofparameters are not known. By contrast, SUFI provides amethod for situations where the parameters are treated withuncertainty.
The results with the random stratified sampling which was
also performed for this example are shown in Table 10. In thisexample we divided the uncertainty intervals of the parameters
into six strata. This produces an exhaustive sample space of
46,656 combinations. We allowed only 0.3% of the combinations to be simulated, and this resulted in about 140 simulations in each iteration. We stopped after the tenth iteration,
where a goal value of 0.0005 was obtained. The random stratified sampling also worked well for this example, and on thebasis of our experience it is better to use this sampling scheme
when many parameters are involved.

Influence of the Goal Function

Figure 12. A comparison of measured effluent concentrations (data points) with simulated values using parameters
estimated by SUFI assuming a logarithm-based root-meansquare error (LRMSE) objective function and a regular rootmean-square error (RMSE) objective function.

A detailed analysis of different goal functions is beyond thescope of this paper. Still, we would like to show briefly thepotentially important implications of formulating and using
different goal functions. The implications are illustrated for the
first test example which we analyzed also using the logarithmbased (LRMSE) goal function as given by (3). In this case we
obtained the results P 5 57.16 and R 5 1.153. Figure 12
shows that the LRMSE-fitted effluent curve differs substantially from the previously obtained RMSE curve. The pronounced differences between the two cases is caused by a scale
phenomenon. When the RMSE is used, the contributions ofrelatively small concentrations to the value of the goal function
are essentially ignored. By comparison, LRMSE puts moreweight on the smaller concentrations at the expense of thehigher values.
Proper definition of the goal function can become evenmore difficult, and somewhat subjective, when different types

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDUREof data are represented in the goal function, such as was thecase with (16) for the caisson drainage example 2. When we
repeated the second scenario for this example using W rather
than W 2 in (16), the estimates for the three unknowns became
quite different: a 5 0.0224 (cm21), n 5 1.12576, and K s 5
85.5 (cm d21). Given that the measured K s was 12.4 cm d21,
the parameter estimates using W should be relatively inaccurate. Using W instead of W 2 in this case leads to a goal function
whose value is dominated by pressure heads solely because of
their numerically larger values; this occurs in spite of the factthat in situ pressure head measurements often yield less accurate data than water content measurements.

Conclusions
The acronym SUFI was coined to represent a parameter
estimation procedure which is sequential, operates within parameter uncertainty domains, employs only forward calculations, and is iterative in nature. The procedure is general,
stable, always convergent, and has no inherent limitations interms of the number of parameters that can be consideredsimultaneously. The number of function calls can be limited by
invoking a less exhaustive sampling scheme, or the program
can be made to run faster by parallel submission of smallerjobs. The sequential uncertainty domain parameter estimation
procedure performed well in three different examples involving a wide range of flow and transport parameters. The procedure begins with a prior and concludes with a posterior stateof belief of the value of the parameters being estimated. The
states of belief are expressed as uncertainty domains, with theposterior uncertainty domain being much smaller than theprior one. Although not discussed in great detail, the importance of the goal function on the results is demonstrated.

References
Abbaspour, K. C., R. Schulin, E. Schläppi, and H. Flühler, A Bayesian
approach for incorporating uncertainty and data worth in environmental projects, Environ. Model. Assess., 1, 151–158, 1996.
Abeele, W. V., Determination of hydraulic conductivity in crushed
Bandelier Tuff, Tech. Rep. LA-8147-MS, Los Alamos Natl. Lab., Los
Alamos, N. M., 1979.
Abeele, W. V., Hydraulic testing of crushed Bandelier Tuff, Tech. Rep.
LA-10037-MS, Los Alamos Natl. Lab., Los Alamos, N. M., 1984.
Benjamin, J. R., and C. A. Cornell, Probability, Statistics, and Decision
for Civil Engineers, McGraw-Hill, New York, 1970.
Carrera, J., State of the art of the inverse problem applied to the flow
and solute transport equations, in Groundwater Flow and Quality
Modeling, NATO ASI Ser. Ser. C, vol. 224, edited by E. Curtodio, A.
Gurgui, and J. P. Lobo Ferreira, pp. 549 –583, D. Reidel, Norwell,Mass., 1987.
Carrera, J., and S. P. Neuman, Estimation of aquifer parameters undertransient and steady state conditions, 1, Maximum likelihood
method incorporating prior information, Water Resour. Res., 22,
199 –210, 1986.
Cislerova, M., Comparison of simulated water balance for ordinaryand scaled soil hydraulic characteristics, Tech. Rep. 82, Dep. of
Hydraulics and Catchment Hydrol., Wageningen Agric. Univ.,
Wageningen, Netherlands, 1987.
Clifton, P. M., and S. P. Neuman, Effects of kriging and inverse
modeling on conditional simulation of the Avra Valley aquifer insouthern Arizona, Water Resour. Res., 18, 1215–1234, 1982.
Cooley, R. L., Incorporation of prior information on parameters into
nonlinear regression groundwater flow models, 1, Theory, Water
Resour. Res., 18, 965–967, 1982.
Cooley, R. L., Incorporation of prior information on parameters into
nonlinear regression groundwater flow models, 2, Applications, Water Resour. Res., 19, 662– 676, 1983.
Dagan, G., Stochastic modeling of groundwater flow by unconditional1891

and conditional probabilities: The inverse problem, Water Resour.
Res., 21, 65–72, 1985.
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 head data, Soil Sci.
Soc. Am. J., 57, 1167–1175, 1993.
Eching, S. O., J. W. Hopmans, and W. W. Wallender, Estimation of insitu unsaturated soil hydraulic functions from scaled cumulativedrainage data, Water Resour. Res., 30, 2387–2394, 1994.
Hoeksema, R. J., and P. K. Kitanidis, An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling, Water Resour. Res., 20, 1003–1020, 1984.
Hopmans, J. W., and J. N. M. Stricker, Stochastic analysis of soil water
regime in a watershed, J. Hydrol., 105, 57– 84, 1989.
Kitanidis, P. K., Quasi-linear geostatistical theory for inversing, Water
Resour. Res., 31, 2411–2419, 1995.
Kitanidis, P. K., and E. G. Vomvoris, A geostatistical approach to theinverse problem in groundwater modeling (steady state) and onedimensional simulations, Water Resour. Res., 19, 677– 690, 1983.
Kool, J. B., and J. C. Parker, Analysis of the inverse problem fortransient unsaturated flow, Water Resour. Res., 24, 817– 830, 1988.
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.
Marquardt, D. W., An algorithm for least squares estimation of nonlinear parameters, SIAM J. Appl. Math., 11, 431– 441, 1963.
McKay, M. D., R. J. Beckman, and W. J. Conover, A comparison of
three methods for selecting values of input variables in the analysis
of output form a computer code, Technometrics, 21(2), 239 –245,
1979.
Mishra, S., and J. C. Parker, Parameter estimation for coupled unsaturated flow and transport, Water Resour. Res., 25, 385–396, 1989.
National Research Council (NRC), Ground Water Models, Scientificand Regulatory Applications, 303 pp., Comm. on Groundwater
Model. Assess., Water Sci. and Technol. Board, Natl. Res. Counc.,
Natl. Acad. Press, Washington, D. C., 1990.
Parker, J. C., and M. T. van Genuchten, Determining transport parameters from laboratory and field trace experiments, Va. Agric. Exp.
Stn. Bull. 84, 92 pp., 1984.
Parker, J. C., J. B. Kool, and M. T. van Genuchten, Determining soil
hydraulic properties from one-step outflow experiments by parameter identification, 2, Experimental studies, Soil Sci. Soc. Am. J., 49,
1354 –1359, 1985.
Ségol, G., Classic Groundwater Simulations; Proving and Improving
Numerical Models, 531 pp., Prentice-Hall, Englewood Cliffs, N. J.,
1994.
Shouse, P. J., J. B. Sisson, T. R. Ellsworth, and J. A. Jobes, Estimatingin situ unsaturated hydraulic properties of vertical heterogeneous
soils, Soil Sci. Soc. Am. J., 56, 1673–1679, 1992.
Simunek, J., T. Vogel, and M. T. van Genuchten, The SWMS 2D Codefor Simulating Water Flow and Solute Transport in Two-Dimensional
Variably Saturated Media, Version 1.21, 197 pp., U.S. Salinity Lab.,
Agric. Res. Serv., U.S. Dep. of Agric., Riverside, Calif., 1994.
Sisson, J. B., and M. T. van Genuchten, An improved analysis ofgravity drainage experiments for estimating the unsaturated soil
hydraulic functions, Water Resour. Res., 27, 569 –575, 1991.
Sun, N. Z., and W. W.-G. Yeh, Coupled inverse problems in groundwater modeling, 1, Sensitivity analysis and parameter identification,
Water Resour. Res., 26, 2507–2525, 1990.
Van Dam, J. C., J. N. M. Stricker, and P. Droogers, From one-step tomulti-step: Determination of soil hydraulic functions by outflow
experiments, Tech. Rep. 7, Dep. of Hydrol., Soil Phys., and Hydraulics, Wageningen Agric. Univ., Wageningen, Netherlands, 1990.
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., and P. J. Wierenga, Solute dispersion coefficients and retardation factor, in Methods of Soil Analysis, 1, Physicaland Mineralogical Methods, Agron. Monogr. Ser., vol. 9, edited by A.
Klute, 2nd ed., pp. 1025–1053, Am. Soc. of Agron. and Soil Sci. Soc.
of Am., Madison, Wis., 1986.
van Genuchten, M. T., J. C. Parker, and J. B. Kool, Analysis andprediction of water and solute transport in a large lysimeter, in
Modeling Study of Solute Transport in the Unsaturated Zone, Workshop Proceedings, Los Alamos National Laboratory, NUREG/CR-1892

ABBASPOUR ET AL.: A SEQUENTIAL UNCERTAINTY DOMAIN INVERSE PROCEDURE

4615, LA-10730-MS, vol. 2, edited by E. P. Springer and H. R.
Fuentes, pp. 4 –31, U.S. Nucl. Regul. Comm., Washington, D. C.,
1987.
Vogel, T., K. Huang, R. Zhang, and M. T. van Genuchten, The
HYDRUS code for simulating one-dimensional water flow, solutetransport, and heat movement in variable-saturated media, version
5, U.S. Salinity Lab., Agric. Res. Serv., U.S. Dep. of Agric., Riverside, Calif., 1996.
Yeh, W. W.-G., Review of parameter identification procedures inground water hydrology: The inverse problem, Water Resour. Res.,22, 95–108, 1986.
Yeh, W. W.-G., and Y. S. Yoon, Aquifer parameter identifiability with
optimum dimension in parameterization, Water Resour. Res., 17,664 – 672, 1981.
Yeh, T.-C. J., J. Minghui, and S. Hanna, An iterative stochastic inverse

method: Conditional effective transmissivity and hydraulic head
fields, Water Resour. Res., 32, 85–92, 1996.
K. C. Abbaspour and R. Schulin, Department of Soil Protection,
Swiss Federal Institute of Technology, Grabenstrasse 3, CH-8952
Schlieven, Switzerland. (e-mail: abbaspou@baikal.ethz.ch)
E. Schläppi, Columbi Schmutz Dorthe AG., Konsumstrasse 20, 3007
Bern, Switzerland.
M. T. van Genuchten, U.S. Salinity Laboratory, Agricultural Research Service, U.S. Department of Agriculture, 450 West Big Spring
Road, Riverside, CA 92507.

(Received February 12, 1997; accepted April 23, 1997.)