Journal of Hydrology, 91 (1987) 255-293
Elsevier Science Publishers B.V., Amsterdam - - Printed in The Netherlands

255

[3]

PARAMETER ESTIMATION FOR UNSATURATEDTRANSPORT MODELS -- A REVIEWFLOW AND

J.B. KOOL', J.C. PARKER 1 and M.TH. VAN GENUCHTEN 2

1 Virginia Polytechnic Institute and State University, Blacksburg, VA 24061 (U.S.A.)
2 U.S. Salinity Laboratory, 4500 Glenwood Drive, Riverside, CA 92501 (U.S.A.)(Received October 27, 1986; revised and accepted November 14, 1986)

ABSTRACT

Kool, J.B., Parker, J.C. and Van Genuchten, M.Th., 1987. Parameter estimation for unsaturated
flow and transport models - - A review. J. Hydrol., 91: 255-293.
This paper reviews the current status of parameter estimation techniques and their utility for
determining key parameters affecting water flow and solute transport in the unsaturated (vadose)zone. Historically, hydraulic and transport properties of the unsaturated zone have been determined by imposing rather restrictive initial and boundary conditions so that the governing flow
and transport equations can be inverted by analytical or semi-analytical methods. Contrary tothese direct methods, parameter estimation techniques do not impose any constraints on the model,on the stipulation of initial and boundary conditions, on the constitutive relationships, or on thetreatment of inhomogeneities via deterministic or stochastic representations. While parameter
estimation analyses of subsurface saturated flow are increasingly common, their application tounsaturated flow and transport processes is a relatively new endeavor. Nevertheless, a number oflaboratory and field applications currently exist that show the great potential of parameter
estimation techniques for improved designs and analyses of vadose zone flow and transport
experiments. Several practical examples for determining unsaturated soil hydraulic functions andvarious transport parameters are presented, and advantages and limitations of the estimation
process are discussed. Specific research areas in need of future investigation are outlined.
INTRODUCTION

Increasing demands on groundwater resources have greatly accentuated theneed for accurate predictions of subsurface flow and chemical transport toevaluate effects of management practices and alternatives for contaminant
remediation. Computer simulations based on numerical models have been
increasingly used for these purposes, a trend t h a t undoubtedly will continue as
more sophisticated models are being developed and computer costs keep
decreasing. However, with greater model sophistication comes a need for more
intensive data requirements, and real improvements in precision will eventually hinge on our ability to accurately determine the required model parameters. Difficulties in model calibration are nowhere more evident than in
analyses of water and chemical transport in the unsaturated(vadose) zone.
Since most groundwatercontaminantsources originate in the unsaturated
zone, proper understanding
and managemento f t h i s z o n e is e s s e n t i a l t o

0022-1694/87/$03.50© 1987 Elsevier Science Publishers B.V.

256
protecting or improving the quality of groundwater supplies. Model calibration in the vadose zone may be especially difficult because of problems informulating accurate constitutive relationships for multiphase flow.
Hydraulic and transport properties of the unsaturated zone are commonly
determined by imposing rather restrictive initial and boundary conditions sothat the governing equations can be inverted by analytical or semi-analyticalmethods. Such procedures allow direct computation of the specific functional
form of deterministic model coefficients. A multitude of laboratory and field
methods currently exist (Klute, 1972; Bouwer and Jackson, 1974) that havefound great utility in the study of unsaturated flow and transport processes.
Unfortunately, these direct inversion methods also have a number of limitations which restrict their practicality, in particular when used for calibrating
field-scale models. Experimental analyses based on direct methods are generally quite time-consuming and hence costly owing to the need to meet conditions requisite for the explicit calculation of model coefficients. For example,
in some procedures it is necessary to repeatedly achieve steady-state or equilibrium hydraulic conditions for different boundary conditions. Often thegoverning equations themselves also need to be linearized or otherwise approximated to allow their (semi-)analytical inversion. Another limitation results from the need to impose relatively simple initial and boundary conditions.
This is especially problematic for field experiments where accurate control ofthe boundary conditions on a large scale can be difficult and expensive. Given
the inadequacy of most currently available models to describe field-scale behavior, options to deal with spatial heterogeneity in more than a superficial
manner are also severely curtailed by the direct inversion requirement. Finally, information concerning parameter uncertainty is not readily obtained
from direct inversion techniques.
An alternative and more flexible approach to solving the inverse problem is
to employ parameter estimation techniques. In this approach, the direct problem is posed for prescribed but arbitrary initial and boundary conditions, while
the physical problem itself can be solved with any appropriate analytical or
numerical method. The constitutive relations thought to be applicable to thesystem are parameterized based on a-priori knowledge, and the coefficientsdetermined by means of an optimization algorithm that extremizes some obj ective function (e.g., to minimize deviations between observed and modelpredicted output). Contrary to direct inversion methods, the optimization
approach does not put any inherent constraint on the form or complexity of themodel, on the stipulation of the initial and boundary conditions, on theconstitutive relationships, or on the treatment of inhomogeneities via deterministic or stochastic representations. Thus, a major advantage is that experimental conditions can be selected on the basis of convenience and expeditiousness, rather than by a need to simplify the mathematics of the direct inversion
process. Also, if information concerning parameter uncertainty and effects ormodel accuracy is desired, it can be obtained from the parameter estimation
analysis with little additional effort.

257
What one sacrifices for the flexibility of the parameter estimation method isan ability to determine the specific form of the constitutive properties of thesystem. One must assume some model formulation which is presumed to holdto a sufficient degree of approximation. If more than one possible parameterization is deemed feasible, the problem of parameter identification arises which
involves performing the estimation analyses for many possible models, andselecting from those models the most accurate or precise one in terms of some
objective criterion. Here, we will deal primarily with the parameter estimation
problem and only peripherally with the problem of parameter identification.
Clearly, if the adopted parametric model does not accurately represent the
behavior of the system, then results of the parameter estimation analysis will
be of dubious utility.
Other sources of parameter uncertainty are related to the ill-posedness ofmany inverse problems. A problem is ill-posed when it has either no solution
at all, no unique solution, or the solution is unstable (Isaacson and Keller,1966). In practice, ill-posedness is characterized by non-uniqueness and instability (Yeh, 1986). Non-uniqueness occurs when there are multiple parameter vectors t h a t correspond to extrema of roughly equal magnitude in theobjective function, making it impossible to determine the correct solution. As
Yeh (1986) points out, the uniqueness problem is closely related to the notionof parameter identifiability, i.e. whether or not it is possible to obtain accurate
estimates of the parameters in the mathematical model from available data.
Identifiability thus depends on both the assumed model and the experimental
data used. A frequent cause for ill-determined parameters is correlation among
them. When parameters are highly correlated, a change in one parameter isbalanced by a corresponding change in the correlated parameter, with theresult that neither can be determined accurately. On the other hand, evenwhen the parameters in a model are completely independent, the available
experimental data may lead to an objective function that lacks sensitivity to
one or more of the parameters, again with the result that these parameters will
have large estimation variances. Instability, on the other hand, occurs whenthe estimated parameters are excessively sensitive to changes in data. Relatively small measurement errors can then lead to significant errors inestimated parameter values.
The main purpose of this paper is to illustrate the application of parameter
estimation techniques to the determination of key parameters affecting flow
and transport in the unsaturated zone. Before doing so, a brief review of
relevant literature dealing with the parameter estimation process in general
will be given, followed with specific applications to the analysis of water flow
and solute transport in the vadose zone. A few areas in need of further investigation are outlined at the end of this paper.
PARAMETER ESTIMATIONMETHODSThis section briefly reviews the general formulation and solution of theparameter estimation problem. The discussion is by no means exhaustive. With258
the exception of the penalty function method, only unconstrained parameter
estimation is considered. The discussion is further limited to Newton-type
optimization techniques. Other methods, such as steepest descent and conjugate gradient methods are left out.
General formulationof the estimation problem

Many parameter estimation problems can be formulated as a weighted
least-squares minimization problem:minO(b)= ~[q* - q ( b ) ] T W [ q * -- q(b)] + ½(b* - b ) T V ( b * - b)(1)bwhere the objective function (or performance index), O(b), is a function of themodel parameters b, b = {bl . . . . . bin}T; q * = { q * , . . . , q*}Tis the observation
vector whose elements represent measured heads, water contents, concentrations or fluxes; q ( b ) = {ql (b) . . . . . q , (b)} Trepresents the predicted response fora given parameter vector b; b* represents direct estimates or measurements ofthe parameters b; and W and V are symmetric weighting matrices. The coefficient ½in eqn. (1) is purely for notational convenience. The objective is to findthe parameter vector b f t h a t minimizes eqn. (1) or, in other words, results in abest fit between the model and available data. Provided the problem has aunique solution, the final parameter values b f are also the best estimates for the
unknown model parameters b.
The second term of eqn. (1) is sometimes called a plausibility criterion (e.g.,
Carrera, 1984). As Townley (1983) remarks, inclusion of this term in eqn. (1) isequivalent to the penalty function approach for constraining parameters toremain in some feasible region around b*. When the parameter estimates b
move away from b*, the contribution of the second term of eqn. (1) to the total
weighted sum of squares will increase quadratically and minimization of O ( b )can only be achieved if b remains close to b*. The weighting matrices W and
V contain information about measurement accuracy, as well as possible correlations between measurement errors and between parameters. In the absenceof any additional information besides the observations q*, the simplest andrecommended (Beck and Arnold, 1977) approach is to set W equal to theidentity matrix and V to zero:
W = 1,

V = 0(2)

In this case, eqn. (1) reduces to the well-known ordinary least-squares (OLS)
problem:minO(b) = ~[q1 , - q ( b ) ] r [ q * - q(b)]b= ~1[q* - q(b)] 2

(3)

i-1

The OLS formulation has probably been the most popular one for parameter
estimation problems. Its attraction is due to its simplicity and because it259
requires a minimum amount of information. When observation errors arenormally distributed, are uncorrelated and have a constant variance, the OLS
estimates possess optimal statistical properties (Bard, 1974). When these conditions are not met, the OLS method will no longer yield optimal parameter
estimates in terms of precision and minimum variance. When the normality
assumption alone is violated, OLS can still be used with good results (Beck and
Arnold, 1977).
More serious difficulties arise due to violation of the constant variance anduncorrelated errors assumptions. These situations often occur in practicalproblems. For instance, error variances are commonly found to increase withthe magnitude of the property being measured. Unequal error variances also
result when the observation vector q* contains different types of measurements(e.g., pressure heads and water contents) that are expressed in different units.
In some cases the remedy is quite simple. In the case of correlated mean andvariance, a suitable transformation (e.g., logarithmic) of the measurements
may be sufficient to stabilize variances. In general, the weighting is achieved
through the matrices W and V. In what is generally called the weighted
least-squares (WLS) method, W is a diagonal matrix which corrects for unequal
error variances, and V is taken as zero. The elements o f W are chosen such t h a t
the most accurate measurements receive the most weight. The elements of W
and V may be chosen by the modeler based on knowledge and understandingof the physical system. Since subjective procedures may lead to biased parameter estimates, a preferred method is to base the weighting on statisticalconsiderations. This leads to selection of the inverse of the error covariancematrix as the weighting matrix:
w=

C~ ~

Cq = E[(q* - ~)(q* - ~ ) T ]
4

=

(4)

E(q*)

where E denotes expectation. The diagonal elements of W account for unequal
error variances, and off-diagonal elements for correlated errors. In case theoff-diagonal elements of W are also taken into account, the method is classified
as generalized least-squares (GLS). Often, when a method is described as WLS
or GLS, no prior information on the parameters themselves is presumed and theobjective function consists of only the first term of eqn. (1).
Although the formulation of the parameter estimation problem as aweighted least-squares problem can be arrived at by heuristic arguments,
casting the problem in the framework of either maximum likelihood orBayesian estimation in practice also leads to objective function (1) (Townley,
1983; Carrera, 1984). Whereas the underlying concepts are quite different, eachof the latter two methods leads to the same estimator. The maximum likelihood(ML) interpretation is more commonly used, because the ideas it involves are
more straightforward (Beck and Arnold, 1977; Townley, 1983). In the ML

26O
method, model parameters are viewed as unknown but deterministic and theobjective is to find the parameters that maximize the likelihood of obtaining
the measured data, given the joint probability density function of all measurements. A common approach is again to assume errors are normally distributed.
This leads to the convenient form of eqn. (1) for the objective function with W
and V given by:W=

C ; 1,

V

=

(5)

C; 1

