Estimating Uncertain Flow and Transport Parameters Usinga Sequential Uncertainty Fitting Procedure
K. C. Abbaspour,* C. A. Johnson, and M. Th. van Genuchten

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT

and Erofeeva, 2002), and seawater intrusion (Iribar etal., 1997).
Inverse modeling has brought new opportunities, butalso new challenges. Some of the advantages are savingsin time and cost of laboratory and/or field experiments
generally needed to obtain the unknown parameters andattainment of a better fit with available data. Anotheradvantage is the usefulness of inverse modeling in theanalysis of model structure (the invoked conceptualmodel), boundary conditions, and prevailing subsurface
flow and contaminant transport processes. One of thelimitations of inverse modeling is that the fitted parameters are conditioned on the experimental setup and its
set of measured variables, which is usually limited inboth time and space. Other conditioning factors include
the choice of the parameter estimation routine, the objective function, and the weights associated with thedifferent components or parts of the form of objective function.
The main challenge of inverse modeling is to find that“unique” set of effective parameters that best describe
the data and envisioned processes. However, many stochastic inverse methods are based on the philosophy
that there is no such unique set of effective parametervalues. Beven and Binley (1992), among others, explainthat since the invoked hydrological model structure usually contains some error, and since all observations andmeasurements on which inverse modeling is based arealso subject to error, no reason exists why any one setof parameter values will represent the sought after effective parameters in a unique way. Hence, the nonuniqueness problem is an inherent characteristic of inverse
modeling regardless of the calibration procedure or experimental setup.
Furthermore, with field-scale time series data, as opposed to relatively simple well-defined laboratory solute
breakthrough curves for example (e.g., van Genuchten,1981), definition of a good simulation can be very subjective because visually different simulations with very
different parameter combinations may yield identicalvalues of the objective function (Boyle et al., 2000). In
such cases, parameter fitting leads to parameter regionswith similar values of the objective function. These parameter regions comprise the uncertainty in the parameters. Quantifying the degree of uncertainty in the parameters leads to more reliable models (Cooley, 1993).
Uncertainty regions of potential candidates for effective
parameter values can be used to produce prediction
uncertainties (Konikow, 1986, p. 183). Or vice versa, as
is done in our study, if a desired prediction uncertainty

Inversely obtained hydrologic parameters are always uncertain(nonunique) because of errors associated with the measurements andthe invoked conceptual model, among other factors. Quantification ofthis uncertainty in multidimensional parameter space is often difficult
because of complexities in the structure of the objective function. In
this study we describe parameter uncertainties using uniform distributions and fit these distributions iteratively within larger absolute intervals such that two criteria are met: (i) bracketing most of the measured
data (⬎90%) within the 95% prediction uncertainty (95PPU) and (ii)obtaining a small ratio (⬍1) of the average difference between theupper and lower 95PPU and the standard deviation of the measured
data. We define a model as calibrated if, upon reaching these two
criteria, a significant R2 exists between the observed and simulated
results. A program, SUFI-2, was developed and tested for the calibration of two bottom ash landfills. SUFI-2 performs a combined optimization and uncertainty analysis using a global search procedure andcan deal with a large number of parameters through Latin hypercube
sampling. We explain the above concepts using an example in which
two municipal solid waste incinerator bottom ash monofills were successfully calibrated and tested for flow, and one monofill also fortransport. Because of high levels of heavy metals in the leachate,monitoring and modeling of such landfills is critical from environmental points of view.