where Cq was defined previously and Cb is analogously given by:
Cb=

E[(b -

6=E(b)5)(b -

b')T]

i (6)

I

Hence, W and V are inverses of the covariance matrices of measurement errors
and prior parameter estimate errors, respectively.
If a distribution of errors other than the normal distribution is assumed, theresulting objective function will not have the generalized least-squares form ofeqn. (1). Although the true error distribution is frequently unknown, violationof the normality assumption will generally not lead to serious errors (Bard,
1974). Maximum likelihood estimates are asymptotically (i.e. for large samples)
unbiased, normally distributed and have minimum variance (Bard, 1974). To
achieve the full benefit of the ML method, the covariance matrices Cq and Cb
need to be specified. Indeed, Beck and Arnold (1977) recommend simple OLS
estimation if these matrices are completely unknown. In general, the weighting
matrices will not be known a-priori. In some cases the appropriate weights inthe objective function can be estimated along with the model parameters.
Examples are given by Bard (1974), Beck and Arnold (1977), Sorooshian andDracup (1980), and Troutman (1985). A generally applicable procedure is thetwo-stage least-squares method recommended by Cooley (1982) and Sadeghipour and Yeh (1984). The first stage consists of OLS estimation which utilizes
only the available observations q*. The OLS solution will generally not lead to
optimal parameter estimates, but it will provide unbiased estimates of the error
covariances. Let e be the residual vector at the OLS minimum:O(b f) : [q* - q ( b f ) ] T [ q * -- q(bf)] = e T e(7)
An unbiased estimate of Cq is then given by:
Cq-1
n - 1 eeT(8)

where n is the number of observations. To further improve the parameter
estimates, the second stage of the procedure solves the generalized leastsquares problem using the estimated Cq to provide weights. Based on a numberof numerical experiments involving estimation of aquifer transmissivities fromhydraulic head data, Sadeghipour and Yeh (1984) concluded that the two-stage
approach as compared with OLS estimates led to significantly better parameter
estimates.

261

In cases where prior information on the model parameters is available, thetwo-stage procedure may be extended to include a plausibility criterion. The
parameter error covariance matrix Cb may be estimated at the end of the firststage from:
Cb -

eTe

- n

--

H -1

(9)

m

where m is the dimension of the parameter vector and H is the Hessian, or
matrix of second derivatives, of the objective function with respect to theparameters with components:
Hij -

620
~bi~bj(10)

Depending on the algorithm for solving eqn. (1), H will already be available orcan be easily computed at the end of the first stage estimation• Note thatregardless of whether or not the second stage estimation is carried out, eqn. (9)may be used to provide an estimate of parameter uncertainty. As an overall
index of inversion precision, a norm such as the trace or maximum eigenvalueof Cb may be employed (Yeh, 1986).
Solution methods for least-squares problems
In this section Newton's and related methods for solving least-squares minimization problems will be reviewed. The starting point is the objective function, given by eqn. (1) in the previous section. In order to simplify the notation
it is assumed t h a t the weighting matrices W and V are diagonal. This allows
a new vector of residuals r to be written, where:Wl[q~l - - q l ( b ) ]rlrn

w~ [qn* -- qn(b)]
=

(n)

r~ ÷

vl (b* - bl)r,+vm(b* - bin)

with:
w,

=

(w~,) '~2,

vj

=

(vJ

~

Unequally weighted measurements and prior parameter information are all
lumped into r, so t h a t eqn. (1) becomes:

262NminO(b) = ~rTr = ½~ ~
b

(12)

i=l

with N = n + m. If W and V are not diagonal, they must be propagatedexplicitly in subsequent equations. While this is a straightforward procedure(e.g., Beck and Arnold, 1977; Townley, 1983) the equations become lengthierand salient features of the different solution methods are less clear. Unless themodel is linear in all parameters, minimization of eqn. (12) must be carried out
iteratively. Initial values, b °, for the parameters must be supplied by the user.
At every iteration, a linear system of equations is constructed and solved togive a correction Ab for the current parameter vector b i, such that:O(b i + Ab) ~< O(b ~)(13)

The iteration continues until some suitable convergence criterion is satisfied.
Commonly used stopping criteria are:O(b ~+1) - O ( b ~) <~ ~(14)

and/or:Ab i <~ "c2where T~ and z2 are set to some small value.
Newton's method can be derived by writing a 3-term Taylor series expansion
for O(b) around bi:O(b~ + Ab) = O(b ~) + VO(bi)TAb(15)+ ~AbTV20(bi)Ab + "small"where "small" contains terms of the order of IIAbll3 and smaller. We wish to
select Ab so t h a t O(b i + Ab) is approximately minimized. The correction Ab isdetermined such that:(l)(Ab) = VOTAb + ½AbTV2OAb(16)

has a minimum. A necessary condition for Ab to be a minimizer is that thegradient of ¢(Ab) is zero:
VO + V2OAb = 0(17)

and thus:
Ab =-(V20) 'VO(18)

First and second derivatives VO and V20 are obtained by differentiating eqn.
(12) with respect to b:
N

VO

and:

=

E riVri
i 1

=

JTr

(19)

263
N

V20 = H=~(VriVr T + riV 2ri) = j T j

+ S

(20)

i=l

where J is called the Jacobian or sensitivity matrix and H is the Hessian (seeeqn. 10). Substitution of eqns. (19) and (20) into eqn. (18) gives Newton'salgorithm for updating the unknown parameters b:
Ab

=

-

(jwj

+ S)-IjTr(21)

In the so-called full-Newton method, the Hessian (eqn. 20) is evaluated directly.
This method is (at least locally) quadratically convergent (Dennis and Schnabel, 1983) and will quickly converge for most problems. A major disadvantageof the full-Newton method is the cost of evaluating the second derivatives:
V2rk -~2rk~bi~bj(22)

These second derivatives are generally not available analytically. Because
approximation by finite differences requires (rn2 + 3m)/2 + 1 function evaluations for m parameters, this can greatly add to the overall computational
expense. The full-Newton method is for this reason, not commonly used forsolving non-linear least-squares problems. More popular are modifications of
Newton's method t h a t do not calculate eqn. (20) directly, but use someapproximation thereof. These modifications are known as quasi-Newton,
Gauss-Newton and Levenberg-Marquardt, respectively. Closest to the fullNewton method is the quasi-Newton method which approximates the Hessian
from available first derivative information using a secant update method. A
positive definite matrix is constructed and updated at every iteration of the
optimization algorithm so t h a t it yields an approximation to H or H - 1 in far
fewer steps t h a n a full Hessian evaluation requires. Positive definiteness of the
updated matrix ensures that the optimization will proceed in a descending
direction. Well-known general methods for updating the matrix include theDavidon-Fletcher-Powell (DFP) and Broyden-Fletcher-Goldfarb-Shanno(BFGS) methods. Equation (20) shows t h a t the least-squares Hessian consistsof a first derivative term j T j and a second derivative term S. Dennis et al.
(1981) describe a quasi-Newton algorithm t h a t takes into account this special
structure of the least-squares Hessian by approximating only the second derivative term in eqn. (20) with a secant update, and adding this to the directly
evaluated first derivative term.
In the Gauss-Newton method it is simply assumed that the first term j T jabove is a sufficient approximation to the Hessian. Comparison with eqn. (19)shows t h a t this is a reasonable assumption when the residual vector r is small.
In fact when the residuals are zero at the minimum, the Gauss-Newton method
is asymptotic to Newton's method near the minimum and thus locally quadratically convergent. On the other hand, when the solution is far from the minimum or when the residuals remain large at the minimum, the Gauss-Newton
method may fail to converge at all. Another problem with the Gauss-Newton

264

method is its sensitivity to rank deficiency in the Jacobian J which occurs
when one or more of the columns of J are linear combinations of other columns.
This rank deficiency is related to the identifiability problem mentioned in the
introduction. Identifiability of the model parameters requires t h a t the Hessian,given by eqn. (20), is non-singular. Otherwise, eqn. (18) will not have a uniquesolution. For the Gauss-Newton method where H is approximated by j T j , anon-singular Hessian is equivalent to a full rank Jacobian. Exact rank deficiency in the Jacobian occurs when two or more parameters in the model arelinearly dependent. Such overparameterization can usually be detected andcorrected by inspection of the model that is fitted. Because of finite precision
arithmetic, problems will often arise in practice when parameters are highly
correlated, e.g., with correlation coefficients exceeding 0.90. This leads to apoorly conditioned Jacobian and slow convergence of the optimization routine.
Numerical problems can be alleviated to some extent by using robust techniques (such as orthogonalization) for solving the linear least-squares problem(21), but the final parameter estimates may still have large variances. Thislatter problem can only be eliminated by choosing a simpler model and/or by
including additional information in the objective function.
In the Levenberg-Marquardt method the Hessian is approximated as:
H

~ jTj

+ 4DTD(23)

where 4 is a positive scalar and D a diagonal scaling matrix. The elements of
D are usually set equal to the norms of the corresponding columns of J (Mor6,
1977). The Levenberg-Marquardt method has a number of advantages over the
standard Gauss-Newton method. For non-zero 4, the Hessian approximation isalways positive definite which ensures the descent property of the algorithm,
even if initial parameter estimates are poor. The Levenberg-Marquardt method can be viewed as an interpolation between the steepest descent and Gauss-Newton methods. When ,t is large, the result will be a small step in thesteepest-descent direction, while the method degenerates to the Gauss-Newton
method as 4 approaches zero. The usual strategy is to set 4 initially to somelarge value, and then to decrease its value as the solution approaches theminimum. In principle, the method is also numerically more robust than theGauss-Newton method. For non-zero 4, j T j + 4DTD will have full column
rank even if j T j itself does not. In the Gauss-Newton and Levenberg-Marquardt methods, the approximations j T j and are often used to compute H - 1 in
eqn. (9). This first-order approximation will be accurate to the extent that thesecond term in eqn. (20) is negligible. The alternative of course is to compute
H in eqn. (9) directly.
All of the methods discussed are expected to work well for many problems.
Since parameter estimation problems involving transient flow and/or transport
are computationally expensive to solve, it will in many cases be worthwhile tofind the most efficient method rather t h a n just use one that is already available.
Different optimization algorithms are widely available through software
libraries such as IMSL and MINPACK. Comparisons of different algorithms (e.g.,
Bard, 1974; Beck and Arnold, 1977; Hiebert, 1981; Cooley, 1985) indicate that
results can be very problem dependent. While the quasi-Newton method is265
preferred for general function minimization, there appears to be little
difference between it and Levenberg-Marquardt methods for least-squares
problems (Hiebert, 1981; Dennis and Schnabel, 1983). In the problems considered by Hiebert (1981) the Gauss-Newton method performed least well.
Hiebert further concludes that the specific implementation of an algorithm cansignificantly affectperformance, and discusses factors that affect efficiency androbustness.
Besides using an efficient optimization algorithm, i.e., one that converges in
the least number of iterations, a reduction in overall computational cost canbe achieved by reducing cost per iteration. The most expensive part in thesolution process is often the evaluation of the Jacobian, i.e., calculation of theparameter sensitivities. These sensitivities are usually evaluated by finite
difference approximation. Even when one-sided differences are used, the direct
problem must be evaluated m + 1 times, in order to obtain the Jacobian for m
parameters. Recently, a number of researchers (Neuman, 1980; Carrera and
Neuman, 1984, 1986; Townley and Wilson, 1985) have applied adjoint techniques for obtaining the gradient of the objective function in cases where thedirect problem is given by a differential equation (e.g., the non-steady flow
equation). The technique involves formulation of an adjoint equation to thedirect problem. The gradient VO sensitivities can then be obtained by solving
both direct and adjoint problems once, at a cost that is approximately the same
as solving the direct problem twice. Net savings thus result as soon as there ismore than one parameter to be estimated. A disadvantage of the adjoint method
is that one obtains the gradient VO, but not the sensitivity matrix J. Consequently, no direct approximation of the Hessian is available, and one either has
to resort to less powerful optimization methods than the Gauss-Newton orLevenberg-Marquardt method, or combine the adjoint method with a secant
Hessian approximation to yield a quasi-Newton algorithm. The adjoint method
has thus far been applied to saturated flow problems. A first application tounsteady, unsaturated flow problems is described by Wittmeyer and Neuman(1985).
APPLICATIONS TO UNSATURATEDFLOW
Several applications of the parameter estimation method to the determination of unsaturated soil hydraulic properties are discussed in this section.
First, the direct problem governed by the saturated-unsaturated flow equationand suitable constitutive relationships for the hydraulic properties will bedefined. This is followed by discussion of laboratory and in-situ methods.
The direct problem

The partial differential equation governing transient one-dimensional unsaturated flow in a rigid porous medium is taken as:~hC(h) ~t-~ [ K ( h ) ~ ( h + z).]

~x

~x

(24)

266

where x is distance; t is time; K ( h ) is the hydraulic conductivity as a functionof the pressure head h; C ( h ) is the soil water capacity, being the slope dO/dh of
the soil water retention curve 0(h), where 0 is the volumetric water content;
and z is the gravitational head such t h a t d z / d x = 0 for horizontal flow, + 1 forvertical flow with x positive upwards and - 1 for vertical flow with x positive
downwards. The unsaturated hydraulic properties O(h) and K ( h ) are strongly
non-linear functions of the pressure head. It is assumed that suitable analyticalexpressions for these functions are available. Strictly, the problem is one of themodel identification as well as parameter estimation since the correct functional forms of 0(h) and K ( h ) are generally unknown. Frequently used expressions
for the hydraulic functions are those developed by Brooks and Corey (1964) and
Van Genuchten (1980). The expressions derived by Van Genuchten are:
So =

g

=

1 + Ie h l n ) ~'

h < 0

1,

h >~ 0K8S~/2[1 - (1 - Se,/mm
) ]

(25a)

m

=

1 -

1In

(25b)

where Se is the effective saturation (0 - 0r)/(0s - 0r); 08 and 0r are saturated
and residual water contents, respectively; Ks is the saturated conductivity andand n are curve shape parameters. Note t h a t m = 1 - 1/n is not an independent parameter. An expression for C ( h ) = dO/dh is readily obtained by differentiating eqn. (25a).
As written here, eqn. (25) is valid for monotonic wetting or drying only.
When the flow process involves both wetting and drying, hysteresis in the C(h)relation will have to be taken into account. Kool and Parker (1987) recently
combined Van Genuchten's expressions (25) with an empirical hysteresis model
to obtain a set of concise, closed-form expressions describing hysteretic soil
hydraulic properties with only one additional parameter. The model wasdeveloped with an eye towards application in parameter estimation problems(Kool et al., 1986) where parsimony in the number of model parameters isdesirable.
Parameters in eqn. (25) are commonly determined by measuring Ks directly
and estimating the remaining parameters by fitting eqn. (25a) to measured (0, h)data. The (0, h) measurements are obtained from equilibrium desorption orabsorption experiments on laboratory samples or can be measured in-situ using
tensiometers and neutron probe. The parameter estimation problem in this casebecomes very simple. Disadvantages of such an approach are generally timeconsuming data collection and the fact that parameters are fitted to (0, h) data
only, so t h a t any inaccuracy in the assumed hydraulic relationships, as well as
effects of measurement error are forced into the predicted K ( h ) (Parker et al.,1985). Rather than directly measuring only pairs of (0, h) data, a number ofauthors have recently tried to estimate the hydraulic properties from transient
flow experiments by numerical inversion of eqn. (24). Parameter estimation
studies of this type have generally been limited to relatively short laboratory

267
soil columns, mostly by measuring outflow volumes from initially saturated
columns. A comparatively few field-scale experiments involving in-situ
measurements of water contents and pressure heads during transient flow havealso been attempted. Laboratory approaches will be discussed first.

Laboratory methods
Probably the first attempts to estimate the hydraulic properties from transient flow experiments in the laboratory were those by Zachmann et al. (1981,1982) who considered hypothetical gravity drainage experiments involving
initially saturated columns of sand. Cumulative drainage outflow data wereused to estimate by ordinary least-squares two unknown coefficients in afour-parameter model for the hydraulic properties. The direct problem wassolved using finite differences, while the objective function for the 1982 study
was given by eqn. (3) with the observation vector q* representing cumulative
outflow volumes (Zachmann et al., 1982). Their study shows the importance of
selecting correct parametric forms for K(h) and O(h) in the estimation process.
When incorrect expressions are used, it may still be possible to obtain an
acceptable solution to the inverse problem but the hydraulic properties corresponding to that solution may be in error.
This last problem is illustrated in Fig. 1, adapted from Zachmann et al.
(1982), in which both the "true" and estimated O(h) and K(h) relationships areshown. The "true" parametric functions were represented by Brooks andCorey-type expressions, while in the inverse problem an exponentional func10 ~

.C

E

I

*

/-

o -'

IO -z

10

-z

10

-3

•

,o-,

/a

°

O. O 0O. ;258(m3m -3)

t
0. 50

b

l O -~,11111

O. O 0O. 2 5

8

O. 5 0(m3m -3)

Fig. 1. Actual and predicted hydraulic properties for a hypothetical porous medium (after
Zachmann et al., 1982). (a) water retention curve 0(h); (b) hydraulic conductivity K(O).

268
tion was used for 0(h) and a power function for K(O). Although our own
experience, as well as other published results (e.g., Van Genuchten andNielsen, 1985), suggest that Van Genuchten's expressions (25) are reasonably
accurate for many soils, the appropriateness of this or any other model for aparticular case will not be known beforehand, and it may be desirable tocompare the performance of different functional forms for the hydraulic functions. An obvious criterion for model discrimination is to select the model thathas the minimum residual error, under the subjective constraint that the model
should predict realistic properties for the medium. A more rigorous criterion
involves evaluation of how well a fitted set of hydraulic properties can describe
system behavior under initial and boundary conditions t h a t are different fromthose used in model calibration.
Another common cause for erroneous parameter estimates, even when acorrect parametric model is used, occurs when different combinations of parameters lead to roughly the same minimum in the objective function. This typeof uniqueness problem was investigated by Hornung (1983) and Kool et al.
(1985). In both studies, soil hydraulic properties were represented by Van
Genuchten's model. Hornung assumed the parameters Ks, 0s and 0r to be known
and attempted to determine a and n from a column drainage experiment. A
vertical column, initially at hydrostatic equilibrium, was subjected to a constant infiltration flux at the surface and the column outflow rate measured. No
unique solution could be found using only these measurements. The nonuniqueness problem could be resolved when additional information was included, in Hornung's case the final steady pressure head at some fixed position
inside the column. Kool et al. (1985) estimated the three parameters 0 , a andn in eqn. (25) also from a column drainage experiment. Unlike the gravity
drainage studies of Zachmann et al. (1982) and Hornung (1983), the experimental set-up in their study was that of a "one-step" pressure outflow method where
desorption was initiated by a step change in air pressure at the upper column
boundary. Kool et al. (1985) concluded that the ability to uniquely solve thethree-parameter estimation problem is a function of the difference between the
average initial and final water content in the column. Uniqueness was not a
problem when this difference was at least 50% of the difference between 0s and
0r- This suggests that in gravity drainage experiments the desired resolution
can only be achieved by using unrealistically long soil columns, except for
extremely coarse materials such as those used by Zachmann et al. (1982). The
one-step outflow method is much more convenient to use since the watercontent difference depends on the imposed step change in potential, which is
directly controlled by the experimenter.
Parker et al. (1985) used the one-step outflow procedure to determine hydraulic properties of undisturbed soil cores taken from four different field soils
ranging from sandy loam to clay. Soil cores were assembled in pressure cells,
saturated and then subjected to a step change in the pressure head, Ah, of 10m
by increasing the external gas pressure. The parameters 0s and Ks were independently measured, while a, n and 0r were estimated from the cumulative

269

outflow versus time measurements. Predicted O(h) and K(h) relationships werethen compared with independently determined hydraulic properties for thesame soil cores. Good agreement was found for all soils in the range
l0 < h < 0 m, i.e. the same range of pressure heads over which the outflow
measurements were taken, but extrapolation to lower pressure heads outside
the measurement ranges proved to be less reliable. Predicted O(h) curves
generally overestimated water contents at pressure heads < - 10 m, indicatingpoor estimates of 0r. This problem could be alleviated by including additional
information in the objective function, in this case a direct measurement of 0 at
h = -150m. Thus, their objective function was of the form:-

N

M

O(b) = w ~ [q*- qi(b)]2 + v ~ [0"- Oj(b)]2
i-1

(26)

j-1

where q ' i s measured cumulative outflow at times ti (i = 1 , . . . , N); qi(b)represent outflow at times t~, computed by solving eqn. (24) for the trial parameter vector b with b having elements (~, n, 0r)T; 0* is the measured watercontent at pressure head hj; and Oj(b) is the corresponding predicted watercontent obtained from eqn. (25a). Measurement errors were assumed to be
uncorrelated and to have constant variance. Weights w and v account for thedifferent scales of measurement for outflow volumes and water contents. Since
only the ratio w/v is of interest, w was set equal to 1.0 and v calculated as:NM

v = M ~ q*/N~ 0*
i-1

j

(27)

1

The one-step outflow method is a convenient procedure for obtaining input
data for the inverse problem because it is experimentally simple and quick,
Jennings et al. (1985) describe another parameter estimation procedure basedon a pressure outflow test. Their experiment involves a stepwise desorption testin which each step is not necessarily continued until equilibrium. Jennings et
al. use spline functions to describe the retention and relative permeability
functions, with " k n o t s " on the spline acting as parameters obtained by matching the observed stepwise outflow curve.
A disadvantage of these and most other parameter estimation procedures
based on transient, unsaturated flow data, is the computational expense associated with the numerical solution of eqn. (24). During iterative optimization,
the direct problem given by eqn. (24) may have to be solved 30 or more times.
As discussed in the previous section, a large part of total computational costis associated with the finite difference approximation of the Jacobian. A reduction in computer costs can likely be obtained by using adjoint-state techniques
for evaluating the gradient of the objective function (Wittmeyer and Neuman,
1985). An alternative is to use an experimental design that permits a (semi-)
analytical solution for the direct problem. For horizontal flow eqn. (24) can be
written as:

270

0-~

D h

~xx

-

~t

(28)

where Dh = Kdh/dO is the hydraulic diffusivity. Philip (1955) gave a simple
iterative solution of eqn. (28) for adsorption into a semi-infinite column subject
to a first-type inlet boundary condition. The solution was recently modified by
Van Grinsven et al. (1985) for desorption. Using these solutions, eqn. (28) canbe solved much more cheaply than eqn. (24). Input data for the inverse problem
could consist of measurements of cumulative absorption (desorption), and/or
moisture content distribution inside the column at one or more times. W.
Bouten (University of Amsterdam, Netherlands, unpublished results) applied
the modified Philip solution to estimate D h from cumulative evaporation measurements. Using Van Genuchten's expressions, the parametric form for Dh (0),written for simplicity in terms of the effective saturation Se, is:D~(Se)(1 - m)K~ 1/2 1/m.ql/m)-m.ql/m~m~m~s---~Se[(1--_e ,+ ( 1 - - _ o i -- 2](29)

with m = 1 - 1/n. Note that Ks and a appear in eqn. (29) as a ratio andtherefore cannot be determined simultaneously from horizontal flow experiments only. For laboratory soil columns, direct measurement of K~ is a simple
procedure, so t h a t this should not present a problem. Richie et al. (1987) outline
a parameter estimation procedure for determining hydraulic properties frommeasured pressure head distributions along an unsaturated horizontal soilcolumn during steady-state flow, and using additional measurements of total
mass of water in the column, steady flow rate and saturated water content.
Depending on the functional relationships for O(h)and K(h), the direct problem
is given in closed or semi-closed form, and the inverse problem could thus besolved very cheaply. Richie et al. (1987) obtained very good results for two
hypothetical soils but did not test the procedure on real soils. A disadvantageof their experimental set-up is that it requires a considerable number of tensiometers to be installed in the column. In the hypothetical experiments 16
tensiometers were installed along a 30 cm long column. As compared with thetransient outflow procedures of Kool et al. (1985) and Jennings et al. (1985), themethod of Richie et al. (1987), as well as those based on eqn. (30), sacrifice someexperimental convenience for reduced computational expense.

In-situ methods
Whereas laboratory experiments have the advantage of being easy, quick
and precise, a major disadvantage is that they lead to soil properties that areoften non-representative of field conditions. Main reasons for this are thetypically small size of laboratory samples, and the fact that collection of soilcores invariably introduces some disturbance that may affect flow and transport properties. Since analyses of the soil hydraulic functions are ultimately
directed towards field-scale processes, determination of in-situ properties is271
more relevant t h a n data obtained from laboratory analyses. The application ofparameter estimation techniques to determination of in-situ soil hydraulicproperties is, at least in principle, a straightforward matter. Restricting ourselves to vertical one-dimensional flow, input data for the inverse problem canconsist of measured water contents and/or pressure heads at different depths
and times during infiltration and/or drainage events. Solution of the direct
problem (24) further requires the stipulation of suitable initial and boundary
conditions. Having instrumented the site with equipment for measuring watercontents and/or pressure heads, initial conditions can be determined directly,
while boundary conditions can be treated as additional unknowns in theestimation problem as needed.
To date, only one study has been published (Dane and Hruska, 1983) in whichparameter estimation methods were employed to determine the water retention
and hydraulic conductivity functions simultaneously from a transient field
experiment. In this study, Dane and Hruska studied gravity drainage from aclay loam soil. Initial water content profiles and profiles after 7 and 25 days ofdrainage were measured with a neutron probe. The soil surface was covered toprovide a zero-flux upper boundary condition, while a first-type lower boundary
condition was determined from tensiometer readings at a depth of 0.9 m. Parameters ~ and n in Van Genuchten's model (25) were determined by numericalinversion of eqn. (24) while other parameters in the model were assumed
known. 0s and K S were determined from the maximum water content andsteady-state infiltration rate, respectively, both measured during wetting of thesoil prior to drainage. A guessed value, obtained by inspection of available 0 (h)data, was assigned to 0r. Dane and Hruska obtained good agreement between
the predicted and independently obtained O(h) relation for the clay loam soil,
but found t h a t the predicted K(O) relation overestimated independently determined conductivities by approximately one order of magnitude. Much better
agreement was achieved when the value of Ks was arbitrarily taken to be 10
times lower than its measured value; this produced only small changes in the
optimum values for ~ and n. Use of a much lower than measured Ks value wasjustified by arguing that macropore flow probably occurred during ponded
infiltration, leading to an inflated value for Ks.
A more general procedure than that followed by Dane and Hruska (1983)
would be to not fix 0 r and Ks, but to take these as additional unknowns in the
inversion problem. As indicated by the results of Dane and Hruska and others(see the discussion by Van Genuchten and Nielsen, 1985), a reliable value forK~ is difficult to obtain, especially in structured media. A value for 0~ could beobtained by determining the water content at some suitably low pressure head,e.g., - 150 m. This is a somewhat arbitrary definition of Or. However, experience
shows that when input data for the inverse problem do not include measurements at very low pressure heads, as is the case for gravity drainage, inversion
results are not very sensitive to the value of 0r. On the other hand, it also
implies t h a t predicted O(h) and K(h) relations cannot be reliably extrapolated
to dry conditions.

272
Next, an example of par a m et er estimation for a drainage problem similar tothe study of Dane and H r us ka (1983) is considered, but with Ks and 0r t aken as
additional unknowns. The example involves drainage from a 3 m diameter by
6 m deep lysimeter, filled with crushed Bandelier Tuff, a material with silty
sand texture. The in-situ experiment was carried out at the Los Alamos
National L a b o r a t o r y and is described in detail by Abeele (1984). The lysimeter
was instrumented with n e u t r o n probe access tubes and tensiometers at depths
of 0.4, 1.16, 1.91, 2.71, 3.47 and 4.23m. The average water content following
ponded infiltration for more t h a n one month provided a value for 0s of 0.331.
The saturated conductivity, Ks, computed from the steady-state flux exciting
the bottom of the lysimeter was 0.1244 m day 1. After saturation, the lysimeter
was allowed to drain for 100 days during which period water content andpressure head profiles were monitored. The surface of the lysimeter wascovered during the drainage period. In the present analysis of the problem Van
Genuchten's model (25) is used to describe hydraulic properties with 0s the only
known parameter.
Preliminary analyses of hypothetical data sets showed t hat the parameter
estimation problem for the two parameters a and n can be solved uniquely using
only information on water cont ent profiles during drainage, but t hat thesimultaneous estimation of three or more parameters requires additional information. For this reason, measured water contents, and additionally, measured pressure head at one depth (0.4 m) were used as input data for the inverse
problem. Observed water contents and pressure head at t = 1, 4, 10, 20, 40 and
100 days after the start of drainage were used in the objective function. Observed pressure heads at depth x = 4.23m were used to provide the lower
boundary condition for the solution of the direct problem. Observed watercontents at this depth were not used, so that simultaneous information onwater contents and heads at x = 0.4m only was employed. The initial condition was h(x) = 0, while the upper boundary condition was t hat of a zero flux.
The direct problem (24) was solved using an efficient fully implicit, mass-lumped
Galerkin-type, linear finite-element code with a variable time step and constantnode spacing Ax of 5 cm. As a check on numerical accuracy, a computer run forthe same problem with Ax = 2.5 cm was made, which gave essentially identicalresults. The inverse problem was formulated as a weighted least-squares
problem with suitable weighting to account for the different scale of measurement of water contents and pressure head:
5

O(b)

=

6

6

~ ~ [0" - Oij(bl] 2 + v ~ [ h * - hi(bl] 2
i-lj

1

j

(30)

1

where 0* represents measured water contents at depths xi and times tj, h* is themeasured pressure head at depth 0.4 m and times tj, and b is the four-parameter
vector (a, n, 0r, Ks) T. After P a r k e r et al. (1985), relative weights for pressure
heads are computed as the ratio of mean observed water content to meanobserved pressure head:

273
5

6

6

i=lj~lj=l

v = 6 ~ ~ 0*(xi, tj)/30 ~ h*(ti)(31)

With pressure heads expressed in m H20, this gave a weighting of v = 0.160.
The non-linear least-squares problem for the four parameters (~, n, 0r, Ks) was
solved using the Levenberg-Marquardt method. Error variances were computed with eqn. (8) as a check on the correctness of the weighting coefficient(31). Weighting the residuals in the objective function by reciprocals of error
standard deviations would have given a relative weight of v = 0.124 to thepressure head observations. This was considered sufficiently close to the actual
weight so that the analysis was not repeated.
Estimated parameter values and their statistics are given in Table 1. Since
the estimated value for 0r is identically zero, there are no error estimates forthis parameter. Of the remaining parameters, the estimated Ks is the leastaccurate with a 95% confidence interval of _+50% of the estimated value.
Observed and fitted water content profiles at t = 1, 10, 40 and 100days areshown in Fig. 2. For pressure heads at 0.4 m depth, the mean deviation between
observed and fitted values was 0.07 m H20. Fig. 2 shows good overall agreement
between observed and fitted water contents, with poorer results for t = 1 day
and for the 3.47 depth at t = 4 days. Computed water content profiles show amore rapid initial decrease than observed. This is related to the high estimatedvalue for K s (Table 1) which is about 2 times higher than the observed value of0.124 m day -1. A too high Ks value will result in overly rapid initial drainage,
but will have less of an effect where the soil has become unsaturated. Consequently, the estimated value of Ks is determined mainly by observations at
small times. The wide 95% confidence region further indicates the poor identifiability of Ks. As discussed previously, the remedy for this is to include priorinformation on Ks in the least-squares function. In practice, determining a goodvalue for in-situ Ks is quite difficult. Alternatively, the early time observations
could have been assigned higher weights to obtain a more accurate estimate for
Ks. If an experiment contains no observations at saturated or near-saturated
conditions, Ks will be very ill-determined, as for example in the experiment ofDane and Hruska (1983) who observed little sensitivity to K s when it changedby an order of magnitude. To verify accuracy of parameter estimates, predicted
TABLE 1
V a n G e n u c h t e n p a r a m e t e r s for Bandelier Tuff estimated from in-situ d r a i n a g e
Parameter(m -1)n0r
K s (m day -1)

Initial
value

2.00
2.00
0.05
0.25

Estimated
value

1.433
1.506
0.0
0.250

Std.
dev.

0.15
0.052
0.063

95% confidence limits
lower

upper

1.13
1.40
0.124

1,74
1.61
0.376

274
0. O
o

1.O

v
Ii
~

2"O

[]

3.0

4.0
0.2

O.

O. 3

8

0.4

(m3m -3)

Fig. 2. Observed and fitted water content profiles during drainage of Bandelier Tuff. Observed datafrom Abeele (1984).

10l O -1

O

10 'j
10 -$

v

1'-

10

b

o

X
]O

-1
10

-4

10

-5

b
1O
O.C)

(7]

0.2

8

(rn3rn-3)

0.3

0.4

o.n

0.1

0.2

B

O.N

0.4

(m~m -3)

Fig. 3. Hydraulic properties of Bandelier Tuff. (a) water retention curve 0(h); (b) hydraulic conductivity K(O). Data points from Abeele (1984), * is the measured saturated conductivity.

275
and independently determined O(h) and K(O) relations for the Bandelier Tuffare compared in Fig. 3. The 0(h) data points in Fig. 3a represent in-situ measurements from the lysimeter as well as from the laboratory at the lower
pressure heads. The K(O) data points shown in Fig. 3b were calculated from thelysimeter drainage experiment (Abeele, 1984) using the instantaneous profile
method (Watson, 1966). Both figures show excellent agreement.
The results of Fig. 3 as well as those by Dane and Hruska (1983) clearly
illustrate the potential of parameter estimation techniques for determinationof in-situ soil hydraulic properties. We acknowledge that the lysimeter experiment was carried out in a very homogeneous medium, whereas field soils may
exhibit marked heterogeneity (e.g., layering) at this scale of observation.
Nevertheless, the type, intensity and accuracy of measurement in the lysimeter
example are those of a typical field situation. In this case, measurements weretaken over a 100 day period. While this is quite long, the frequency with whichmeasurements need to be taken should decrease quickly as the duration of theexperiment increases. The lysimeter data changed little after approximately 20
days. In the present study no attempt was made to determine the minimum
required experimental duration. Experience with laboratory outflow tests indicate that the important factor is change in water content and pressure head
between beginning and end of the experiment rather than duration of theexperiment per se. This means that measurements would have to be taken overa longer period of time for a slowly draining heavy-textured soil as comparedwith a more quickly draining light-textured soil. Alternatively, the desired
effect may be obtained more rapidly by not covering the soil surface duringdrainage and allowing simultaneous evaporation and gravity drainage. Unfortunately, unless the generally variable evaporation rate is known a-priori, thiswill lead to additional unknowns in the parameter estimation process. A
practical solution would be to take periodic measurements of the surface
moisture content during the experiment and use this information to obtain aprescribed first-type surface boundary condition. Another alternative would be
to use transient infiltration data in the inversion problem, either instead of orin addition to drainage data. An infiltration test, especially when the soil is
initially dry and infiltration is continued until the medium is saturated, should
lead to better resolution in the input data, and hence to more accurate estimates for Ks. The simultaneous use of both infiltration and drainage data also
leads to a better definition of the parameters. As shown by Kool et al. (1986) for
laboratory data, an additional advantage of the latter approach is that it canbe used quite effectively to determine the soil's hysteretic hydraulic properties.
APPLICATIONS TO SOLUTE TRANSPORTThe common representation of subsurface transport employs the convection~lispersion equation t h a t includes terms for solid-liquid partitioning andirreversible kinetics. Thus, the general one-dimensional transport equation istaken to be of the form:

276
~t +

8t

-

8x

OD-~x - qcr- Opcr + 07

(32)

where Cr is the resident solution concentration, s is the adsorbed concentration,p is the dry soil bulk density, D is a dispersion coefficient, q is the Darcian fluid
flux density, tt is a first-order degradation coefficient and y is a zero-orderproduction coefficient. Other symbols are the same as before. Most parameter
estimation studies of solute t r a ns por t thus far have assumed steady-state flow
or some equivalent steady-state system using a transformation t hat avoids thesolution of the t r a ns i ent flow equation (eqn. 24) along with eqn. (32). For
steady-state flow in a homogeneous profile, the transport equation reduces to:~ + -~- = D~Zx~ - VTx - pc, + y(33)

where v = q/O is the mean pore water velocity.
Application of par a m et er estimation techniques to the determination ofv a r i o u s u n k n o w n coefficients in eqn. (33) are discussed below. For convenience,a distinction will be made between t r a ns port models involving local equilibrium and non-equilibrium conditions. Spatial variability of medium properties often limits the applicability of deterministic approaches based on eqn. (33)
to mass tr an s p o rt predictions at the field-scale. Thus, a stochastic formulationt h at explicitly considers the effects of areal variations in hydraulic fluxes onfield-scale solute transport will also be discussed.
Equilibrium transport

Consider the special case where adsorption is described by a linear isothermof the form s = kcr, where k is an empirical distribution coefficient. Equation(33) reduces then to:
R ~Cr _
0t

D 02cr~Cr~ x2 - v~--~x - pcr + 7(34)

where R = 1 + pk/O is the r e t a r d a t i o n factor. Inspection of eqn. (34) showst h at the equilibrium transport model contains five unknow n parameters: R, D,v, p and 7. Of these five parameters, a maximum of four can be fitted simultaneously to a set of observed cr(x, t) data. This dependency of at least one of
one of the coefficients follows immediately by noting t hat division of eqn. (34)by a constant permits one of the coefficients to be eliminated. Consequently, atleast one of the non-zero coefficients R, D, v, g or ~ must be known independently. Because values of D, p and 7 are not easily measured independently, inpractice either v or R (or both) must be known beforehand. For non-adsorbing
chemicals, R = 1 and v can be fitted to the data. This may be opportune if, forexample, u n c e r t a i n t y exists about the hydraulic flux or the effective waterc o n ten t in the system. For adsorbing chemicals, R can at least in principle beestimated directly using batch equilibration techniques.

277
P r o p e r s t i p u l a t i o n of b o u n d a r y c o n d i t i o n s for the s o l u t i o n of eqn. (34) can be
critical to the p a r a m e t e r e s t i m a t i o n problem. Choice of b o u n d a r y conditions
depends on the m e t h o d of solute i n t r o d u c t i o n into the system and on the modeof o b s e r v a t i o n of c o n c e n t r a t i o n s in time a n d / o r space. If solute of c o n c e n t r a tion Co is i n j e c t e d at a c o n s t a n t flow rate, t h e n the a p p r o p r i a t e b o u n d a r y
c o n d i t i o n is:( V C r - - D ~Cr~
= IVCo,~X ]x=o+(0,0 <. t <~ tot > to(35)

w h e r e to is the d u r a t i o n of pulse Co. T h e s o l u t i o n of eqn. (34) subject to eqn. (35)and to the a d d i t i o n a l conditions:cr = cix >~ 0,

t = 0

(36a)

~Cr/~X = 0

X --* ~ ,

t > 0

(36b)

m a y be used to describe r e s i d e n t c o n c e n t r a t i o n s in effectively semi-infinite
systems or in finite systems of l e n g t h L w h e n P = vL/D ~> 5 (Van G e n u c h t e n
and P a r k e r , 1984). S u b j e c t to some a d d i t i o n a l constraints, the s o l u t i o n of eqns.
(34)-(36) m a y suffice for smaller P in a finite system w h e n D is due to m o b i l e immobile zone i n t e r a c t i o n s (Parker, 1984; P a r k e r and Valocchi, 1986). Detection of r e s i d e n t solutions m a y be o b t a i n e d by sectioning soil cores and extracting the s o l u t i o n via a s u c t i o n a p p a r a t u s or, for non-adsorbed species, byr e p e a t e d w a s h i n g with tracer-free solution.
A n o t h e r mode of solute detection, p e r t i n e n t to c e r t a i n e x p e r i m e n t a l conditions, involves flux c o n c e n t r a t i o n s which r e p r e s e n t h y d r a u l i c flux-weighted
m e a n c o n c e n t r a t i o n s on a cross section of the porous medium. Such c o n c e n t r a tions arise n a t u r a l l y from o b s e r v a t i o n s of effluent from a soil column, alysimeter, or an e x t r a c t i o n well. M a c r o s c o p i c a l l y , the flux c o n c e n t r a t i o n cf m a y
be defined as a r a t i o of solute flux density to h y d r a u l i c flux density which
i n d i c a t e s (Kreft and Zuber, 1978):D ~c~cf = c r - - - - v ~x

(37)

Using eq. (37) we find t h a t eqns. (34) and (36) t r a n s f o r m to expressions ofi d e n t i c a l form in cf while eqn. (35) becomes:

cf =

c0,

0 ~< t ~< to

(0,

t >

I

to

(38)

If b o u n d a r y c o n d i t i o n s are i m p r o p e r l y formulated, appreciable e r r o r s inp a r a m e t e r estimates m a y arise for the inverse problem. As an example, consider a problem discussed by P a r k e r and V a n G e n u c h t e n (1984a). Resident
solutions were g e n e r a t e d using a n a n a l y t i c solution to eqns. (34)-(36) form e d i u m coefficients given in T a b l e 2. The i n i t i a l l y solute-free m e d i u m (ci = 0)
is subjected to a pulse i n p u t of Co = 1 0 0 # g c m -3 for a d u r a t i o n to = 5 days.

278
TABLE 2Comparison of fitted transport parameters using correct and erroneous detection modes for hypo.
thetical problem

True values for cr
Fit to cr solution
Fit to cf solution
Fit to c~ solution

D

R

#

100.00100.2885.6
161.7

2.50002.50062.52
3.00

0.25000.24960.34(0.25)

0.50000.49710.76(0.50)

D in cm2day 1, 7 in ttm cm Sday 1, values in parentheses was assumed known.

Resident solutions were predicted for 11 depths from 0-100 cm at t = 5 and 10
days. Parameters were estimated by a non-linear least-squares regression
procedure with equal weights for all observations and using only the t = 5days data. Initial guesses for all parameters were taken to be unity. Using thecorrect solution for resident concentrations, parameter estimates were foundto be very close to their true values whereas for the erroneous flux concentration boundary conditions, significant errors in parameter estimates occurred.
When # and y were fixed at their correct values, a two-parameter fit greatly
overestimated D and R. This overestimation is a characteristic occurrence
when resident concentrations are misinterpreted as flux concentrations (Parker and Van Genuchten, 1984b). The converse will be found when flux concentrations are misinterpreted as resident concentrations.
Non-equilibriumtransport

Although the linear equilibrium convection-dispersion model represented
by eqn. (34) is the most common approach to describing transport, many
circumstances arise for which it may be inadequate. Limitations arise due tokinetic control of adsorption reactions disregarded by eqn. (34). One approachto the extension of eqn. (31) has been to consider two types of adsorption sites,
one governed by equilibrium adsorption and one by reversible first-order kinetics. Such two-site models have been discussed by Selim et al. (1976), Cameron
and Klute (1977), Rao et al. (1979), De Camargo et al. (1979), Van Genuchten(1981), Flfihler and J u r y (1983) and others. If type-1 sites are regarded as
equilibrium-controlled, then R in eqn. (34) may be considered the retardation
due to adsorption on type-1 sites (R1). For kinetically-controlled type-2 sites,
the term:~C29--[ =~(k2cr -

c2)

(39)

must be added to the LHS of (34), where ~2 is the amount of solute adsorbed on
type-2 sites expressed as an equivalent solution concentration, ~ is a rate
coefficient, and k 2 is the slope of the equilibrium isotherm for type-2 sites

279
TABLE 3

Initial and fitted parameter values and associated residual SSQ's for the two-site model fit to boron
effluent data from Glendale clay loam

Initial values
Final values
Initial values
Final values

D

Rt

fl

~

1.0
216.0
2.00
47.7

1.00
3.58
10.00
4.30

0.500
0.564
0.200
0.600

0.200
14.2
0.200
0.424

SSQ

0.142
0.53

Units of D in cm 2 day -1, other values dimensionless.

d~2/dc r. A n a l y t i c s o l u t i o n s for the two-site model may be obtained for bothresident and flux d e t e c t i o n modes if a d s o r p t i o n / e x c h a n g e is linear (Coats andSmith, 1964; V a n G e n u c h t e n and Wierenga, 1976b; Parker and Van G e n u c h t e n ,1984a), while for non-linear adsorption numerical s o l u t i o n s may be employed(e.g., V a n G e n u c h t e n and Wierenga, 1976a; Fliihler and Jury, 1983).
As an application of the parameter estimation method to the two-site model,
consider an experimental b r e a k t h r o u g h curve for boron in a Glendale clay
loam described by V a n G e n u c h t e n (1974). A boron tracer of pulse duration
to = 5.06 days and h a v i n g a c o n c e n t r a t i o n co = 20 gg cm -~ was leached t h r o u g han initially solute-free 30 cm long column. The measured pore water velocity
was 38.5cm day -1. Four parameters were fitted to the data: the dispersion
coefficient D, total retardation factor Rt = 1 - R I + R2 where R2 = 1 + k2, asite distribution factor fl = R1/Rt,
and a reduced rate c o n s t a n t(o = ~(1 - fl)RtL/v where L is the c o l u m n l e n g t h (for more details see V a n
I.D

o

0.8

Z
0
I~

~c

oo

o

~ °~
O. 6

W

F ITTEO

Z
WO.
U
Z
C3

4

0.2

0.0

0

r

,

5

I0

PORE

15

20

VOLUMES

Fig. 4. Experimental boron effluent curve for Glendale clay loam and fitted curve using linear
two-site kinetic model (from Parker and Van Genuchten, 1984a).

280
Genuchten, 1981 or Parker and Van Genuchten, 1984a). Final parameter estimates and sums of squared residuals (SSQ) for two different sets of initial
parameter estimates using the solution for flux concentrations are shown in
Table 3. The results of the first example yielding the lowest SSQ are comparedwith the experimental data in Fig. 4. For the second set of initial parameter
estimates, the inversion program converged to a set of parameter values having
a larger SSQ using the stopping criteria that all parameters change by a factor
of less than 0.0005 from the previous iteration and allowing a maximum of 50
trials with no residual decrease at an iteration. In this instance a local minimum or plateau occurs in the response surface, which prevents convergence tothe global minimum. This points to the need for realistic initial estimates ofparameter values. Some specific guidelines for selecting initial values for thetwo-site model are discussed by Parker and Van Genuchten (1984a).
As pointed out by Van Genuchten (1981) and Nkedi-Kizza et al. (1983), amathematical problem of identical form to that of the two-site kinetic model
may be obtained when apparent non-equilibrium in the system is attributed todiffusional limitations between mobile and immobile pore regions, provided
that diffusion is approximated by an apparent first-order exchange process. In
non-dimensionalized form, the two-site and two-region models become identical, thus leading to difficulty in ascribing proper physical significance to fitted
model parameters (e.g., fl and o~ in the preceding example). Further complications relevant to the parameter identification problem arise owing to relationships between the first-order two-site/two-region model, and to interrelationships between the first-order two-site/two-region model and the equilibrium
convection-dispersion model on the simpler side, or a true diffusional kinetic
model on the more complex side. It may be demonstrated that under certain
conditions, the true diffusion kinetic analysis will degenerate to the form of thefirst-order kinetic model or to the local equilibrium model (Passioura, 1971;
Raats, 1981; Bolt, 1982; De Smedt and Wierenga, 1984; Van Genuchten, 1985;
Valocchi, 1985; Parker and Valocchi, 1986). In such circumstances, the parameter identification problem will be essentially indeterminant and model
selection should be governed by parsimony (e.g., Parker, 1984).
Further difficulties in the parameter estimation problem may arise due tonon-linearity and/or non-uniqueness in the adsorption isotherm, and hence inR (or R1 and R2), and to similar complications for the production and decaycoefficients g and y. A number of modelling approaches have been developed
which incorporate more complex equilibrium chemistry, ranging from the useof empirical non-linear isotherms with or without hysteresis (Van Genuchten
and Cleary, 1979) to more detailed geochemical analyses (Rubin and James,
1973; Jennings et al., 1982; Cederberg et al., 1985). The use of overly simplified
chemical models to describe complex systems can lead to difficulties when
predictions are extrapolated to situations that are different from those used in
the calibration process. For example, Jardine et al. (1985) found that a linear
two-site non-equilibrium transport model closely described observed A1
breakthrough curves for a soil in which Ca-A1 exchange and A1 polymerization

281
TABLE 4
Effect of influentiAl concentration 0~g ml 1) on calculated and fitted parameters for a linear
two-site kinetic model (after Jardine et al., 1985)
Influent
concentration
0.73
1.51
2.37
4.85
7.75

R

fl

fitted~(h 1)

calc.

fittedcalc.

fitted

40.2
25.2
16.8
13.6
*

51.8
42.2
49.6
17.1
12.0

0.632
0.415
0.216
0.325
0.329

0.620
0.417
0.261
0.363
0.385

0.132
0.124
0.127
0.262
0.330

* Data not available from isotherm.
c o n t r o l l e d t r a n s p o r t . An a c c e p t a b l e fit was o b t a i n e d in spite of highly nonl i n e a r a d s o r p t i o n - p r e c i p i t a t i o n relations. F o r experiments r u n at different
influent A1 c o n c e n t r a t i o n s , values of D, Rt, fl and ~ for the l i n e a r two-site model
were e s t i m a t e d from b r e a k t h r o u g h c u r v e s using a least:squares r e g r e s s i o nanalysis. I n d e p e n d e n t estimates of linearized d i s t r i b u t i o n coefficients, allowing
c a l c u l a t i o n of Rt and fl, were made from b a t c h equilibrium data. Comparisonsof fitted and c a l c u l a t e d v a l u e s i n d i c a t e d good c o r r e s p o n d e n c e for fl, while R t
v a l u e s were less c o n s i s t e n t (Table 4). N o t a b l e v a r i a t i o n s in the a p p a r e n t r a t e
c o n s t a n t ~ o c c u r r e d with influent c o n c e n t r a t i o n , w h i c h should not arise if ther e a c t i o n s t r i c t l y followed first-order chemical kinetics or was diffusionallycontrolled. H i g h e r - o r d e r chemical kinetics m a y be suggested. An equally plausible e x p l a n a t i o n is t h a t at lower influent Al-levels, i n c r e a s i n g n o n - l i n e a r i t y in
the i s o t h e r m produces g r e a t e r tailing in the b r e a k t h r o u g h curves. This tailingis t h e n s i m u l a t e d by the model as a k i n e t i c effect yielding lower fitted ~.
In m o r e r e c e n t studies, potassium t r a n s p o r t subject to n o n - l i n e a r K - C a
e x c h a n g e in a h e t e r o g e n e o u s medium with equilibrium and k i n e t i c a d s o r p t i o nsites has b e e n simulated. To e v a l u a t e p a r a m e t e r i n t e r a c t i o n s , the l i n e a r twosite k i n e t i c model was fitted to a series of h y p o t h e t i c a l b r e a k t h r o u g h curvest h a t were g e n e r a t e d n u m e r i c a l l y with a similar k i n e t i c model t h a t also acc o u n t e d for n o n - l i n e a r adsorption. While k e e p i n g the r a t e coefficients ~ fixed at
the o r i g i n a l values used in the n o n - l i n e a r model, the fitted (linear) r e t a r d a t i o nfactors were found to d e c r e a s e significantly with i n c r e a s i n g ~. This i n t e r a c t i o no c c u r s b e c a u s e the l i n e a r model will simulate tailing, due a c t u a l l y to nonl i n e a r i t y in the isotherms, t h r o u g h a d j u s t m e n t s in the r e l a t i v e p r o p o r t i o n ofe q u i l i b r i u m v e r s u s non-equilibrium a d s o r p t i o n sites, as well as t h r o u g h
c h a n g e s in the dispersion coefficient. I n t e r a c t i o n s b e t w e e n D and observedc h e m i c a l - k i n e t i c effects were also n o t e d by V a n G e n u c h t e n (1981) for nonl i n e a r b o r o n d i s p l a c e m e n t t h r o u g h a clay loam soil. Thus, w h e n a s t r o n g l yn o n - l i n e a r d i s p l a c e m e n t model is linearized or w h e n i m p o r t a n t k i n e t i c mechanisms are n e g l e c t e d in the p a r a m e t e r e s t i m a t i o n analysis, the results m a y282

become physically meaningless, even for parameters t h a t are not directly
affected by the approximation.

Field-scale transport
The question of field-scale variability effects on both the formulation of thedirect problem and on the inverse problem solution will now be considered. To
evaluate the effects of field-scale heterogeneities, a stochastic approach must
generally be taken. Here, discussion will be restricted to a relatively simple
model discussed by Parker and Van Genuchten (1984a) which is similar to the
one-dimensional stochastic transport models of Bresler and Dagan (1981),
Amoozegar-Fard et al. (1982) and Simmons (1982). Conceptually, the transport
region is regarded as being composed of numerous independent parallel soilcolumns (denoted as the "local" scale), each having specific properties andbeing subject to specific local boundary conditions. Additionally, it is assumedt h a t transport within each column can be described by the one-dimensional
convection-dispersion eqn. (34) with constant coefficients. Lateral flow, transverse dispersion, and vertical inhomogeneities are conveniently ignored.
To approximate transient flow, hydraulic fluxes and water contents areaveraged over the time and space domains of interest. At the local scale, mean
local surface hydraulic fluxes q0 and water contents 0 are defined as (Parker
and Van Genuchten, 1984a):tm

qo = ~

qo(t)dttm

Xm

xmtm f f O(x, t)dxdt
0

(40)

(41)

0

where (0, tin) and (0, xm) are the time and distance intervals for averaging andqo is the hydraulic flux at the soil surface. The time-averaged local velocity isthen taken as:= 0o10(42)

At the global- or field-scale, means over the areal domain A of the timeaveraged velocity are defined as:@} = ~

~dA

(43)

(A)

of the time-averaged surface flux qo:
1
(0o)

(/odA

=
(n)

and of the instantaneous surface flux qo at time t as:

(44)

283
(qo)(t)= ~(/o(t)dA(45)(A)
The equivalent steady-state time variable is then defined by the transformation:tt*(t) = | (qo)(z)dz/(~l)(46)
0Variations in local pore water velocities 0 are assumed to be lognormally
distributed with probability density function p(O) = N(#ln, a~n) where #~nand a~
are the mean and variance of In 0, respectively. The local dispersion coefficientis assumed to be perfectly correlated to 0 such that D = eO for deterministic
apparent dispersivity e. Ignoring R, # and y leads then to a three-parameterfield-scale transport model described by e, aln and #in.
Field-scale resident concentrations Cr are given as areal means over domain
A as:
Cr(X, t*)t*)dA
= (~) Cr (X'
(~)dA

(47)

where Cr is the local-scale concentration. Since fi is the only random variable
on A, p(~)dfi may be substituted for dA:
t*; ~)p(~)d~Cr(X, t*) = ~ Cr(X,
oo

(48)

I p(~)dfi
Ignoring variations in 0 in the field as small compared to q, field averaged flux
concentrations are obtained as:~f(x, t*) =I Ce(X, t*; ~)00(5)dO

i

(49)

00(0)d0

The integrals in eqns. (48) or (49) may be evaluated by numerical quadrature
once local boundary conditions are stipulated. As appropriate, these may be
defined by invoking a constant pulse duration to for the field or by imposing a
constant local mass loading which leads to local values of to varying inverselywith ~. For purposes of model calibration, note that for n random measurements of local concentration cr or Cr in the field, the field-average values ~ or
8f may be calculated by replacing the integrals in eqns. (48) or (49) by sums withweighting coefficients p(0) = 1/n assigned to each observation. To obtain ~f,local velocities must be known.
Application of the above model to a field experiment described by J u r y et al.

284
TABLE 5
P a r a m e t e r e s t i m a t e s a n d t h e i r s t a n d a r d e r r o r s o f e s t i m a t i o n ( ± ) for t h e field-scale b r o m i d e t r a c e r
e x p e r i m e n t of J u r y et al. (1982), v a l u e s fixed o n i n p u t a r e s h o w n in p a r e n t h e s e s
Trial

e
(mm)

(v)
( m m d a y - 1)

aln

SSQ

1
2
3
4
5
6

1.0 ± 227
(0.01)(1.0)(10.0)(100.0)123 ± 2230.5 ± 21.0
30.5 ± 1.8
30.5 + 1.8
29.7 __+ 1.7
24.7 + 1.6
23.6 ± 1.4

0.800 __+ 0.943
0.803 + 0.060
0.800 + 0.060
0.763 ± 0.063
0.373 ± 0.143(0.0)

0.0005177
0.0005185
0.0005177
0.0005182
0.000560.00061(1982) involving measurements of bromide transport in a 0.64 ha field undertransient hydraulic conditions is considered. Concentrations were determinedwith solution samplers at various depths taken at different times after additionof the tracer. Areal-averaged resident concentrations at a depth of 30 cm versus
transformed time t* were used as input to the inverse problem to estimate E,(note t h a t (~) is a known function of aln and #in). The results indicate large
uncertainty in ~ (Table 5). Using a two-parameter fit with e fixed at various
values, little sensitivity of the model to ~ for e ~< 10 mm was found. For largere, the mean error increases slightly, but compensation between s and a~ isobserved as the fitted values of a~, decrease for larger e. Thus, field-scaledispersion is reduced as local-scale dispersion increases. Note from Table 5 thateven when a~, = 0, which corresponds to the deterministic convection-dispersion model, SSQ is only moderately larger than for the three-parameterstochastic model.
Fitted concentration-time curves using trial 1 and trial 6 parameter estimates from Table 5, corresponding to the three-parameter stochastic model andthe two-parameter deterministic model, are compared to observed data for
x = 30cm in Fig. 5a. Using the same parameter estimates, breakthroughcurves were predicted for 60-90 cm depths and are compared to observed data
in Fig. 5b and c. It is apparent that while the two models are virtually indistinguishable at the calibration depth, they yield increasingly divergent predictions at greater depths. J u r y and Sposito (1985) presented a detailed comparison of the deterministic convection~iispersion model (denoted as the CDEmodel) and a "transfer function model" (TFM, Jury, 1982) which corresponds
to the stochastic model described above with e -* 0 (local piston flow). Resultsof their analyses for calibration of the two models to breakthrough curves foreach depth are shown in Fig. 6. Error ellipses computed by partitioning squared
deviations between model prediction s and observations by the X2 theorem
indicate better clustering for the TFM than the CDE model, but both models
are clearly too simplistic to fully describe the observed transient transport ofbromide. While the CDE tends to underpredict the rate of spreading of thesolute distribution due to the assumption of homogeneity, the TFM tends to

285
I

0.08

I

I

_

x = 30 cm
0.06

0.04
(..)

<~

~

o.o2

Observed~%/

itted ( I )

z

0
I-'<[
re
I-.
Z
IJJ
U
Z
0

0.00
0.06

U

0.04

0

0.02

hi
>
-I-<%
.J
hi

0.00
0.06

U

I

/'~'\
Observedg"-n¢~

//
,Y'%

i i/

V

%/

I
x =60cm

Pr.O,c,.o,,,\\,',.ic,ed

J.
I
_

I

I

/~ObservedI
x = 90cm _0.04~.~.~ \ \ . \/I/~.~

j Predicted ( , )'K..Predicted (2)

0.02

0.00
0

J//'/ I

I

I~'/''':='~'-

20

40

60

80

TIME

(t*),

'00

days

Fig. 5. Experimental areally-averaged resident concentrations at three depths as a function oftransformed time t* for bromide transport in a field soil (data from Jury et al., 1982). (a) data for
30 cm depth and fitted curve using the deterministic convection~lispersion model and a threeparameter stochastic model; (b) data for 60 cm depth and predicted curves for the two models usingparameters fitted to the 30 cm depth data; (c) same as (b) for the 90 cm depth data.

overpredict spreading due to the assumption of zero lateral transport and of
perfect correlation of velocity distributions with depth. These results
emphasize the caution which must be exercised when a calibrated model is used
to extrapolate to conditions and particularly to travel distances which differ
greatly from those employed for the inverse problem solution.
Data specification and data error
Given a suitable parametric model, the most critical factor for successful
parameter estimation analysis becomes specification of a data set which willlead to a unique and well-behaved solution of the inverse problem. In practice,experimental design will be a compromise between optimum conditions for the

286
0.015

I

0.010

!

I

I

I

I

I

I

'... ",,.,,.
\.'..,

0.005

A.

o

.,.,...'

I

i

0.01

0.02

I

V(md

I .00

--

0.:50

-

I

I

0.06

0.07

-I )

I

I
_I#*~

cm
6 0 cm30---0.67

I
0.05

Depth

-0

c

I

0.0:5 0 . 0 4

/

/%

//

o

......... 9 0 cm

, (oil/+..... 1 2 0 c m

t~..~

b
( P = .05 )

.
I

.

./

B.
i
I

I
2

Fig. 6. Mean and 0.05 probability error elipses for parameters fitted to breakthrough curves atdifferent depths for (a) CDE model, and (b) TFM model (fromJury and Sposito, 1985).
inverse problem and technological and economic constraints. In other words,
we would like to obtain a well-behaved inverse problem but not pay too dearlyin computational and experimental costs. The assumed model will dictate tosome extent the type of data which must be collected. Data sufficient to define
the hydrologic conditions will be necessary, although the requirements may
range from stipulation of (quasi-)steady-state hydraulic fluxes or water balances for simple hydrologic models, to complete specification of soil hydraulicproperties and flow boundary conditions for explicit modelling of transientvariably-saturated flow. Irrespective of the flow model, observations of tracer
concentrations must be made at various locations and/or times. Frequency andspatial distribution of observations will impact the feasibility of the inverse
problem.
As an example of the effects of data set composition on the inverse problem,
consider a hypothetical data set generated with the linear two-site/two-region
model consisting of 12 concentration observations versus time at two depths
each. Such a situation is not atypical of field experiments where the data may
be obtained from suction solution samplers (Parker and Van Genuchten,

287
TABLE 6
S t a t i s t i c s for i n v e r s e p r o b l e m s on 240 r e a l i z a t i o n s of a d a t a set s u b j e c t to r a n d o m e r r o r (afterW a g n e r a n d Gorelick, 1986)

True values

D

v

tt

7

100.0

25.00

0.2500

0.500

100.9
0.9
7.6

24.98
0.02
0.29

0.24950.00050.0164

0.487
0.013
0.227

101.7
1.7
12.4

25.06
0.06
0.73

0.24710.00290.0258

0.389
0.111
0.941

95.3
4.7
27.9

25.15
0.15
0.44

0.25990.00990.0635

0.901
0.401
1.926

Time and distance data
Mean
Bias
Standard deviation

Distance-only data
Mean
Bias
Standard deviation

Time-only data
Mean
Bias
Standard deviation

D in cm 2 d a y -1, v in cm d a y -1, # in d a y -1, 7 in p g c m -z d a y 1.

1984a). Using data for both depths, a five-parameter inversion for v, D, R, fl,and to successfully converged to the true values. Using concentration data for
only one depth, convergence could not be obtained for the five-parameter
inversion. Convergence for the smaller data set was feasible if to was assumed
known and only four parameters were estimated.
In actuality, exact correspondence between model predictions and observed
data will never occur, due in part to simplifications inherent in the parametricmodel as well as to measurement errors. Wagner and Gorelick (1986) studied
the effects of the latter via a Monte Carlo approach. Twenty concentration
observations were generated for two times each at 25 depths using the simple
one-dimensional convection-dispersion model for the same problem as discussed previously in conjunction with Table 2. Each concentration value wasassumed to be subject to random measurement error such t h a t the probabilityof a given observed value is normally distributed with standard deviation equal
to 20% of the true concentration. Parameter estimation analyses for 240 realizations of the data, with random error overlain, were carried out to determine
D, v, tt and 7. Means, variances and biases of the parameter estimates from the
240 replications were computed using a weighted least-squares procedure withweights inversely proportional to the estimated concentration values. The
results (Table 6) indicate t h a t parameter estimates are essentially unbiased.
This may be attributed in part to the adopted weighting scheme which corresponded to the assumed model for distributions of measurement error. Sensitivity to data error varied markedly with parameters. For v, the standard deviation of estimates for the 240 realizations was only 1% of the mean value, for288
D and ~ they were approximately 7%, and for p the standard deviation was 45%of the mean. The difficulty in accurately determining # would be evident inindividual inversion analyses from a large standard error of estimation in theparameter.
To compare the suitability of concentration observations in time or space,similar Monte Carlo analyses were performed using either observations in time
at single points in space or observations versus distance at single times. The
results (Table 6) indicate less satisfactory parameter estimates than when bothtime and distance data were used. Use of spatially distributed observations led
to results with markedly lower biases and standard errors than the use of
temporally distributed observations. Additionally, better agreements between
parameter confidence limits, obtained from individual inversion analyses andthose estimated from the Monte Carlo analyses, were obtained using spatiallydistributed data, indicating that the linear error analysis is more suitable forspatial data.
Similar preference of spatially over temporally distributed data appear to be
indicated also by analyses of J u r y and Sposito (1985) for field-scale tracer
movement. Parameter estimation variances for their transfer function model
calibrated using soil core data (spatial distributions) were lower than whenmodel calibration was carried out using suction solution sampler data (temporal distributions) at single depths. The interpretation of these results is,however, confounded by differences in experimental procedures. In the soilcore experiments, 36 replicates were available. The solute was moved into theground by daily sprinkler irrigation, and the entire experiment took only two
weeks. In contrast, the solution sampler experiment was irrigated by rainfallwith long periods of evaporation in between and had only 14 sites for solution
sampler analysis (W.A. Jury, pers. commun., 1986).
The results above indicate that the form and accuracy of data for the inverse
problem can incur varying degrees of estimation error in the parameters andin some instances impede convergence. Data uncertainty may affect the sensitivity of the results to different parameter estimation methods as well. In
weighted least-squares regression, the weighting coefficients may be very significant in this respect. In the parameter estimation studies of J u r y and Sposito(1985), three different objective functions were studied based on an ordinary
least-squares procedure, a least-squares minimization of observed and predicted time moments to fourth-order, and a maximum likelihood analysis. Significant differences in the parameter estimates by the three methods were observedwith the moments method giving greater weight than the other methods to
low-concentration observations at longer times and the ordinary least-squares
procedure giving them the least weight.
FUTURE DIRECTIONS
As unsaturated zone flow and transport simulations become increasingly
common due to heightened concerns with the prevention and remediation of

289
groundwater contamination, the need for efficient and accurate methods ofmodel calibration will become more crucial. Also, because of inherent limitations of any modelling effort, it is important that probable errors associatedwith model parameters and ultimately with model predictions, be quantifiable.
Parameter estimation methods offer the most suitable means of meeting theserequirements. Such methods have been widely applied in surface hydrology
and lately also in saturated subsurface flow and transport modelling; they have
only recently begun to be applied to vadose zone flow and transport model
calibration. In this paper, we have summarized the basic methodology involvedin parameter estimation and reviewed recent applications pertinent to vadosezone flow and transport modelling. It is clear that while certain progress has
been made, a great deal of work remains to be done.
To date, the majority of studies have involved laboratory-scale experiments
and rather simple parametric models. The few reported field-scale parameter
estimation studies have involved simple model formulations, boundary conditions and/or soil conditions. Efforts are needed to extend parameter estimation methods to more complex field conditions using models capable of accommodating soil heterogeneity, variable and uncertain boundary conditions,
simultaneous flow and transport, complex biochemical processes and otherphenomena. Parameter identification and parameter estimation studies willneed to be carried out in conjunction to meet these needs. Parametric formulations for the direct problem incorporating different processes in varying degrees of sophistication will need to be studied to determine the most efficient
representations which do not excessively compromise accuracy and precision.
For example, to accommodate effects of soil heterogeneity, investigations must
address the development of concise parametric representations of soil variability, perhaps employing scaling analyses in conjunction with geostatisticalmethods and incorporating prior information to facilitate blocking of statistically homogeneous zones.
Given that a variety of approaches to the direct problem will be possible,
explicit methods of dealing with the parameter identification problem will need
to be addressed. To do this, accuracy and precision associated with different
parametric representations must be considered simultaneously to evaluate thepoint at which increased complexity and accuracy of the model is counteracted
by loss in precision due to an inability to estimate a larger number of parameters with confidence.
The analysis of optimum experimental design and optimum input data requirements for unsaturated flow and transport model inversion via parameter
estimation methods has not been investigated in a systematic manner forfield-scale problems. Attention must be given to the selection of the mostsuitable measurements (e.g., water contents, heads, solution concentrations,solid phase concentrations), the most appropriate methods of measurement(e.g., concentrations via suction samples or cores), and optimum distributionsof observations in time and space. Consideration of optimum initial and boundary conditions for experiments when these are wholly or partially controlled,
suitability of ambient environmental conditions, and sensitivity to uncertain-

290
ties in b o u n d a r y and initial c o n d i t i o n s will be i m p o r t a n t in the d e v e l o p m e n t offield c a l i b r a t i o n methods.
Finally, the i m p o r t a n c e of i m p r o v i n g the efficiency and r o b u s t n e s s ofn u m e r i c a l m e t h o d s for b o t h t h e d i r e c t problem and the inverse problem isemphasized. B e c a u s e p a r a m e t e r e s t i m a t i o n involves r e p e a t e d s o l u t i o n of thed i r e c t problem, efficiency of the l a t t e r becomes i n c r e a s i n g l y crucial to thep r a c t i c a l i t y of p a r a m e t e r e s t i m a t i o n especially w h e n spatial v a r i a b i l i t y effects
are explicitly c o n s i d e r e d in the direct problem and d i m e n s i o n a l i t y and modelc o m p l e x i t y increase. B e c a u s e p a r a m e t e r values are u n d e r a u t o m a t i c c o n t r o l in
the i n v e r s i o n analysis, r o b u s t n e s s of the d i r e c t problem solution, in the senseof i n s e n s i t i v i t y of model a c c u r a c y to p a r a m e t e r values, will be of f u n d a m e n t a limportance. A l t h o u g h m u c h r e s e a r c h on optimization m e t h o d s has been done,
this is still an a c t i v e field and gains m a y be a n t i c i p a t e d from improved
efficiency and r o b u s t n e s s of optimization techniques. The d e v e l o p m e n t anda p p l i c a t i o n of adjoint s t a t e m e t h o d s to e v a l u a t e p a r a m e t e r g r a d i e n t vectors
would a p p e a r to be a p a r t i c u l a r l y fruitful endeavor, w h i c h should yield subs t a n t i a l benefits especially w h e n a large n u m b e r of p a r a m e t e r s must be estimated. B e c a u s e the c h o i c e of objective f u n c t i o n f o r m u l a t i o n can h a v e signific a n t effects on the c o m p u t e d results, p a r t i c u l a r l y w h e n d a t a u n c e r t a i n t y islarge, the objective f u n c t i o n selection must be given suitable a t t e n t i o n . F o rconclusion, m u c h w o r k r e m a i n s to be done in the d e v e l o p m e n t of p a r a m e t e r
e s t i m a t i o n m e t h o d s for u n s a t u r a t e d zone flow and t r a n s p o r t models, but as thequest for field-scale p r e d i c t i v e capability and the need for efficient c a l i b r a t i o nbecomes i n c r e a s i n g l y great, these m e t h o d s will be of ever i n c r e a s i n g utility.
ACKNOWLEDGEMENTS
The a u t h o r s wish to t h a n k Dr. C.A. Beattie, D e p a r t m e n t of M a t h e m a t i c s ,VPI & SU, for his r e v i e w and c o n s t r u c t i v e c o m m e n t s on the section c o n c e r n i n g
p a r a m e t e r e s t i m a t i o n methods.
REFERENCES
Abeele, W.V., 1984. Hydraulic testing of crushed Bandelier Tuff. Rep. LA-10037-MS,Los AlamosNat. Lab., 21 pp.
Amoozegar-Fard, A., Nielsen, D.R. and Warrick, A.W., 1982. Soil solute concentration distributions for spatially varying pore water velocities and apparent diffusion coefficients. Soil Sci.
Soc. Am., J., 46: 3-9.
Bard, Y., 1974. Nonlinear Parameter Estimation. Academic Press, New York, N.Y., 341 pp.
Beck, J.V. and Arnold, K.J., 1977. Parameter Estimation in Engineering and Science. Wiley, New
York, N.Y., 501pp.
Bolt, G.H., 1982. Movement of solutes in soil: principles of adsorption/exchange chromatography.
In: G.H. Bolt (Editor), Soil Chemistry: B. Physicochemical Models. Elsevier, Amsterdam, 2nded., pp. 28~348.
Bouwer, H. and Jackson, R.D., 1974. Determining soil properties. In: J. van Schilfgaarde (Editor),
Drainage for Agriculture. Am. Soc. Agron., Madison, Wisc., Monogr., 17: 611-672.
Bresler, E. and Dagan, G., 1981. Convective and pore scale dispersive solute transport in unsaturated heterogeneous fields. Water Resour. Res., 17: 1683-1693.

291
Brooks, R.H. and Corey, A.T., 1964. Hydraulic properties of porous media. Civ. Eng. Dep., Colorado
State Univ., Fort Collins, Colo., Hydrol. Pap. No. 3.
Cameron, D.A. and Klute, A., 1977. Convective~ispersive solute transport with a combined
equilibrium and kinetic adsorption model. Water Resour. Res., 13: 183-188.
Carrera, J., 1984. Estimation of aquifer parameters under transient and steady-state conditions.
Ph.D. Dissertation, Dep. Hydrol. Water Resour., Univ. of Arizona, Tucson, Ariz., 257 pp.
Carrera, J. and Neuman, S.P., 1984. Adjoint state finite element estimation of aquifer parameters
under steady-state and transient conditions. Proc. 5th Int. Conf. on Finite Elements in Water
Resour., Springer, New York, N.Y., pp. 199-211.
Carrera, J. and Neuman, S.P., 1986. Estimation of aquifer parameters under transient and steadystate conditions: 2. Uniqueness, stability, and solution algorithms. Water Resour. Res., 22:
211-227.
Cederberg, G.A., Street, R.L. and Leckie, J.O., 1985. A groundwater mass transport and equilibrium
chemistry model for multicomponent systems. Water Resour. Res., 21: 1095-1104.
Coats, K.H. and Smith, B.D., 1964. Dead-end pore volume and dispersion in porous media. Soc. Pet.
Eng., J., 4: 73-84.
Cooley, R.L., 1982. Incorporation of prior information on parameters into nonlinear regression
groundwater flow models. 1. Theory. Water Resour. Res., 18: 965-976.
Cooley, R.L., 1985. A comparison of several methods of solving nonlinear regression groundwaterflow problems. Water Resour. Res., 21: 1525-1538.
Dane, J.H. and Hruska, S., 1983. In situ determination of soil hydraulic properties during drainage.
Soil Sci. Soc. Am., J., 47: 61~624.
De Camargo, O.A., Biggar, J.W. and Nielsen, D.R., 1979. Transport of inorganic phosphorus in an
Alfisol. Soil Sci. Soc. Am., J., 43: 884-890.
Dennis, J.E. and Schnabel, R.B., 1983. Numerical methods for unconstrained optimization andnonlinear equations. Prentice-Hall, Englewood Cliffs, 378 pp.
Dennis, J.E., Gay, D.M. and Welsch, R.E., 1981. An adaptive nonlinear least-squares algorithm.
Assoc. Comput. Mach., Trans. Math. Software, 7: 348-383.
De Smedt, F. and Wierenga, P.J., 1979. A generalized solution for solute flow in soils withmobile-immobile water. Water Resour. Res., 15: 1137-1141.
De Smedt, F. and Wierenga, P.J., 1984. Solute transfer through columns of glass beads. Water
Resour. Res., 20: 225-232.
Fliihler, H. and Jury, W.A., 1983. Estimating solute transport using nonlinear, rate dependent, twosite adsorption models. Swiss Fed. Inst. For. Res., Birmensdorf., Rep. no. 245.
Gill, P.E., Murray, W. and Wright, M.H., 1981. Practical Optimization. Academic Press, New York,
N.Y., 401 pp.
Hiebert, K.L. 1981. An evaluation of mathematical software that solves nonlinear least squares
problems. Assoc. Comput. Mach., Trans. Math. Software, 7: 1-16.
Hornung, U., 1983. Identification of nonlinear soil physical parameters from an input-output
experiment. In: P. Deuflhardt and E. Hairer (Editors), Workshop on Numerical Treatment of
Inverse Problems in Differential and Integral Equations. Birkhauser, Boston, Mass., pp.
227-237.
Isaacson, E. and Keller, H., 1966. Analysis of numerical methods. Wiley, New York, N.Y., 541 pp.
Jardine, P.M., Parker, J.C. and Zelazny, L.W., 1985. Kinetics and mechanisms of aluminum adsorption on kaolinite using a two-site nonequilibrium transport model. Soil Sci. Soc. Am., J., 49:
867472.
Jennings, A.A., Kirkner, D.J. and Theis, T.L., 1982. Multicomponent equilibrium chemistry ingroundwater quality models. Water Resour. Res., 18: 1089-1096.
Jennings, J.W. Jr., McGregor, D.S. and Morse, R.A., 1985. Simultaneous determination of capillarypressure and relative permeability by automatic history matching. 60th Annu. Tech. Conf., Las
Vegas, SPE 14418.
Jury, W.A., 1982. Simulation of solute transport with a transfer function model. Water Resour.
Res., 18: 363-368.

292
Jury, W.A. and Sposito, G., 1985. Field calibration and validation of solute transport models for theunsaturated zone. Soil Sci. Soc. Am., J., 49: 1331-1341.
Jury, W.A., Stolzy, L.H. and Shouse, P., 1982. A field test of a transfer function model for predicting
solute transport. Water Resour. Res., 18: 369-375.
Klute, A., 1972. The determination of the hydraulic conductivity and diffusivity of unsaturated
soils. Soil Sci., 113: 164-176.
Kool, J.B. and Parker, J.C., 1987. Development and evaluation of closed-form expressions forhysteretic soil hydraulic properties. Water Resour. Res., 23: 105-114.
Kool, J.B., Parker, J.C. and Van Genuchten, M.Th., 1985. Determining soil hydraulic properties
from one-step outflow experiments by parameter estimation. I. Theory and numerical studies.
Soil Sci. Soc. Am., J., 49: 1348-1354.
Kool, J.B., Parker, J.C. and Van Genuchten, M.Th., 1986. The inverse problem for hysteretic
unsaturated flow. Proc. 6th Int. Conf. on Finite Elements in Water Resour., Springer, New
York, N.Y., 337-346pp.
Kreft, A. and Zuber, A., 1978. On the physical meaning of the dispersion equation and its solution
for different initial boundary conditions. Chem. Eng. Sci., 33: 1471-1480.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. SIAM
J. Appl. Math., 11: 431~441.
Mor~, J.J., 1977. The Levenberg-Marquardt algorithm: Implementation and theory. In: G.A.
Watson (Editor), Lecture Notes in Mathematics 630. Springer, New York, N.Y., pp. 105-116.
Neuman, S.P., 1980. A statistical approach to the inverse problem of aquifer hydrology. 3. Improved
solution method and added perspective. Water Resour. Res., 16: 331-346.
Nkedi-Kizza, P., Biggar, J.W., Van Genuchten, M.Th., Wierenga, P.J., Selim, H.M., Davidson, J.M.
and Nielsen, D.R., 1983. Modeling tritium and chloride 36 transport through an aggregatedoxisol. Water Resour. Res., 19: 691-700.
Parker, J.C., 1984. Analysis of solute transport in column tracer studies. Soil Sci. Soc. Am., J., 48:
719-724.
Parker, J.C. and Van Genuchten, M.Th., 1984a. Determining transport parameters from laboratory
and field tracer experiments. Va. Agric. Exp. Stn., Bull., 84-3: 92pp.
Parker, J.C. and Van Genuchten, M.Th., 1984b. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour. Res., 20: 866-872.
Parker, J.C. and Valocchi, A.J., 1986. Constraints on the validity of equilibrium and first-orderkinetic transport models in structured soils. Water Resour. Res., 22: 399-409.
Parker, J.C., Kool, J.B. and Van Genuchten, M.Th., 1985. Determining soil hydraulic properties
from one-step outflow experiments by parameter estimation: II. Experimental studies. Soil Sci.
Soc. Am., J., 49: 1354-1360.
Passioura, J.B., 1971. Hydrodynamic dispersion in aggregated media, I. Theory. Soil Sci., 111:339-344.
Philip, J.R., 1955. Numerical analysis of equations of the diffusion type with diffusivity concentration dependent. Trans. Farraday Soc., 51: 885-892.
Rao, P.S.C., Davidson, J.M., Jessup, R.E. and Selim, H.M., 1979. Evaluation of conceptual models
for describing non-equilibrium and adsorption-desorption of pesticide during steady flow insoils. Soil Sci. Soc. Am., J., 43: 22-28.
Raats, P.A.C., 1981. Transport in structured porous media. Proc. Euromech. 143, Sept. 2-4, Delft,
The Netherlands, pp. 221 226.
Richie, E.B., Schnabel, R.R. and Hoover, J.R., 1987. Estimation of soil-water parameters fromunsaturated steady-state flow. Proc. Int. Conf. Infiltration Dev. Appl., Univ. of Hawaii, Honolulu, Hawaii.
Rubin, J. and James, R.V., 1973. Dispersion-affected transport of reacting solutes in saturated
porous media: Galerkin method applied to equilibrium-controlled exchange in unidirectional
steady water flow. Water Resour. Res., 9: 1332-1356.
Sadeghipour, J. and Yeh, W.W.-G., 1984. Parameter identification of groundwater aquifer models:A generalized least squares approach. Water Resour. Res., 20: 971-979.

293
Selim, H.M., Davidson, J.M. and Mansell, R.S., 1976. Evaluation of a two-site adsorption-desorption model for describing solute transport in soils. Proe. 1976 Summer Comput. Simulation Conf., Washington, D.C., pp. 444-448.
Simmons, C.S., 1982. A stochastic-convective transport representation of dispersion in onedimensional porous media systems. Water Resour. Res., 18: 1193-1214.
Sorooshian, S. and Dracup, J.A., 1980. Stochastic parameter estimation procedures for hydrologic
rainfall-runoff models: Correlated and heteroscedastic error cases. Water Resour. Res., 16:
430-442.
Townley, L.R., 1983. Numerical models of groundwater flow: Prediction and parameter estimationin the presence of uncertainty. Ph.D. Dissertation, Dep. Civ. Eng., Mass. Inst. Teeh., 354pp.
Townley, L.R. and Wilson, J.L., 1985. Computationally efficient algorithms for parameter estimation and uncertainty propagation in numerical models of groundwater flow. Water Resour. Res.,
21: 1851-1861.
Troutman, B.M., 1985. Errors and parameter estimation in precipitation-runoff modeling. 1.
Theory. Water Resour. Res., 21: 1195-1213.
Va]occhi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solutetransport through homogeneous soils. Water Resour. Res., 21: 808~820.
Van Genuchten, M.Th., 1974. Mass transfer studies in sorbing porous media. Ph.D. Thesis, New
Mexico State Univ., Las Cruces, N.M., 148pp.
Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils. Soil Sei. Am., J., 44: 892-989.
Van Genuchten, M. Th., 1981. Non-equilibrium transport parameters from miscible displacement
experiments. U.S. Salinity Lab. Res. Rep., no. i19, Riverside, Calif.
Van Genuchten, M.Th., 1985. A general approach for modeling solute transport in structured soils.
Proc. 17th Int. Congr., Int. Assoc. Hydrogeol., Mem., 17 (2): 513-526.
Van Genuchten, M.Th. and Wierenga, P.J., 1976a. Numerical solution for convective dispersion
with intra-aggregate diffusion and non-linear adsorption. In: C.G. Vansteenkiste (Editor),
System Simulation in Water Resources. North-Holland, Amsterdam, pp. 275-291.
Van Genuchten, M.Th. and Wierenga, P.J., 1976b. Mass transfer studies in sorbing porous media:
I. Analytical solutions. Soil Sci. Soc. Am., J., 40: 473-480.
Van Genuchten, M.Th. and Cleary, R.W., 1979. Movement of solutes in soil: computer-simulated
and laboratory results. In: G.H. Bolt (Editor), Soil chemistry, B. Physiochemieal Models.
Elsevier, Amsterdam, pp. 349-386.
Van Genuchten, M.Th. and Parker, J.C., 1984. Boundary conditions for displacement experiments
through short laboratory soil columns. Soil Sci. Soc. Am., J., 48: 703-708.
Van Genuchten, M.Th. and Nielsen, D.R., 1985. On describing and predicting the hydraulicproperties of unsaturated soils. Ann. Geophys., 3: 61~628.
Van Grinsven, J.J.M., Dirksen, C. and Bouten, W, 1985. Evaluation of the hot air method formeasuring soil water diffusivity. Soil Sci. Soc. Am., J., 49: 1093-1100.
Watson, K.K., 1966. An instantaneous profile method for determining the hydraulic conductivityof unsaturated porous materials. Water Resour. Res., 2: 709-715.
Wagner, B.J. and Goreliek, S.M., 1986. A statistical methodology for estimating transport parameters: Theory and applications to one-dimensional advective~lispersive systems. Water
Resour. Res., 22: 1303-1316.
Wittmeyer, G.W. and Neuman, S.P., 1985. Inverse solution of variably saturated flow using adjointstate finite elements. Am. Geophys. Union, Washington, D.C., Abstr. EOS Trans., 66 (96): p. 886.
Yeh, W.G., 1986. Review of parameter estimation procedures in groundwater hydrology: The
inverse problem. Water Resour. Res., 22: 95-I08.
Zachmann, D.W., Duchateau, P.C. and Klute, A., 1981. The calibration of the Richards flow
equation for a draining column by parameter identification. Soil Sci. Soc. Am., J., 45: I01~I015.
Zachmann, D.W., Duchateau, P.C. and Klute, A., 1982. Simultaneous approximation of watercapacity and soil hydraulic conductivity by parameter identification. Soil Sci., 134:157 163.