Inverse modeling is becoming increasingly popularin many branches of the earth and environmental
sciences. One major reason to apply inverse modeling is
to estimate parameters that cannot be measured directly
because of a variety of reasons, including scale issues.
What makes inverse modeling attractive is that various
output variables are often much easier to measure thansome of the required input parameters. Numerous applications of inverse modeling now exist in the literature.
These include estimation of soil hydraulic parameters(e.g., Pan and Wu, 1998; Simunek et al., 1999; Abbaspour et al., 2001; Young et al., 2002; Wang et al., 2003),
hydrologic parameters (e.g., Duan et al., 1992; Bevenand Binley, 1992; Lahkim et al., 1999; Savenkoff et al.,2001), and transport parameters (e.g., Schmied et al., 2000;
von Gunten and Furrer, 2000; McKenna et al., 2001;
Bell et al., 2002), as well as parameters in many otherapplications such as flood discharge prediction (Sulzer
et al., 2002), barotropic ocean tide simulations (Egbert
K.C. Abbaspour and C.A. Johnson, Swiss Federal Inst. for Environmental Science and Technology (EAWAG), Ueberlandstrasse 133,
P.O. Box 611, CH-8600 Dübendorf, Switzerland; M. Th. van Genuchten, USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 West Big
Springs Road, Riverside, CA 92507. Received 12 Dec. 2003. OriginalResearch Paper. *Corresponding author (abbaspour@eawag.ch).

Abbreviations: EC, electrical conductivity; GLUE, Generalized Likelihood Uncertainty Estimation; MSWI, municipal solid waste incinerator; SSQ, sum of square errors; SUFI, Sequential Uncertainty FItting;
95PPU, 95% prediction uncertainty.

Published in Vadose Zone Journal 3:1340–1352 (2004).
© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

1340

1341

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

level is obtained, all parameter distributions producing
this uncertainty are potential solutions. In this paper,
we outline a procedure that uses a stopping rule basedon a certain definition of the prediction uncertainty.
Inverse modeling usually reduces to an optimization
problem in which a vector of unknown parameters is estimated by minimizing an assumed goal or objective function. Many procedures have been developed to solve theresulting nonlinear minimization problem (Brad, 1974;
Beck and Arnold, 1977; Yeh, 1986; Kool et al., 1987;
Duan et al., 1992; Yapo et al., 1998; Vrugt et al., 2003a,b).
Among soil scientists and vadose zone hydrologists thegradient approach of Levenberg-Marquardt (Marquardt,
1963) method has become the method of choice (e.g.,van Genuchten 1981; Kool et al., 1985, 1987), whilemodelers of larger-scale problems prefer using moreglobal approaches and multiobjective optimization functions (e.g., Beven and Binley, 1992; Duan et al., 1992,
2003; Yapo et al., 1998; Gupta et al., 1998).
A common problem with gradient methods is thatthey are very sensitive to the initial values of the optimized parameters and often lead to a local minimum.
Also, their application to problems having a large number of parameters is sometimes problematic. Moreover,gradient methods do not provide a reliable estimate ofparameter uncertainty and, hence, prediction uncertainty. Global procedures such as the shuffled complex(Duan et al., 1992) and genetic algorithms (Wang, 1997)
generally have the drawback of being inefficient andrequiring too many iterations. As an example, Gupta
et al. (2003) reported that a complete single-criterion
optimization problem using the shuffle complex model
as reported in Sorooshian et al. (1993) would requireas many as 10 000 function calls, while 100 Pareto solutions will require approximately one million function
calls. In a multiobjective formulation, a Pareto solution
is a solution that cannot be improved without disadvantaging at least one of the objectives. To overcome someof these problems we developed a procedure that consists of a sequence of steps in which the initial (large)
uncertainty in the model parameters is progressively
reduced until certain calibration criteria for prediction
uncertainty are met. The program uses an efficient sampling procedure (Latin hypercube), along with a global
search algorithm that examines the behavior of an objective function by analyzing the Jacobian and Hessian
matrices. The current algorithm is henceforth referred
to as SUFI-2, as it is the second version of a previously
developed Sequential Uncertainty FItting (SUFI) program (Abbaspour et al., 1997).
In this study, the SUFI-2 procedure was applied to
two municipal solid waste incinerator (MSWI) bottom
ash monofills. Landfills of this type have recently been
investigated by, among others, Hjelmar (1996), Belevi
et al. (1992), Johnson et al. (1998), and Hartmann et al.
(2001). Bottom ash in these landfills is typically composed of equal amounts of fine ash material and melted
components (of which one-half have crystallized) andsmall quantities of metallic components, ceramics, andstones (Kirby and Rimstidt, 1993). The Lostorf landfill,a MSWI bottom ash monofill in Switzerland, was stud-

ied in detail by Johnson et al. (1999). Chemical analysesof leachate from this landfill at discrete time intervals
between 1994 and 1996 showed average total concentrations of 44.5, 47.1, 11.8, 0.63, 8.2, and 12.4 mM for Na,
Cl, K, Mg, Ca, and SO4, respectively. Many other metals,
such as Cu, Zn, Sb, Cr, Cd, Mo, V, Mn, and Pb, werealso detected. While the leachate composition was foundto be relatively constant during dry periods, considerable dilutions occurred during rain events. The relatively
good reproducibility of the experimental observations
in response to rain events motivated us to use transportmodeling to predict heavy metal concentrations.
Different models have been used to simulate flow
through MSWI bottom ash landfills (Guyonnet et al., 1998;
Hartmann et al., 2001; Johnson et al., 2001). In the studyof Johnson et al. (2001), flow through the Lostorf landfill
was modeled using the programs MACRO (Jarvis, 1994),
HYDRUS5 (Vogel et al., 1996), a neural network approach (Schaap and Bouten, 1996), and a linear storage
model adapted from Huwe et al. (1994). They foundthat flow was dominated by preferential paths. The best
simulation results were obtained with the variably saturated dual-permeability model MACRO of Jarvis (1994).
We extend the work of Johnson et al. (2001) by simulating both flow and transport at the Lostorf landfill site.
Solute transport was simulated in terms of the electricalconductivity (EC). A second bottom ash landfill, Seckenberg, was also calibrated for flow. This was done bylinking SUFI-2 with MACRO. Our objectives were to
test the SUFI-2 procedure by calibrating and testing the
MACRO model using hourly discharge and EC data.
THEORY
Simulation Program MACRO
MACRO is a dual-permeability model based on the conceptof having two different domains: a micropore or matrix domain, in which flow is governed by the Richards equation,
and a macropore domain that drains by gravity flow. The
Richards equation for one-dimensional isothermal flow in avariably saturated rigid porous medium is given by

冤 冢

冣冥

⳵␪
⳵
⳵h
⫽
⫺ 1 ⫾ Sw
K
⳵t
⳵z
⳵z

[1]

where ␪ is the volumetric water content (L3 L⫺3), K is thehydraulic conductivity (L T⫺1), h is the pressure head (L), Sw
is a term signifying a source or a sink (T⫺1), and z is depth(L) positive downward. For the macropore region, a simplified
approach is used in which fluxes are predicted by assuming
no change in h with depth (i.e., ⳵h/⳵z ⫽ 0), thus reflecting
laminar flow under gravity.
The source–sink term Sw in Eq. [1] is used here to accountfor water exchange between the micropores and macropores.
The flux between the macropores to the micropores is given bySw ⫽冢3Dd ␥ 冣(␪ ⫺ ␪ )
w w
2

b

mi

[2]

where d is the effective diffusion path length (L), ␪b is thesaturated water content (L3 L⫺3), ␪mi is the actual water content
in the micropores (L3 L⫺3), Dw is the effective water diffusivity(L2 T⫺1), and ␥w is a scaling factor introduced to match the approximate and exact solutions of the diffusion problem. The

1342

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

冢

冣

⳵c
⳵␪c
⳵
D␪ ⫺ qc ⫾ U
⫽
⳵t
⳵z
⳵z

[8]

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

where c is the solute concentration of the liquid phase (ML⫺3), q is the Darcian water flux (L T⫺1), U is a source–sinkterm (M L⫺3 T⫺1), and D is the dispersion coefficient (L2 T⫺1)given byD ⫽ Dvv ⫹ D0 f*[9]2T⫺1), f*where D0 is the diffusion coefficient in free water (Lis a constant impedance factor, Dv is the dispersivity (L), andv is the pore water velocity (L T⫺1). Dispersion in the macropores is set to zero since solute transport is assumed to be
dominated by advection. The source–sink term U accounts
for diffusive and advective mass transfer between the two flow
domains as follows:
Fig. 1. Water retention and hydraulic conductivity curves as usedin MACRO.

scaling factor varies somewhat with the initial water content
and hydraulic properties, but not strongly (Gerke and van
Genuchten, 1993). This parameter is set to an average value of
0.8 in MACRO. The model assumes that water is immediatelyrouted from the micropores to the macropores when watercontent exceeds saturation, ␪b. Dispersion in the macroporesis neglected since advection is assumed to dominate solutetransport.
Different functions are used in MACRO to describe thewater retention characteristics of the micropore and macropore regions, as illustrated in Fig. 1. For the micropores,
the water retention function, hmi(␪), is given by the modified
Brooks and Corey (1964) function:
1/␭hmi(␪) ⫽ hbS⫺mi

␪r ⱕ ␪ ⱕ ␪ b[3]

␪ ⫺ ␪r␪b ⫺ ␪r

␪r ⱕ ␪ ⱕ ␪ b

[4]

where

Smi ⫽

and hb is the critical boundary head (L) at which microporesdrain (also known as the air-entry value), ␪r is the residual
water content, and ␭ is the pore size distribution index. The
hydraulic conductivity function for the micropore region, Kmi,
is derived from the retention curve using Mualem’s (1976)
model, yieldingn⫹2⫹2/␥
Kmi ⫽ Ks,miS mi[5]

where Ks,mi is the hydraulic conductivity (L T⫺1) at the criticalboundary head hb, and n is the tortuosity factor of the micropore region.
Because flow in the macropore domain is assumed to be
driven by gravity only, no retention curve is required. The hydraulic conductivity of the macropores, Kma, is given by asimple power function of the degree of saturation in the macropores (Sma):n
Kma ⫽ Ks,maS ma*[6]

where ␪s is the saturated water content, n* is an empiricalparameter (unit-less) describing the tortuosity in macropores,
and Sma is defined as

Sma ⫽

␪ma
␪s ⫺ ␪ b[7]

where ␪ma is the water content in macropores.
Solute transport in MACRO is described with the standard
advection–dispersion equation for nonsorbing tracers:U⫽冢3Dd ␪ 冣(c ⫺ c ) ⫹ S c⬘
e mi
2

ma

mi

w

[10]

where the prime indicates solute concentration in either themacropores or micropores depending upon the direction offlow (the sign of Sw), and De is an effective diffusion coefficient(L2 T⫺1) given byDe ⫽ D0 f*Sma[11]

Infiltration from the soil surface into the micropores is assumed to be limited by the saturated conductivity, Ks,mi, of the
micropore region. The portion of rain or irrigation exceeding
Ks,mi flows directly into the macropores. We calculated thesolute concentration of water infiltrating into the macropores,
c*ma , by assuming complete mixing of the infiltrating water withwater stored in a shallow “mixing depth,” zd:

c*ma ⫽

Qd(t⫹⌬t) ⫹ Rcr
R ⫹ zd␪1[12]where Qd is the amount of solute stored in the mixing depth(M L⫺2), R is the amount of rain reaching the surface withintime interval ⌬t (L), cr is the concentration in the rain water(M L⫺3), and ␪1 is the water content in the layer zd. The amountof solute added or removed from the micropores in the top
layer within time interval ⌬t, Qmi(l), is calculated as the difference:

Qmi(l) ⫽ Rcr ⫺ Imac*ma[13]

where the Ima denotes infiltration in macropores (L).

Inverse Procedure SUFI-2
We assume that inversely obtained parameters are always
uncertain and that an inverse optimization approach shouldseek the smallest possible parameter uncertainty or parameter
range that satisfies a desired prediction uncertainty. The SUFI-2
procedure described below uses elements of Generalized Likelihood Uncertainty Estimation (GLUE) (Beven, 1989, Bevenand Binley, 1992) with elements of a gradient approach (Kooland Parker, 1988), modified to enable a global search. Similar
to GLUE, the procedure is applied to parameter sets, rather
than to individual parameter values, so that any interactions
between parameters are taken into account explicitly.
Also similar to GLUE, our final objective is not necessarily
to find a set of best-fit parameters, even though such a bestfit parameter set is calculated. Unlike GLUE, however, thekey output of SUFI-2 is a “best range” for each parameter.
Parameter combinations within the parameter ranges are ensured to produce high quality simulations because of the two
calibration criteria discussed below.

1343

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 2. An example showing the partitioning of a response signal into different compartments: base flow (c1), the recession or rise (c2), intermediatepeaks (c3), and large peaks (c4). The weights are calculated such that each compartment contributes equally to the objective function.

Starting with some initial parameter value ranges, SUFI-2
is iterated until (i) the 95% prediction uncertainty (95PPU) between the 2.5th and 97.5th percentiles bracket more than 90%of the measured data and (ii) the average distance between
the 2.5th and 97.5th percentiles is smaller than the standard
deviation of the measured data. A model is considered calibrated if, upon reaching the above two criteria, a significant
R2 exists between the best simulation and calibration data.
Further testing of the calibrated parameter ranges, as performed in this study, should increase model confidence. If
calibration cannot be attained with the above criteria, themodel structure should be reexamined.

could be given a larger weight. If two variables such as discharge and concentration are considered simultaneously, thena multiobjective, multicompartment formulation with eight
RMSE values is defined, in which case the objective function becomes

Step-by-Step Application of SUFI-2

where the average (avg) is computed for all simulations withina Latin hypercube sampling round (see Step 5), w1 is always
equal to 1, andStep 1In this first step an objective function is defined. The literature shows many different ways of formulating an objectivefunction (e.g., Legates and McCabe, 1999; Gupta et al., 1998).
Each formulation may lead to a different result; hence, the final parameters and parameter ranges are always conditionedon the form of the objective function. To overcome this problem, some studies (e.g., Yapo et al., 1998) combine different
types of statistics together (e.g., RMSE, absolute difference,logarithm of differences, R2, ␹2) to yield a “multicriteria” formulation. The use of multiobjective formulation (Duan et al.,
2003; Gupta et al., 1998), where different variables, including different variables at different locations, are embedded in
the objective function, is also important to reduce the nonuniqueness.
Since there is no unique formulation of the objective function that can be applied universally to all situations, the choiceof an objective function must correspond to the goal of the
project. In this study, we required a good overall simulation of
the peaks as well as the base flow; hence, we used the formulation of Boyle et al. (2000) and partitioned the response
function into four compartments, as illustrated in Fig. 2. This
formulation, henceforth referred to as “multicompartment”
formulation, depicts base flow (c1 in Fig. 2), rise and recession(c2), intermediate peaks (c3), and large peaks (c4). For eachcompartment, the RMSE is calculated and weighted so they
have equal contributions to the objective function. This formulation has proven to provide an overall good fit if the model
has no conceptual limitation to simulating each compartmentof the response function (i.e., not being able to simulate large
peaks due to the absence of a macropore flow component).
If simulation of the large peaks is important, the peak regionIg(b) ⫽ 兺 wi RMSEi

[14]

i⫽1

where b is the parameter vector, I (here equal to 8) is thetotal number of compartments, and wi is calculated here aswi ⫽avg(RMSE)1
avg(RMSE)i

RMSEi ⫽1 ki ⫽ 1, …, I

冪k 兺 (q ⫺ q )
j⫽1

o

s 2
j

[15]

[16]

where q is a measured variable, k is the number of observationsin the ith compartment, and the superscripts “o” and “s” refer
to observed and simulated values, respectively. Equation [15]implies that the weights are different for each sampling round.
Step 2The second step establishes physically meaningful absolute
minimum and maximum ranges for the optimized parameters.
There is no theoretical basis for excluding one specific distribution. However, because of lack of information, we assume thatall parameters were uniformly distributed within the region
bounded by minimum and maximum values. Because the parameter ranges play a constraining role, they should be aslarge as possible, yet physically meaningful:

bj : bj,abs_min ⱕ bj ⱕ bj,abs_max

j ⫽ 1.... m[17]

where bj is the jth parameter and m is the number of parameters to be estimated.
Step 3This step involves an optional, yet highly recommended
sensitivity analysis for all parameters. We pose that no automated optimization routine can replace the insight from physical understanding and knowledge of the effects of parameters
on system response. The sensitivity analysis is carried out by

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1344

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

keeping all parameters constant to realistic values, while varying each parameter within the range assigned in Step 1. For
each parameter, about 10 simulations are conducted by simply
dividing the absolute ranges in equal intervals and allowing the
midpoint of each interval to represent that interval. Plotting
results of these simulations along with the observations on
the same graph gives insight into the effects of the parameters
on observed signals.

where b*j is the current best estimate of parameter b (i.e., thatparameter value that produces the smallest value of the objective function), and ␯ is the degrees of freedom (n ⫺ m).
Parameter correlations can be assessed using the diagonal andoff-diagonal terms of the covariance matrix as follows:

Step 4

It is important to note that the correlation matrix A quantifies the change in the objective function as a result of a changeof parameter i, relative to changes of the other parameters j.
Parameter sensitivities, S, are calculated next by averaging
the columns of the Jacobian matrix as expressed by

Initial uncertainty ranges are assigned to parameters forthe first round of Latin hypercube sampling; that is,

bj : [bj,min ⱕ bj ⱕ bj,max]j ⫽ 1, m[18]

In general, the above ranges are smaller than the absolute
ranges, are subjective, and depend on experience. The sensitivity analysis in Step 3 can provide a valuable guide for selecting appropriate ranges. Although important, these initial estimates are not crucial since they are updated and allowed tochange within the absolute ranges.
Step 5A Latin hypercube (McKay et al., 1979) sampling resultsin n parameter combinations, where n is the number of desired
simulations. This number should be relatively large (≈1000–
2000). The simulation program is run n times and the simulated
output variable(s) of interest, corresponding to the measurements, are saved.

Aij ⫽

Cn

Cij
√Cii√Cjj兩 兩

1 2 ⌬gi
S j ⫽ bj n 兺
C 2 i⫽1 ⌬bj

j ⫽ 1,…, m

[25]

[26]

where bj is the average value of the jth parameter. We emphasize that the measures of sensitivity given by Eq. [26] aredifferent from the sensitivities calculated in Step 3. The sensitivities given by Eq. [26] are estimates of the average changes
in the objective function resulting from changes in each parameter, while all other parameters are changing. Therefore, thesensitivities in Eq. [26] are relative sensitivities. By contrast,
the sensitivities in Step 3 define the absolute sensitivity of aparameter that can change when other parameters take on
different optimized values.
Step 8

Step 6
As a first step in assessing the simulations, RMSEs arecalculated followed by calculation of the weights and the objective function for each simulation.
Step 7In this step a series of measures is calculated to evaluateeach sampling round. First, the sensitivity matrix or Jacobianof g(b) is computed using

Jij ⫽

⌬gi
i ⫽ 1,…, C 2n, j ⫽ 1,…, m
⌬bj

[19]

where C2n is the number of rows in the Jacobian (equal to all
possible combinations of two simulations), and j is the numberof columns (number of parameters). Following the Gauss–
Newton method and neglecting the higher-order derivativesof the Hessian matrix, H, the Hessian matrix of g(b) is calculated as

H ⫽ JTJ[20]

According to the Cramer–Rao theorem (Press et al., 1992),
an estimate of the lower bound of the parameter covariancematrix, C, is calculated fromC ⫽ s2g(JTJ)⫺1

[21]

2
g

where s is the variance of the objective function. The estimated standard deviation and 95% confidence interval of aparameter bj is calculated from the diagonal elements of C
(Press et al., 1992) fromsj ⫽ √Cjj[22]

bj,lower ⫽ b*j ⫺ t␯,0.025 sj[23]bj,upper ⫽ b*j ⫹ t␯,0.025 sj[24]In this step the 95PPUs of all predicted variables are computed using the n Latin hypercube simulations. Distributionof the predicted variables is not always Gaussian, and could
be highly skewed. Thus, the usual calculation of uncertainty
limits as a function of the variance of the predicted values is
not applicable (Beven and Binley, 1992). In this study we
represent the prediction uncertainty as the 2.5th (qL) and
97.5th (qU) percentiles of the cumulative distribution of every
simulated point.
The goodness of fit is calculated for each variable q fromthe percentage of measured data that fall within the 95PPU
region (prediction uncertainty), the R2 between the optimized
and observed data, and the average distances d between theupper and lower 95PPU as determined fromd⫽1 k
兺 (qU ⫺ qL)l
K l⫽1

[27]

in which l is a counter, and K is the total number of observations for variable q. The best outcome is that 100% of themeasurements fall within the 95PPU, R2 is close to 1, and thed value is close to zero. However, because of measurement
errors and model uncertainties, generally this will not be
the case.
We consider a model satisfactorily calibrated if approximately 90% of the measured data are within the 95PPU whend is smaller than the standard deviation of the measured data
and when R2 ⬎ 0.8.
Step 9Because in general parameter uncertainties are large initially, the value of d tends to be quite large in the first samplinground when about 100% of the data are within the 95PPU.
Hence, further sampling rounds are needed with updated parameter ranges calculated from:1345

www.vadosezonejournal.org

Table 1. Data periods for calibration and testing.
Landfill

Calibration†

Testing(yr-mo-d-h)

Lostorf
Seckenberg

1996-04-22-13 to 1996-11-10-20
1998-04-01-00 to 1998-12-03-23

1995-01-02-00 to 1995-03-31-13
1999-01-01-00 to 1999-05-26-20

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

† The year with larger amount of data was chosen for calibration.

b⬘j,min ⫽ bj,lower ⫺ Max

冤(bj,lower⫺ bj,min) (bj,max ⫺ bj,upper)
,
2
2

冤

冥

(b
⫺ bj,min) (bj,max ⫺ bj,upper),b⬘j,max ⫽ bj,upper ⫹ Max j,lower
22

冥 [28]

where b⬘ indicates updated values. The above criteria, while
producing smaller parameter ranges, ensure that the updated
parameter ranges are always centered on the current best
estimates. If the best estimates are close to their limits, parameter ranges are increased while not exceeding the absoluteboundaries. In the final step, parameters are ranked according
to their sensitivities, and highly correlated parameters are
identified. Of the highly correlated parameters, those with thesmaller sensitivities should be fixed to their best estimates andremoved from additional sampling rounds.
This procedure is continued until the two stopping rules
for model prediction uncertainty are satisfied. The resulting
parameter ranges are considered the best-fit parameters forthe problem being investigated. If the stopping rules cannot
be met (i.e., a small value of d, but too many observed data
are excluded from the 95PPU), the model cannot be calibrated
and the structure of the model and its boundary conditions
need to be reexamined.

APPLICATION
Sites and Data
The SUFI-2 sequential uncertainty fitting procedure wasapplied to the two MSWI landfills in Switzerland. The Lostorf
landfill, a MSWI bottom ash monofill near Buchs AG, Switzerland, is situated in an unused gravel pit. The Seckenberg landfill is situated on sandstone and dolomite formations on ahillside above Frick, Switzerland. For the Lostorf landfill, measured discharge and EC data were available for several months
in 1995 and 1996. For the Seckenberg landfill only flow data
were available for a few months in 1997 and 1998. We refer
the reader to Johnson et al. (1999, 2001) for more complete
descriptions of the two landfills.
Hourly values of rainfall, global radiation, relative humidity,
air pressure, wind velocity at 2 m, and temperature were obtained from the Swiss Meteorological Institute at the Buchs/
Suhr and Frick/Rheinfelden sites for application to Lostorf
and Seckenberg, respectively. Potential evapotranspiration
rates from landfill surfaces were estimated using the hourly
form of the FAO Penman–Monteith equation (Allen et al.,1994), but modified to obtain actual evaporation rates for bothlandfill systems according to procedures described by Johnson
et al. (2001).

Table 2. Estimated flow and transport parameters for the Lostorf and Seckenberg landfills
Parameter†

Pseudovalue

Initialrange

Optimized
range

Optimized
value

Sensitivity
ranking

ci–1, mS mm⫺1
ci–2
ci–3
␪i–1, cm3 cm⫺3
␪i–2
␪i–3
cr, mS mm⫺1
Dv, mmd-1, mmd-2
d-3
hb–1, cmhb–2
hb–3
Ks–1, mm h⫺1
Ks–2
Ks–3
Ks,mi–1
Ks,mi–2
Ks,mi–3
␪r–1, cm3 cm⫺3
␪r–2
␪r–3
␪s–1
␪s–2
␪s–3
␪b–1, cm3 cm⫺3
␪b–2
␪b–3
␭-1
␭-2
␭-3
n-1
n-2
n-3
n*-1
n*-2
n*-3

0.767
0.974
1.39
14.23
12.14
16.56
0.5
97.4
111.55
126.52
3.70
18.92
7.13
15.28
276.20476.43142.31
1.075
0.017
0.132
14.23
12.14
16.56
30.34
27.92
27.66
25.84
22.27
24.65
0.183
0.174
0.426
0.071
0.152
0.097
3.691
2.767
0.682

0.1–2
0.5–1.5
1.4–1.5
17–23
17–23
19–25
0.2–1
10–100
20–100
20–100
20–100
5–25
5–25
5–25
150–500
150–500
150–500
0–2
0–1
0–1
12–16.5
12–16.5
12–18.5
27–30
27–30
26–30
23–26.9
23–26.9
23–25.9
0.1–0.4
0.1–0.4
0.1–0.4
0.01–0.15
0.01–0.15
0.01–0.15
1–4
1–4
1–4

0.57–0.67
0.83–1.09
1.38–1.40
14.7–15.3
13.1–15.4
16.4–17.6
0.47–0.54
95.4–102.3
106.4–114.9
124.7–134.2
2.7–3.9
18.1–19.2
6.0–7.0
14.4–15.5
266.7–300.2
432.3–468.9
138.8–147.6
0.99–1.16
0.014–0.02
0.12–0.154
12.0–14.0
12.1–12.6
16.7–17.8
29.9–31.3
27.5–29.9
27.5–28.3
25.5–26.3
22.1–23.4
24.1–24.9
0.17–0.20
0.17–0.20
0.47–0.51
0–0.040.15–0.164
0.084–0.095
3.5–3.9
2.7–2.84
0.62–0.69

0.634
1.06
1.39
14.86
14.4
16.63
0.5
95.7
113.5
129.6
3.3
18.69
6.31
15.42
280.44
468.7
140.4
1.019
0.019
0.151
13.2
12.2
17.2
30.2
27.7
27.9
25.62
22.90
24.69
0.181
0.192
0.492
0.024
0.160
0.093
3.77
2.83
0.67

31
35
1
3
16
8
29
14
18
17
37
11
30
19
25
20
12
33
36
34
32
6
13
7
21
2
4
10
5
26
28
15
38
24
27
22
9
23

† Numbers 1, 2, and 3 refer to soil layers.

1346

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Table 3. Estimated flow and transport parameters for the Lostorf
and Seckenberg landfills.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Optimized calibration range
Parameters†

Lostorf

Seckenberg

ci–1, mS mm⫺1
ci–2
ci–3
␪i–1, cm3 cm⫺3
␪i–2
␪i–3
cr, mS mm⫺1
Dv, mmd-1, mmd-2
d-3
hb–1, cmhb–2
hb–3
Ks–1, mm h⫺1
Ks–2
Ks–3
Ks,mi–1
Ks,mi–2
Ks,mi–3, mm h⫺1
␪r–1, cm3 cm⫺3
␪r–2
␪r–3
␪s–1
␪s–2
␪s–3
␪b–1
␪b–2
␪b–3
␭-1
␭-2
␭-3
n-1
n-2
n-3
n*-1
n*-2
n*-3

0.6–1.2
0.6–1.0
1.3–1.4
0.18–0.20
0.16–0.19
0.20–0.22
0.22–0.52
57–108
84–112
94–1372–27
14.8–20.3
2.9–8.1
13.6–19.1
248–409
400–510
135–265
0.88–1.47
0.0–0.16
0.08–0.26
0.13–0.15
0.12–0.13
0.15–0.18
0.29–0.31
0.28–0.29
0.27–0.28
0.26–0.27
0.22–0.24
0.24–0.25
0.15–0.25
0.10–0.26
0.34–0.45
0.04–0.07
0.15–0.22
0.08–0.17
2.4–3.7
2.6–3.2
0.46–1.3–––
0.24–0.26
0.22–0.25
0.22–0.24
–
–
66–74
38–50
33–49
19.9–21.8
10.6–15.3
10.1–13.8
318–352
183–238
421–494
1.1–1.3
0.0–0.003
0.11–0.26
0.13–0.15
0.18–0.19
0.16–0.18
0.28–0.29
0.28–0.30
0.27–0.28
0.25–0.26
0.23–0.25
0.23–0.24
0.11–0.13
0.29–0.39
0.30–0.36
0.07–0.08
0.037–0.59
0.08–0.09
3.1–3.5
2.4–2.6
1.1–1.25† Numbers 1, 2, and 3 refer to soil layers.

Calibration and Testing of the Flowand Transport Model
The Lostorf landfill was calibrated using both dischargeand EC in the objective function given by Eq. [14] and [15].
The objective function for the Seckenberg landfill included
only discharge data. Time periods used for calibration andmodel testing are given in Table 1.
Based on the filling history of the landfills, a three-layer
profile was modeled by treating the parameters listed in Table
2 as unknowns. A common problem of historic landfill databases is that the initial conditions for both the flow and transport simulations are poorly defined. To determine the initialwater content and electrical conductivity profiles, we examined two options: (i) simulating the flow regime by letting
the landfills drain from a nearly saturated condition until thesimulated discharge equaled the observed discharge at time
zero (the water content and concentration profiles at this time
were assumed the correct initial water content and solute
distributions) and (ii) fitting the initial conditions along withthe other parameters. Although both options have problems,
we opted for the second option since it was more straightforward and produced better results. The main problem with thefirst option was that flow parameters defining flow rates werestill unknown. To reduce the effect of the initial condition,
we started the simulations a few months before time zero.
Including the initial water content and concentration distributions produced a total of 33 unknown parameters when only
flow is considered, and 38 parameters when transport is modeled (Table 2). Transport simulations were performed by using
electrical conductivity as a tracer.
The bottom boundary conditions were assumed to haveconstant hydraulic gradients, while the surface boundaries
were assumed to be atmospheric. Rain or irrigation exceeding
Ks,mi was allowed to flow directly into the macropores.

RESULTS AND DISCUSSION
The absolute parameter ranges and the starting initial
ranges were selected based on our experience with theFig. 3. Calibration results for flow and transport at Lostorf. The blue curves represent the 95% prediction uncertainty (95PPU), the red curve
represents the measured data, and the green curve shows the best fit of optimization.

1347

www.vadosezonejournal.org

Table 4. Statistics of measured variables and the corresponding modeling results.
Variable
Lostorf discharge calib. test
Lostorf EC calib. test

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Seckenberg discharge calib. test

Number ofdata points

Average and SD ofmeasured variables

Percent of measured
data in the 95PPU†R2 ‡Average width of 95PPU and its
ratio over SD of measured data

4855
1600
1600
2125
5927
3500

0.062 (0.077)
0.15 (0.14)
1.15 (0.21)
1.02 (0.24)
0.105 (0.16)
0.19 (0.33)

90
91
90
90
93
89

0.86
0.85
0.88
0.82
0.94
0.88

0.06 (0.78)
0.16 (1.1)
0.23 (1.1)
0.40 (1.7)
0.07 (0.44)
0.12 (0.36)† 95PPU, 95% prediction uncertainty.
‡ All regressions were significant at P ⬍ 0.0001.

parameters, literature review (e.g., Dubus and Brown,2002), and personal communication with the author of
MACRO (Jarvis, 1994). We conducted the sensitivity
analysis of Step 3 after the first sampling round. We
concluded that (i) some parameters affected the earlystages of the discharge and concentration (e.g., initialwater contents and concentrations), while others affected the later stages (e.g., dispersivity, Dv; effective
diffusion path length, d), (ii) some parameters had astrong effect on the peak discharge (e.g., hydraulic conductivities Ks,mi and Ks,ma), while others affected mainlybase flow (especially the diffusion path length, d), and(iii) some parameters affected the recession or timingof the peak discharges (e.g., the pore size distribution
index, ␭, and the macropore tortuosity index n*), whileothers affected all compartments as defined in Eq. [15](especially the saturated hydraulic conductivity, Ks,mi).
We used this sensitivity analysis to obtain better estimates of the initial parameter ranges and to improve
our understanding of the effect of all parameters ondischarge and EC.
On the basis of some preliminary runs we selected a
set of pseudo parameters for the Lostorf landfill andgenerated a synthetic set of discharge and EC data. We

used these output data in SUFI-2 to examine how well
we could recover the pseudo parameters. Table 2 summarizes the list of pseudo parameters, the initial parameter uncertainties, the optimized parameter uncertainties,
the optimized best parameter values, and the relative
sensitivity ranking. In the fourth Latin hypercube sampling round, 99% of the measured EC and discharge
data were bracketed by the 95PPU, while the ratio of
average difference between the upper and the lower
95PPU over the standard deviation of generated measured data was 0.83 for EC and 0.53 for discharge, indicating excellent fits. As shown in Table 2, however, someof the insensitive pseudo parameters (i.e., ci–1) wereoutside the optimized range, indicating that parametersinsensitive to the objective function could not be fitted.
Table 3 summarizes the final optimized parameter
uncertainties for the Lostorf landfill, while Fig. 3 presents calibration results for flow and transport at Lostorf. Plots in Fig. 3 show the best simulation results andthe upper and lower limits of the 95PPU of the calibrated parameter ranges. Results were obtained after
five sampling rounds for a total of 5000 simulations.
The 95PPU contained 90% of the measured dischargeand EC data (Table 4), while the ratio of the averageFig. 4. Test results of flow and transport at Lostorf. The blue curves represent the 95% prediction uncertainty (95PPU), the red curve represents
the measured data, and the green curve shows the best fit of optimization.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1348

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Fig. 5. (A) Calibration and (B) test results for flow at Seckenberg. The blue curves represent the 95% prediction uncertainty (95PPU), the red
curve represents the measured data, and the green curve shows the best fit of optimization.

difference between the upper and the lower 95PPU overthe standard deviation of the measured data was ⬍1 fordischarge and only slightly above 1 for EC. With R2values for discharge and EC of 0.86 and 0.88, respectively, almost all calibration requirements for the calibration data set of Lostorf were satisfied.
We tested the above calibration parameters using
data from a previous year; results are shown in Fig. 4.
The statistics in Table 4 indicate that the three calibra-tion requirements were again satisfied for discharge, butthat the ratio requirement was not met for the EC data.
As will be discussed below, this indicates that the calibrated parameter ranges produced several relatively
poor simulations for EC. This is a strong indication thata higher degree of multiobjective formulation is neededin the objective function to limit the number of bad simulations.
To test the SUFI-2 procedure against previous results

Fig. 6. The d-3–hb-3 response surface obtained after 2000 simulations. All other parameters were kept constant.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

1349

Fig. 7. The d-3–␪b-3 response surface obtained after 2000 simulations. All other parameters were kept constant.

by Johnson et al. (2001) for this same landfill, we compared the sum of square errors (SSQ) between the measured and simulated results for both calibration and test
cases. A direct comparison is not possible because theobjective function previously did not include a transportcomponent. The previously obtained SSQs for discharge
calibration and test results were 10.1 and 41.1, respectively. The values obtained with SUFI-2 were 4.1 and
7.5, respectively. This shows a significant improvement
in best-fit simulations obtained with our new sequential
uncertainty fitting procedure.
Figure 5 presents the calibration (A) and test (B)results for flow through the Seckenberg landfill. The
plots show both the best simulation results and the 95PPU
based on the calibrated parameter ranges of Table 3.
The results were obtained after four sampling rounds.
The 95PPU contained 93 and 89% of the measured discharge data for calibration and testing, respectively (Table 4). The ratios of the average difference between theupper and the lower 95PPU over the standard deviationof the measured data were substantially less than one forboth the calibration and test data seta. The R2 requirement was also met with highly significant values. One
reason for the better calibration results for the Secken-

Fig. 8. The d-3–n* response surface obtained after 2000 simulations. All other parameters were kept constant.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1350

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Fig. 9. The d-3–n* response surface obtained after 10 000 simulations. All other parameters were kept constant.

berg case is that the objective function contained only
discharge data.
A comparison of the parameters for both landfills(Table 3) shows substantial overlaps, with the parameters for Seckenberg generally having smaller ranges.
Applying the Lostorf parameter ranges to Seckenberg(except for the initial values which were those of theSeckenberg calibration) and vice versa produced relatively good simulation results, even though not all threecalibration criteria were met. This suggests that it may
be possible to use identical parameters for similar applications.

To investigate the progression toward convergenceand limitations of the calibrated parameter ranges, weplotted the response surfaces for some of the most sensitive parameters. In doing so, all parameters of the Seckenberg example were fixed at their final best values.
In the first example diffusion path length, d, and theboundary soil water tension hb of the third layer werevaried by drawing 2000 samples using the Latin hypercube procedure. The response surface of the objectivefunction is plotted in Fig. 6. Superimposed on the response surface are the parameter ranges from the last
three sampling rounds, delineated by three boxes num-Fig. 10. The d-3–␪b response surface obtained after 10 000 simulations. Parameters d-3, ␪b-3, hb-3, and n* were varied simultaneously, while allothers were kept constant.

1351

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

bered 1, 2, and 3. The two-parameter d–hb response
surface shows a clear valley of small objective function values.
The results in Fig. 6 reflect an apparent limitation of
expressing the parameters in terms of ranges. The twoparameter surface plot suggests that the uncertainty
ranges of both parameters could be much larger, albeit
located along a narrow valley. However, if the size ofthe parameter ranges is increased, the prediction uncertainty may include many poor simulations as was thecase with the EC of Lostorf, leading to a high value ofthe second criterion. For this reason, several (e.g., Yapoet al., 1998) have used the concept of Pareto optimality,in which the simulations are divided into sets of relatively good and poor simulations. The final solution
is then not specified by means of a set of parameter
distributions, but in terms of parameter combinations
that produce good simulations. Still, this approach also
does not guarantee that all good simulations can be
captured by means of a limited number of function calls.
Additional calculations showed that the two-parameter response surface may be somewhat misleading whenassessing the region with minimum values in a multidimensional parameter system. Figure 7 shows the response surface when only d and ␪b are allowed to vary.
Cleary, the region with minimum values of the objectivefunction now has a different shape, where many disconnected minima islands exist. Another example with d–n*(Fig. 8) shows an even more complicated surface, withnumerous local minima.
To verify that 2000 samples will give an accurate response surface, we repeated the d–n* simulations with
10 000 samples. Figure 9 shows a plot similar to that ofFig. 8 for the 2000 simulations. Notice that the response
surface again contains many minima. If more than two
parameters are allowed to vary simultaneously, plots ofany combination of two-parameter surfaces will be evenmore complicated than those in Fig. 8 or 9. Figure 10
shows plots of the surface of the d–␪b response function
when d, n*, ␪b, and hb were allowed to vary simultaneously. This graph again illustrates the difficulty of using
procedures based on gradient methods to solve optimization problems.

CONCLUSIONSThe SUFI-2 procedure we developed for inverse modeling uses a sequence of steps in which the initial (large)
uncertainties in the model parameters are progressively
reduced until a certain calibration requirement basedon the prediction uncertainty is reached. From a seriesof two-parameter response surfaces we observed that
also the SUFI-2 program is limited to finding one ofmany regions with local minima in a multidimensionalparameter space. An encouraging observation, however,
was that none of the detailed runs produced a region
with smaller values of the objective function. The advantages of SUFI-2 are that it combines optimization withuncertainty analysis and can handle a large number ofparameters. Moreover, the procedure is quite efficientin terms of the number of function calls, and it is easy
to implement.
Finally, we recall that one of the main purposes of
inverse modeling is to use easily measurable model outputs to obtain parameters that are much more difficult
and time-consuming to measure directly. In other words,
we prefer to feed our inverse optimization procedure
with data of relatively low value to obtain data of higher
value. Unfortunately, the notion of getting much forlittle is inherently problematic. Despite the attractions,nature does not allow for such easy transformations.
For this reason one must accept that our SUFI-2 procedure, or similar Bayesian inverse methods, will only
produce a certain solution space, rather than a single
unique solution. Hence, one may very well have to livewith the idea that inverse methods will lead to manysolutions regardless of the type of the objective functionused in the optimization. To limit the nonuniqueness,
however, a step forward would be to make the description of the objective function as constraining to parameters as possible by using a combination of multicriteria,
multiobjective, and multicompartment formulation.
ACKNOWLEDGMENTS
The authors thank the Gemeindeverband Abfallbeseitigung
Oberes Fricktal, and in particular Erwin Roth, for their active
cooperation, Gholamabbas Sayyad for data processing andinitial phases of modeling, Dr. Mark Borsuk for reviewing apreliminary version of the paper, and Dr. Peter Reichert for
many fruitful discussions.

REFERENCES
Abbaspour, K.C., M.Th. van Genuchten, R. Schulin, and E. Schläppi.
1997. A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters. Water Resour.
Res. 33:1879–1892.
Abbaspour, K.C., R. Schulin, and M.T. van Genuchten. 2001. Estimating unsaturated soil hydraulic parameters using ant colony optimization. Adv. Water Resour. 24:827–841.
Allen, R.G., M. Smith, L.S. Pereira, and A. Perrier. 1994. An update
for the calculation of reference evapotranspiration. ICID Bull. 43:
35–92.
Beck, J.V., and K.J. Arnold. 1977. Parameter estimation in engineeringand science. John Wiley and Sons, New York.
Belevi, H., D.M. Stämpfli, and P. Baccini. 1992. Chemical behaviourof municipal solid waste incinerator bottom ash in monofills. Waste
Manage. Res. 10:153–167.
Bell, L.S.J., P.J. Binning, G. Kuczera, and P.M.H. Kau. 2002. Rigorous
uncertainty assessment in contaminant transport inverse modeling:
A case study of fluoride diffusion through clay liners. J. Contam.
Hydrol. 57:1–20.
Beven, K.J. 1989. Interflow. p. 191–219. In H.J. Morel-Seytoux (ed.)
Unsaturated flow in hydrologic modeling. D. Reidel, Norwell, MA.
Beven, K., and A. Binley. 1992. The future of distributed models:Model calibration and uncertainty prediction. Hydrol. Processes
6:279–298.
Boyle, D.P., H.V. Gupta, and S. Sorooshian. 2000. Toward improved
calibration of hydrological models: Combining the strengths ofmanual and automated methods. Water Resour. Res. 36:3663–3674.
Brad, Y. 1974. Nonlinear parameter estimation. Academic Press, New
York.
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous
media. Hydrology Paper 3. Colorado State University, Fort Collins.
Cooley, R.L. 1993. Exact Scheffe-type confidence interval for output
from groundwater flow models. 1. Use of hydrogeologic information. Water Resour. Res. 29:17–33.
Duan, Q., S. Sorooshian, and V. Gupta. 1992. Effective and efficientReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1352

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

global optimization for conceptual rainfall-runoff models. Water
Resour. Res. 28:1015–1031.
Duan, Q., S. Sorooshian, V. Gupta, A.N. Rousseau, and R. Turcotte.
2003. Advances in calibration of watershed models. AGU, Washington, DC.
Dubus, I.G., and C.D. Brown. 2002. Sensitivity and first-step uncertainty analyses for the preferential flow model MACRO. J. Environ. Qual. 31:227–240.
Egbert, G.D., and S.Y. Erofeeva. 2002. Efficient inverse modeling ofbarotropic ocean tides. J. Atmos. Ocean Technol. 19:183–204.
Gerke, H.H., and M.Th. van Genuchten. 1993. Evaluation of a firstorder water transfer term for variably-saturated dual-porosity flow
models. Water Resour. Res. 29:1225–1238.
Gupta, H.V., L.A. Bastidas, J.A. Vrugt, and S. Sorooshian. 2003.
Multiple criteria global optimization for watershed model calibration. p. 125–132. In Q. Duan et al. (ed.) Calibration of watershed
models. Water Science and Application 6. AGU, Washington, DC.
Gupta, H.V., S. Sorooshian, and P.O. Yapo. 1998. Toward improved
calibration of hydrologic models: Multiple and noncommensurable
measures of information. Water Resour. Res. 34:751–763.
Guyonnet, D., B. Didier-Guelorget, G. Provost, and C. Feuillet. 1998.
Accounting for water storage effects in landfill leachate modeling.
Waste Manage. Res. 16:285–295.
Hartmann, F., H.P. Bader, R. Scheidegger, and P. Baccini. 2001. Watertransport in a bottom ash landfill from a municipal solid waste(“MSW”) incinerator. J. Solid Waste Technol. Manage. 27:76–81.
Hjelmar, O. 1996. Disposal strategies for municipal solid waste incineration residues. J. Hazard. Mater. 47:345–368.
Huwe, B., H. Goelz-Huwe, and J. Eberhardt. 1994. Parameterschaetzung und Modellrechnungen zum Gebietswasserhaushalt kleiner,
heterogener Einzugsgebiete mit einfachen Modellkonzepten. Mitteilgn. Dtsch. Bodenkundl. Gesellsch. 74:273–276.
Iribar, V., J. Carrera, E. Custodio, and A. Medina. 1997. Inverse
modeling of seawater intrusion in the Llobregat delta deep aquifer.
J. Hydrol. (Amsterdam) 198:226–244.
Jarvis, N.J. 1994. The Macro model (version 4.0). Technical description
and sample simulations. Swedish University of Agriculture Sciences, Department of Soil Science, Uppsala.
Johnson, C.A., M. Kaeppeli, S. Brandenberger, A. Ulrich, and W.
Baumann. 1999. Hydrological and geochemical factors affecting
leachate composition in municipal solid waste incinerator bottom
ash. Part II: The geochemistry of leachate from Landfill Lostorf,
Switzerland. J. Contam. Hydrol. 40:239–259.
Johnson, C.A., G.A. Richner, T. Vitvar, N. Schittli, and M. Eberhard.
1998. Hydrological and geochemical factors affecting leachate composition in municipal solid waste incinerator bottom ash. Part I:
The hydrology of landfill Lostorf, Switzerland. J. Contam. Hydrol.
33:361–376.
Johnson, C.A., M.G. Schaap, and K.C. Abbaspour. 2001. Model comparison of flow through a municipal solid waste incinerator ash
landfill. J. Hydrol. (Amsterdam) 243:55–72.
Kirby, C.S., and J.D. Rimstidt. 1993. Mineralogy and surface properties of municipal solid waste ash. Environ. Sci. Technol. 27:652–660.
Konikow, L.F. 1986. Predictive accuracy of a ground-water model—
Lessons from a post audit. Ground Water 24:173–184.
Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem fortransient unsaturated flow. Water Resour. Res. 24:817–830.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. ONESTEP:
A nonlinear parameter estimation program for evaluating soil hydraulic properties from one-step outflow experiments. Bull. 85-3.
Virginia Agric. Exp. Stn., Blacksburg.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1987. Parameter
estimation for unsaturated flow and transport models—A review.
J. Hydrol. (Amsterdam) 91:255–293.
Lahkim, M.B., L.A. Garcia, and J.R. Nuckols. 1999. Stochastic modeling of exposure and risk in a contaminated heterogeneous aquifer.
2: Application of Latin hypercube sampling. Environ. Eng. Sci. 16:
329–343.
Legates, D.R., and G.J. McCabe. 1999. Evaluating the use of “good-ness-of-fit” measures in hydrologic and hydroclimatic model validation. Water Resour. Res. 35:233–241.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of
non-linear parameters. SIAM J. App. Math. 11:431–441.
McKay, M.D., R.J. Beckman, and W.J. Conover. 1979. A comparisonof three methods for selecting values of input variables in the analysis
of output from a computer code. Technometrics 21:239–245.
McKenna, S.A., L.C. Meigs, and R. Haggerty. 2001. Tracer tests in a
fractured dolomite. 3. Double-porosity, multiple-rate mass transfer
processes in convergent flow tracer tests. Water Resour. Res. 37:
1143–1154.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
Pan, L.H., and L.S. Wu. 1998. A hybrid global optimization method
for inverse estimation of hydraulic parameters: Annealing-simplex
method. Water Resour. Res. 34:2261–2269.
Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling.
1992. Numerical recipe, the art of scientific computation. 2nd ed.
Cambridge University Press, Cambridge, UK.
Savenkoff, C., A.F. Vezina, P.C. Smith, and G. Han. 2001. Summer
transports of nutrients in the Gulf of St. Lawrence estimated byinverse, modeling. Estuarine Coastal Shelf Sci. 52:565–587.
Schaap, M.G., and W. Bouten. 1996. Modeling water retention curvesof sandy soils using neural networks. Water Resour. Res. 32:3033–
3040.
Schmied, B., K.C. Abbaspour, and R. Schulin. 2000. Inverse estimationof parameters in a nitrogen model using field data. Soil Sci. Soc.
Am. J. 64:533–542.
Simunek, J., O. Wendroth, and M.Th. van Genuchten. 1999. Estimating unsaturated soil hydraulic properties from laboratory tension
disc infiltrometer experiments. Water Resour. Res. 35:2965–2979.
Sorooshian, S., Q. Duan, and V.K. Gupta. 1993. Calibration of rainfall
runoff models: Application of global optimization to the Sacramento soil moisture accounting. Water Resour. Res. 29:1185–1194.
Sulzer, S., R. Rutschmann, and W. Kinzelbach. 2002. Flood discharge
prediction using two-dimensional inverse modeling. J. Hydrol. Eng.
ASCE. 128:46–54.
van Genuchten, M.Th. 1981. Non-equilibrium transport parametersfrom miscible displacement experiments. Res. Rep. 119. U.S. Salinity Laboratory, USDA-ARS, Riverside, CA.
Vogel, T., K. Huang, R. Zhang, and M.Th. van Genuchten. 1996. The
HYDRUS 5 code for simulating one-dimensional water flow, solutetransport, and heat movement in variably-saturated media. Version
5. U.S. Salinity Laboratory, USDA-ARS, Riverside, CA.
von Gunten, U., and G. Furrer. 2000. Steady-state modeling of biogeochemical processes in columns with aquifer material: 2. Dynamics
of iron-sulfur interactions. Chem. Geol. 167:271–284.
Vrugt, J.A., W. Bouten, H. V. Gupta, J. W. Hopmans. 2003a. Toward
improved identifiability of soil hydraulic parameters: On the selection of a suitable parametric model. Available at www.vadosezone
journal.org. Vadose Zone J. 2:98–113.
Vrugt, J.A., H.V. Gupta, W. Bouten, and S. Sorooshian. 2003b. A
shuffled complex evolution metropolis algorithm for optimization
and uncertainty assessment of hydrologic model parameters. Water
Resour. Res. 39(8)1201. DOI:10.1029/2002WR001642.
Wang, Q.J. 1997. Using genetic algorithms to optimise model parameters. Environ. Model. Software 12:27–34.
Wang, W., S.P. Neuman, T. Yao, and P.J. Wierenga. 2003. Simulationof large-scale field infiltration experiments using a hierarchy ofmodels based on public, generic, and site data. Available at www.
vadosezonejournal.org. Vadose Zone J. 2:297–312.
Yapo, P.O., H.V. Gupta, and S. Sorooshian. 1998. Multi-objective
global optimization for hydrologic models. J. Hydrol. 204:83–97.
Yeh, W.W.-G. 1986. Review of the parameter identification procedures in groundwater hydrology: The inverse problem. Water Resour. Res. 22:95–108.
Young, M.H., A. Karagunduz, J. Simunek, and K.D. Pennell. 2002.
A modified upward infiltration method for characterizing soil hydraulic properties. Soil Sci. Soc. Am. J. 66:57–64.