Inverse Optimization, Calibration and Validation of
Simulation Models at the Field Scale
J. Šimůnek1 & J. A. de Vos2
U.S. Salinity Laboratory, USDA-ARS, 450 West Big Springs Road, Riverside, CA 92507, USA
2
DLO Research Institute for Agrobiology and Soil Fertility, P.O. Box 14, 6700 AA
Wageningen, The Netherlands
1

Abstract
An overview is given of the issues of parameter estimation, model verification, and model
validation as applied to field-scale subsurface flow and transport problems. We briefly review
inverse optimization methods for estimating soil hydraulic parameters from a variety of field
experiments, including tension disc infiltrometry, cone penetrometry, and gravity drainage
experiments. An example is presented showing calibration of the numerical HYDRUS-2D modelusing data of a tile-drainage experiment. The hydraulic characteristics of the layered soil profileat the site were identified based on the joined use of laboratory data, field monitoring data, andthe numerical model. A split sampling technique was used to test applicability of the numericalmodel for this study.
Keywords: Simulation models, parameter estimation, calibration, tile-drained field

Introduction
Computer models based on numerical solutions of the flow and solute transport equations are
increasingly being used for a wide range of applications in research and management.
Application of the models is enhanced by the ever increasing power of personal computers,
and the development of more accurate and numerically stable numerical techniques. Precisionof the obtained predictions depends to a large extent upon the accuracy of the available model
input parameters, and on a successful description of the actual physical system, including soil
heterogeneity and the possible presence of nonequilibrium flow and transport conditions suchas preferential flow. Parameters in the soil hydraulic functions characterizing the water
retention and permeability properties are the most important input variables for models basedon numerical solutions of the variably-saturated flow (Richards) equation. These hydraulicparameters directly influence the mobility of various chemical species through their effect on
pore-water flow velocities and water contents. Transport and chemical parameters
additionally affect the rate of spreading of chemicals and their distribution between mobile
and immobile phases.
Water flow and solute transport models are increasingly being applied to natural
subsurface systems. While a major purpose of their application is to provide future
predictions, they are more often used also to interpret the complex interactions found inlaboratory and field data (Steefel & Van Cappellen, 1998). Models are also often used to
study the sensitivity of natural systems to selected variables and processes, and to determine
the uncertainties in predictions. The results of a sensitivity analysis can provide insight into
the relative importance of individual model parameters, and thus help guide the allocation of
resources for further laboratory and field investigations (Šimůnek & Valocchi, 2000).

431

There have been more efforts lately to test the reliability of model predictions, in part
motivated by the widespread use of numerical models for management purposes. It is still an
open question whether we ultimately can use water flow and solute transport models as true
predictive tools to make with confidence quantitative forecasts necessary for practical
purposes. At present, there is even controversy regarding the predictive capability of simpler
models of saturated water flow and single component transport (e.g., Konikow & Bredehoeft,
1992; de Marsily et al., 1992; Anderson & Woessner, 1992; Mazoszewski & Zuber, 1992;
Oreskes et al., 1994). Attempts to deal with the problem of reliability of model predictions has
spurred the introduction in the literature of a long list of terms that are often only loosely definedand frequently substituted for each other, such as model verification, validation, confirmation,and calibration. In this paper we give first brief definitions of the most relevant terms as used in
the above references, and provide a short summary of ongoing discussions on model validation.
Next we will review inverse optimization as used for estimating soil hydraulic parameters fromfield experiments, and conclude with a calibration of the numerical HYDRUS-2D model using
field data collected as part of a tile-drainage experiment.

Reliability of Model Predictions
A model is 'a representation of a real system or process' (Konikow & Bredehoeft, 1992). The
word model can represent many things, such as conceptual (or process-based) model,
mathematical model (deterministic or stochastic), numerical or analytical model, physicalmodel (analog), and many others. A conceptual model is a 'hypothesis for how a system orprocess operates' (Konikow & Bredehoeft, 1992) or 'a qualitative description of a system andits representation relevant to the intended use of the model' (OECD/NEA, 1990). A
mathematical model is a representation of a conceptual model with objects, forces and events
replaced with mathematical expressions. A majority of the mathematical models used today tosimulate water flow and solute transport in both saturated and variably saturated systems are
deterministic mathematical models that generally solve a set of partial differential equations
representing conservation of mass, momentum, and energy. Since most practical problems
cannot be simplified and/or idealized to the point so that analytical solutions exist, mostmodels use various numerical techniques to solve the governing equations subject toappropriate initial and boundary conditions. Implementation of a numerical algorithm thatsolves the governing equation(s) results in a computer code that can be considered a genericmodel. A generic model applied to a particular geographical area with given (optimized orcalibrated) parameters and initial conditions, represents a site-specific model (Konikow &Bredehoeft, 1992).
The term verification can refer to the governing equation, the code or the model(Anderson & Woessner, 1992). Verification of the governing equation involves demonstrationthat the equation in the model accurately describes the simulated process. While it is generally
accepted that Darcy's law combined with conservation of mass accurately describes saturated
groundwater flow, the applicability of the Richards equation to flow in the unsaturated zone
and the convection-dispersion equation to solute transport at the field scale is often questioned(Jury & Flhler, 1992). Verification of a numerical code consists of showing that the results
generated by the model for simpler problems are consistent with available analyticalsolutions, or are the same, or similar, as results generated with other numerical codes. The
latter procedure is often also called benchmarking. Verification of a code should ensure thatthe equations constituting the mathematical model are correctly encoded and solved.
However, verification of a code does not mean verification of the governing equation, nordoes it imply absolute certainty that the numerical solution is implemented correctly. Also,432

available analytical solutions are often limited to idealized transport domains, homogeneousand isotropic media, and uniform initial and constant boundary conditions. The very reason
for developing numerical models is to go beyond the range of available analytical solutions,
i.e., to allow irregular transport domains, nonhomogeneous and anisotropic media, variableboundary conditions, and nonlinear processes, i.e. to use them for situations or conditions forwhich they can not be verified. Verification in such conditions is often accomplished using
approximate tests of having internal consistency and accuracy, such as mass conservation,
global mass-balance errors, and sensitivity to changes in mesh size and time steps.
Model verification has often been used synonymously with model validation (Anderson& Woessner, 1992). We agree with the suggestions by Narasimhan (1987) to associate model
verification with the accuracy of the invoked numerical solution schemes and the coding of amodel, and model validation with the inherent capability (or the degree of validity) of a modelin describing a set of processes (in our case subsurface flow and transport processes). Many
definitions have been used for the term model validation; still, much controversy remains ofwhat it means and whether model validation can be truly achieved (Oreskes et al., 1994). For
example, the International Atomic Energy Agency defines a validated model as one that gives'a good representation of the actual processes occurring in a real system' (IAEA, 1982).
Similarly, the U.S. Department of Energy defines a validated model as one that 'reflects the
behavior of the real world' (US DOE, 1986). Another definition states that model 'validation
is a process of obtaining assurance that a model is a correct representation of the process or
system for which it is intended' (OECD/NEA, 1990; Mazoszewski & Zuber, 1992).
Regulatory agencies often request scientists and model developers to validate their models or
to show that the models were independently validated. This did lead to several international
projects for the purpose of addressing the issues of verification and validation ofhydrogeological models assessing nuclear waste repositories. These included projects
INTRACOIN (the International Nuclide Transport Code Intercomparison Study, SNPI, 1986),
INTRAVAL (an International Project to Study Validation of Geosphere/Transport Models,
SNPI, 1987), and HYDROCOIN (Hydrologic Code Intercomparison Study, OECD/NEA,
1990).
Model validation remains a contentious topic among soil and groundwater hydrologists
and scientists. For example, the journal Advances in Water Resources devoted in 1992 two
special issues to the topic of 'Validation of Geo-hydrological models' (Hassanizadeh &Carrera, 1992). Konikow & Bredehoeft (1992), in their enlightening contribution entitled'Ground-water models cannot be validated', suggest that since groundwater models areembodiment of scientific hypothesis they cannot be proven or validated, similarly as any
scientific hypothesis or theory, but only tested and invalidated. They also argue that 'the terms
validation and verification are misleading and their use in groundwater science should beabandoned in favor of more meaningful model-assessment descriptors'. In a comment on thelatter article, de Marsily et al. (1992) first disagreed with that notion, but later offered acompromise statement by stating that 'we do not validate our models, but we try to show thatthey are not invalidated by the data!', which does not seem to contradict the arguments ofKonikow and Bredehoeft. Similarly as Konikow & Bredehoeft (1992), Oreskes et al. (1994)
also argue that verification and validation of numerical models of natural systems is
impossible since such systems are never closed and model results are always non-unique.
Model calibration is generally defined as the process of tuning a model for a particular
problem or case by manipulating the input parameters (e.g., hydraulic conductivities), andinitial or boundary conditions, within reasonable ranges until the simulated model results
closely match the observed variables (e.g., pressure heads, water contents, fluxes). The
general approach in model calibration is to select some merit or objective function that is a433

measure of the agreement between measured and modeled data, and that is directly or
indirectly related to the adjustable parameters. The best-fit parameters are obtained byminimizing the objective function. Model calibration, i.e., minimization of the merit function,
can be achieved by trial and error or, as is becoming more popular, by using an automatedminimization or parameter estimation technique. A model is considered calibrated when it
reproduces data within some subjectively acceptable level of precision (Konikow &Bredehoeft, 1992). Model calibration is often also called history matching. When the ultimate
goal of model calibration is not merely to calibrate a model, but rather to optimize unknown
parameters in that model, the process is often called parameter optimization or parameter
estimation.
A two-step calibration scheme or split sampling has sometimes been used for the model
calibration process, or instead of model validation (if one accepts the notion that models
cannot be formally validated). The existing data set is then divided into two parts with thefirst (calibration) part being used to calibrate the model and to estimate all necessary
parameters, while the second (validation or verification) part of the data set serves to compare
predicted and measured data values using the parameters calibrated against the first part of thedata set. Although successful performance of the model during the second part does notconstitute a rigorous validation of the model, it is often accepted that the model is validated if
it predicts the system response during this second phase within acceptable limits (Anderson &Woessner, 1991). A two-step calibration process of this type is sometimes called historicalvalidation.
The most rigorous test of the performance of a model is the postaudit or predictive
validation (Anderson & Woessner, 1992). The calibrated model in this case is first used tomake a prediction of the future behavior of the system; only after the modeling study iscompleted are new data collected and compared against the model predictions. A literature
study by Anderson & Woessner (1992) revealed five modeling studies that underwent arigorous postaudit test. None of the models was found to accurately predict the future
behavior of the system. This was partly due to errors in the conceptual model of thehydrogeological system and partly due to failure to correctly predict future boundary
conditions. Anderson & Woessner (1992) concluded that the issue of model validation ismainly a regulatory one, not a scientific one. A model can never be proven valid from ascientific standpoint because our understanding of a system will always be incomplete.
Accepting the argument that models cannot be validated experimentally and can be verified
mathematically only in a limited sense, we now discuss the problems of parameter estimationand model calibration.

Parameter Estimation
Objective Function
The general approach in parameter estimation is to select a merit or objective function that isa measure of the agreement between measured and modeled data, and which is directly or
indirectly related to the adjustable parameters. The merit function is often based on a
maximum likelihood estimator, which simplifies to a weighted least-squares problem foruncorrelated measurement errors (Bard, 1974):n(  ) =  wi [qi* - qi (  ) ]2

(1)

i=1

434

where  = {1, 2,...,m} is the vector of optimized parameters, m is the number of optimized
parameters, q* = {q1*, q2*,..., qn*} is the vector of observations, q() = {q1, q2,..., qn} is thecorresponding vector of model predictions as a function of the unknown parameters being
optimized, n is the number of observations, and wi is the weight of a particular measured point.
The weighted least-squares estimator (1) is a maximum-likelihood estimator as long as theweights, wi, contain the measurement error information such that
wi =

1



2

=

i

1
variance of measurement error of q*i(2)

The robustness of the least-squares criterion for the estimation of model parameters has
recently been questioned by Finsterle & Najita (1998). They pointed out that the least-squares
criterion causes outliers to strongly influence the final values of optimized parameters. Hence,outliers (e.g., individual data points with large measurement errors as is often the case withfield measurements) can introduce a significant bias in the estimated model parameters.
Finsterle & Najita (1998) studied several other more robust estimators with different error
distributions that reduce the effect of outliers on the optimized parameters. For example, theysuggest using as an alternative the least absolute deviates or L1-estimator for a double
exponential distribution of errors:n(  ) = 

i= 1

1

qi*

i

- q ( )(3)ior the maximum-likelihood estimator for measurement errors following a Cauchy
distribution:

FG
H

n
1
(  ) =  log 1  wi [qi*2i=1

- q ( ) ] IJK
2

(4)

i

The vector of optimized parameters  typically contains parameters describing theunsaturated soil hydraulic properties (i.e., soil water retention and unsaturated hydraulicconductivity functions), and/or such solute transport parameters as the dispersivity and selected
adsorption or decay coefficients. The vector of observations typically contains measured
variables such as pressure heads, positions of the groundwater table, water contents,
concentrations, infiltration rates, and/or drainage fluxes.
The weighted least-squares estimator (1), as well as the presumably more robust estimators(3) and (4), can be easily modified to accommodate additional information, such as priorinformation about the optimized parameters, and/or known data points of the retention andhydraulic conductivity functions. For example, the objective function in the HYDRUS models
describing variably-saturated flow and transport (Šimůnek et al., 1998b, 1999b) is defined as thesum of three components:435 (  , q , p) =

mq

nq j

 v  w [ q ( x , t ) - q ( x , t , ) ] +
j

*
j

i, j

j=1

i=1

mp

np j

2

i

j

i

 v  w [ p ( ) - p ( ,  ) ] +
j

j=1

*
j

i, j

2

i

j

i

(5)

i=1

nb

 v [  -  ]
j

*
j

2

j

j=1

where the first term on the right-hand side represents deviations between the measured andcalculated space-time variables such as observed pressure heads, h, water contents, , and/or
concentrations, c, at different locations and/or time in the flow domain, or actual or cumulative
fluxes versus time across certain boundaries. In this term, mq is the number of different sets ofmeasurements, nqj is the number of measurements within a particular measurement set, qj*(x,ti)
represents specific measurements at time ti for the jth measurement set at location x(r,z), qj(x,ti,)represents the corresponding model predictions for the vector of optimized parameters , and vj
and wi,j are weights associated with a particular measurement set or point, respectively. The
second term on the right-hand side of (5) represents differences between independently measured
and predicted soil hydraulic properties (e.g., retention, (h), and/or hydraulic conductivity, K()or K(h), data), while the terms mp, npj, pj*(i), pj(i, ), v j and wij have similar meanings as forthe first term but now for the soil hydraulic properties. The last term of (5) represents a penalty
function for deviations between prior knowledge of the soil hydraulic parameters, j*, and their
final estimates, j, with nb being the number of parameters with prior knowledge and v jrepresenting pre-assigned weights. The best-fit parameters are obtained by minimizing thisobjective function using any minimization technique, such as the Marquard-Levenberg
method, the simplex method, or others (Šimůnek & Hopmans, 2000).

Water Flow
Water flow in variably-saturated flow models is typically described using a mixed
formulation of the Richards equation: A h  K izA  K (h) K ijt xi 
x j



(6)

where  is the water content (L3L-3), h is the pressure head (L), K is the hydraulic conductivity(LT-1), KA is the anisotropy tensor, t is time (T), and xi is the spatial coordinate (L). The soil
hydraulic properties can be described using the van Genuchten-Mualem model (vanGenuchten, 1980):

Se ( h ) = ( h ) - r
1=n m
 s - r(1+  h )

(7)

436

K (  ) = K s S el [ 1 - ( 1 - S e1 / m ) m ]2(8)

where Se is effective fluid saturation (-), Ks is the saturated hydraulic conductivity (LT-1), r ands denote the residual and saturated water contents (-), respectively; l is the pore-connectivity
parameter (-), and  (L-1), n (-), and m (= 1 - 1/n) (-) are empirical shape parameters. The above
hydraulic functions contain 6 unknown parameters: r, s, , n, l and Ks.
Field Applications

While inverse methods in subsurface hydrology were initially used almost exclusively forsaturated flow problems (e.g., see review of Yeh, 1986); they are now also increasingly applied
to the analysis of unsaturated zone experiments (e.g., Kool et al., 1987). Except for an early study
by Dane & Hruska (1983), parameter estimation techniques in soil hydrology were first usedwith laboratory experiments, such as for one-step or multi-step type outflow, upward infiltration,and evaporation experiments (see Hopmans et al., 2000, for references on these particularmethods) for which the initial and boundary conditions and homogeneity of the soil sample
could be readily controlled. The application of inverse methods to field experiments is more
complex due to soil heterogeneity, possible anisotropy, and uncertainty associated with the initialand boundary conditions. Examples are ponded infiltration (e.g., Russo et al., 1991), theinstantaneous profile method (Dane & Hruska, 1983; Kool et al., 1987; Romano, 1992), tension
disc infiltrometers (Šimůnek & van Genuchten, 1996, 1997; Šimůnek et al., 1998c, 1999c), cone
penetrometers (Gribb, 1996; Gribb et al., 1998; Kodešová et al., 1998, 1999; Šimůnek et al.,1999a), water flow to a ceramic cup tensiometer (Timlin & Pachepsky, 1998), and the extraction
method (Inoue et al., 1998). Contrary to laboratory methods, field experiments do not provide
exact control, or often do not even permit precise knowledge, of the total amount of water orsolute in the selected transport domain. Also, the use of the prior information is much more
complicated. For example, in the laboratory one can generally estimate with relatively great
precision the saturated water content from bulk density or porosity measurement. Suchinformation is much less useful in field studies where the soil seldom if ever reaches fullsaturation, even after development of positive pressure heads.
The Instantaneous Profile Method
The inverse parameter estimation problem was probably first applied to field data by Dane &Hruska (1983), who successfully optimized the parameters in van Genuchten's (1980) soil
hydraulic functions from transient drainage field data. A similar drainage problem was studied
by Kool et al. (1987) for a 6-m deep lysimeter, and later by Romano (1992) in attempts toquantify soil hydraulic spatial variability along a 50 m-long transect.
The Tension Disc Infiltrometer Method
Tension infiltrations are increasingly being used for evaluating saturated and unsaturated
hydraulic conductivities, and for quantifying the effects of macropores and preferential flow
paths on infiltration. A relatively standard way for estimating unsaturated hydraulicconductivities from tension infiltrometer data has been to invoke Wooding's (1968) analyticalsolution. This approach requires estimates of the steady-state infiltration rate for two different
supply pressure heads, and assumes applicability of an exponential function for K(h). Šimůnek &van Genuchten (1996) suggested the combined use of transient infiltration data obtained during a437

0.005-0.1

Hydraulic Conductivity [cm /s]

Cum ulative Infiltration [cm]

25
20
-1.0
15

-3.0
-6.0

10

-9.0
-11.5 cm

5

MeasuredFitted

0
0

1000

2000

3000 4000
Tim e [s]

5000

6000

Wooding's Analysis
0.003

0.002

0.001

0.000

7000

Figure 1. Measured and optimized cumulative infiltration curves for a tension disc
infiltrometer experiment carried out on a
sandy soil in the Sahel region.

Numerical Optimization

0.004

-1

-0.5

0

0.5
1
1.5
log(|h| [cm ])

2

2.5

3

Figure 2. Unsaturated hydraulic conductivities calculated using Wooding's analytical solution, and the complete function
obtained with numerical inversion.

single tension infiltration experiment, and tensiometer or TDR data measured in the soil below
the disc, to estimate the unknown soil hydraulic parameters via parameter estimation. They later
revised this method by using multiple tension infiltration experiments in combination withknowledge of the initial and final water contents (Šimůnek & van Genuchten, 1997). This
modification avoided the cumbersome use of tensiometers and TDRs. An evaluation of thenumerical stability and parameter uniqueness using numerically generated data withsuperimposed stochastic and deterministic errors showed that a combination of multiple
cumulative tension infiltration data, a measured final water content, and an initial condition
expressed in terms of the water content, provided the most promising parameter estimation
approach for practical applications (Šimůnek & van Genuchten, 1997).
Numerical inversion was later used to estimate the soil hydraulic properties of a two-layered
crusted soil system in the Sahel region of Africa (Šimůnek et al., 1998a). Here we will report
only results for the sandy subsoil obtained with a tension disc diameter of 25 cm and for supply
tensions of 11.5, 9, 6, 3, 1, and 0.1 cm. Figure 1 shows measured and optimized cumulative
infiltration curves. The small breaks in the cumulative infiltration curve were caused by brief
interruptions to resupply the infiltrometer with water, and to adjust the tension for a new timeinterval. Very close agreement between the measured and optimized cumulative infiltration
curves was obtained; the largest deviations were generally less than 60 ml, which was only about
0.5% of the total infiltration volume. Figure 2 compares the parameter estimation results against
results obtained with Wooding's analysis. Both methods gave almost identical unsaturated
hydraulic conductivities for pressure heads in the interval between -2 and -10.25 cm. However,
the hydraulic conductivity obtained with Wooding's analysis at the highest pressure headoverestimated the inverse results by a factor of two. Additional applications of inverse modelingto transient tension disc infiltration data are given by Šimůnek et al. (1998c, 1999c).
The Cone penetrometer method
While tension infiltrometer experiments provide relatively quick estimates of the hydraulicproperties, they can be used only at the soil surface. By comparison, a new cone penetrometer
method currently under development (Gribb, 1996; Gribb et al., 1998; Kodešová et al., 1998,
1999; Šimůnek et al., 1999a) can be used at depth. To obtain the hydraulic properties, a modified
cone penetrometer was instrumented with a porous filter close to the penetrometer tip and two438

tensiometer rings 5 and 9 cm above the filter. The device is pushed into a soil to the desired
depth, and a constant positive head is applied to the 5-cm filter. The volume of water imbibed
into the soil is monitored, as are tensiometer ring readings registering the advancement of thewetting front. After measured pressure heads in the soil stabilize, the water source is shut off,
and the redistribution of water in the soil profile is further monitored with tensiometers.
Cumulative infiltration and the measured pressure heads are used to define the objective function
to be minimized.
Gribb (1996) gave a detailed numerical analysis of this experiment, including a study of theidentifiability of the soil hydraulic parameters. The method was subsequently used to estimate
the hydraulic parameters of a sandy soil in a laboratory aquifer system measuring 5 x 5 x 3 m(Gribb et al., 1998; Kodešová et al., 1998) using the wetting part of the experiment. Šimůnek etal. (1999a) and Kodešová et al. (1999) later examined both the wetting and redistribution partsof cone permeameter experiments to find the wetting and drying branches of the soil
hydraulic properties. The study by Kodešová et al. (1999), carried out under field conditions,also evaluated the impact of having one or two steps in the applied pressure head on estimatesof the wetting soil hydraulic properties. The soil hydraulic properties estimated in this way
corresponded reasonably well with those obtained using more standard techniques.

Model Calibration
Variably-saturated water flow and solute transport models are increasingly being applied tonatural systems. Recent applications include studies of flow and transport in tile-drained
fields by de Vos (1997), de Vos et al. (1999a, b), and Mohanty et al. (1998a, b), among manyothers. The soil hydraulic and transport parameters in these studies were either obtained
independently of the model (Mohanty et al., 1998a, b) or were optimized by also considering
the independently measured prior information (de Vos, 1997; de Vos et al., 1999a, b). Herewe will briefly summarize the "split sampling" calibration of the HYDRUS-2D flow/transportmodel (Šimůnek et al., 1999b) to a tile-drained, layered silt loam soil (de Vos, 1997; de Vos etal., 1999a, b).
Case Study

The study was carried out at an experimental farm near Marknesse in the Noordoostpolder in
the Netherlands. The calcareous silt loam profile showed strong stratification (Figure 3) and was
drained using tile drains spaced 12 m apart and located at depths between 90 and 105 cm. The
soil profile was divided into three main horizons (topsoil between 0 and 25 cm, a 25-40 cm
intermediate layer, and a 40-200 cm subsoil layer), and also contained drain trench with slightly
different hydraulic properties. Within the intermediate and subsoil layers, distinctions were madebetween different sub-layers in terms of the saturated hydraulic conductivity. Tensiometers,
observation wells, piezometers, and a drainage water sampling device were installed in the field
to monitor the hydrological conditions, the drain discharge rate, and nitrate concentrations in thedrainage water.
Laboratory analyses included the hanging column method for measuring retention curves,the constant head method for the saturated hydraulic conductivities, and the crust method and the
hot-air method for unsaturated hydraulic conductivities. Measured soil-hydraulic parameters are
given by De Vos (1997) and de Vos et al. (1999b). Field water retention curves were obtainedby matching the field pressure head and water content measurements. Measured retention
curves and unsaturated hydraulic conductivities were used in the simulations to be discussed

439

Figure 3. Distribution of different soil layers in the flow domain, including the location of adrain trench above and below the drain.

Figure 4. Finite element grid in the flow domain representing half the drain spacing. The grid
system consists of 2563 triangles and 1352 nodes.
below. Since the saturated hydraulic conductivity varied by about two orders of magnitude
within any particular layer, values used for model simulations were obtained by using thecalibration procedure described below. A linear distribution of hydraulic conductivities was
assumed between saturation and the pressure head of h= -20 cm.
A flow domain 6 m wide, representing half the drain spacing, and 2 m deep was used in
the numerical simulations. The drain was located at a depth of 97.5 cm, and was described asa half circular hole of 5 cm diameter having a seepage face boundary condition. A stratified
triangular finite element grid was generated in accordance with the layered soil profile (Figure
4).
Drain discharge rate - groundwater level data, q(hGWL) (Figure 5), measured in the field,
were used to calibrate saturated hydraulic conductivities, Ks, of different layers. Values of Ks,
conditioned to remain within the range of measured values, were adjusted for different constant
upper boundary fluxes until the HYDRUS-2D simulation results corresponded with themeasured q(hGWL) data. To estimate Ks for the subsoil, the calibration procedure started with arelatively low uniform infiltration rate at the soil surface (0.5 mm d-1), which corresponded to a

440

Figure 5. Measured and simulated data of the relationship between drain discharge rate andgroundwater level for the leaching period from 17 December 1991 to 25 May 1992 (0<t<160d). The groundwater data represent the position midway between two drains, i.e., at 6 m fromthe drain. The simulated data were obtained for steady state flow conditions.
very low drain discharge rate. By increasing the infiltration rate, the phreatic surface rose tocause shallower soil layers to become saturated, and allowing their Ks values to be calibrated.
The calibration of the shallower layers depended on the results of the underlying layers. Thismethod of calibration assumes that the saturated zone below the phreatic surface dominateswater flow to the drain. We used a trial-and-error approach for calibration, since an automatic
calibration procedure was beyond the computational means of our PC's due to longcomputational times for each simulation. The model was further calibrated against a period of 40days starting December 17, 1991. Additional adjustments were made to the retention curve of the
topsoil on the basis of field-measured water retention curves (de Vos et al., 1999a).
The HYDRUS-2D model was subsequently used to simulate water flow and nitratetransport for a period of several years (up to April 94). Pressure heads and nitrate concentrations
measured on December 17, 1991 were used as initial condition. Figure 6 shows results of thesimulation for the first 160 days, including the first 40 days of the calibration period. Measured
and simulated groundwater levels at 6-m distance from the drain, drain discharge rates, andNO3 concentrations in the drainage water are presented. Notice that the modeling resultscorrespond equally well with the measured data during the calibration and validation parts ofthe experiment.

Conclusions
Recently developed variably-saturated water flow and solute transport models (e.g.,
HYDRUS-2D) can be extremely useful for analyzing a broad range of transient in-situ field flow
and transport experiments, because of their inverse capabilities and generality (in terms of defi-

441

Figure 6. Precipitation amounts (a), measured and simulated groundwater levels 6-m awayfrom the drain (b), drain discharge rates (c), and NO3 concentrations in the drainage water(d), for the leaching period from 17 December 1991 to 25 May 1992 (0<t<160 d).

442

ning the objective function, possible combinations of different initial and boundary conditions,and options for considering multi-layered systems). Such models represent very powerful tools
for analyzing complex environmental problems, and for developing alternative soil and water
management strategies, in spite of difficulties associated with their complete (mathematical)
verification and (experimental) validation. They provide a useful framework for interpreting
experimental results and for understanding qualitative and quantitative trends andrelationships present in observed field data. Transport models are excellent tools for
integrating our knowledge of the most relevant flow and transport processes and soil
properties, carrying out sensitivity analyses, and for designing, implementing and analyzinglaboratory and field experiments (Steefel & Van Cappellen, 1998; de Vos et al., 1999a;Šimůnek & Valocchi, 2000). This synthesis of modeling and experimentation should lead to amore coherent and rigorous understanding of the underlying transport processes.

References
Anderson, M. P. & W. W. Woessner, 1991. Applied Groundwater Modeling, Simulation of
Flow and Advective Transport. Academic Press, San Diego,CA, 381pp..
Anderson, M. P. & W. W. Woessner, 1992. The role of the postaudit in model validation.
Advances in Water Resources 15 [1], 167-173.
Bard, Y., 1974. Nonlinear Parameter Estimation. Academic Press, New York, N.Y., 341pp.
Dane J. H. & S. Hruska, 1983. In-situ determination of soil hydraulic properties during drainage.
Soil Sci. Soc. Am. J. 47 [3], 619-624.
de Marsily, G., P. Combes & P. Goblet, 1992. Comment on 'Ground-water models cannot bevalidated', by L. F. Konikow & J. D. Bredehoeft. Advances in Water Resources 15 [6],367-369.
de Vos, J. A., 1997. Water flow and nutrient transport in a layered silt loam soil. PhD thesis,
Wageningen Agricultural University, The Netherlands.
de Vos, J. A., D. Hesterberg & P. A. C. Raats, 1999a. Water flow and nitrate leaching in alayered silt loam soil. Soil Sci. Soc. Am. J., (submitted).
de Vos, J. A., J. Šimůnek & P. A. C. Raats, 1999b. Identification of the hydraulic properties of alayered silt loam. In: van Genuchten, M. Th., F. J. Leij, and L. Wu (eds.) Characterization andMeasurement of the Hydraulic Properties of Unsaturated Porous Media, University ofCalifornia, Riverside, CA, 783-798.
Finsterle, S. & J. Najita, 1998. Robust estimation of hydrologic model parameters. Water Resour.
Res. 34 [11] 2939-2947.
Gribb, M. M., 1996. Parameter estimation for determining hydraulic properties of a fine sandfrom transient flow measurements. Water Resour. Res. 32 [7], 1965-1974.
Gribb, M. M., J. Šimůnek & M. F. Leonard, 1998. Development of a cone penetrometer method
to determine soil hydraulic properties. ASCE J. of Geotech. and Geoenviron. Eng. 124 [9],820-829.
Hassanizadeh, S. M. & J. Carrera, (eds.) 1992. Validation of Geohydrological Models.
Special issues of Advances in Water Resources 15 [1,4].
Hopmans, J. W., J. Šimůnek & N. Romano, 2000. Inverse Modeling of Transient Water Flow.
In: J. H. Dane & G. C. Topp (eds.), Methods of Soil Analysis, Part 1, Physical Methods, Third
edition, SSSA, Madison, WI, in press.

443

Inoue, M., J. Šimůnek, J. W. Hopmans & V. Clausnitzer, 1998. In-situ estimation of soil
hydraulic functions using a multi-step soil water extraction technique. Water Resour. Rer. 34[5], 1035-1050.
International Atomic Energy Agency, 1982. Radioactive waste management glossary. IAEATECDOC-264, Vienna.
Jury, W. A. & H. Flhler, 1992. Transport of chemicals through soil: mechanisms, models,and field applications. Advances in Agronomy 47, 141-201.
Kodešová, R., M. M. Gribb & J. Šimůnek, 1998. Estimating soil hydraulic properties fromtransient cone permeameter data. Soil Science 163 [6], 436-453.
Kodešová, R., S. E. Ordway, M. M. Gribb & J. Šimůnek. 1999. Estimation of soil hydraulicproperties with the cone permeameter: Field studies. Soil Science, in press.
Kool, J. B., J. C. Parker & M. Th. van Genuchten, 1987. Parameter estimation for unsaturated
flow and transport models - A review. J. Hydrol. 91, 255-293.
Konikow, L. F. & J. D. Bredehoeft, 1992. Ground-water models cannot be validated.
Advances in Water Resources 15 [1], 75-83.
Mazoszewski, P. & A. Zuber, 1992. On the calibration and validation of mathematical models
for the interpretation of tracer experiments in the groundwater. Advances in Water
Resources 15 [1], 47-62.
Mohanty, B. P., R. S. Bowman, J. M. H. Hendrickx, J. Šimůnek & M. Th. van Genuchten,
1998a. Preferential transport of nitrate to a tile drain in an intermittent-flood-irrigated field:
Model development and experimental evaluation. Water Resour. Rer. 34 [5], 1061-1076.
Mohanty, B. P., T. H. Skaggs & M. Th. van Genuchten, 1998b. Impact of saturated hydraulicconductivity on the prediction of tile flow. Soil Sci. Soc. Am. J. 62, 1522-1529.
Narasimhan, T. N., 1987. Some thoughts on model verification. In: D. D. Evans & T. J.
Nicholson (eds.), Flow and Transport Through Unsaturated Fractured Rock, Geophysical
Monograph 42, AGU, Washington, DC.
OECD/Nuclear Energy Agency, 1990. The International Hydrocoin Project, Level 2: Model
Validation. OECD/Nuclear Energy Agency Publications No. 72957, Paris, France.
Oreskes, N., K. Shrader-Frechette & K. Belitz, 1994. Verification, validation, andconfirmation of numerical models in the earth sciences. Science 263, 641-646.
Romano, N., 1992. Use of an inverse method and geostatistics to estimate soil hydraulicconductivity for spatial variability analysis. Geoderma 60, 169-186.
Russo, D., E. Bresler, U. Shani & J. C. Parker, 1991. Analysis of infiltration events in relation todetermining soil hydraulic properties. Water Resour. Res. 27, 1361-1373.
Steefel, C. I. & P. Van Cappellen, 1998. Reactive transport modeling of natural systems. J. of
Hydrology 209, 1-7.
Šimůnek, J. & M. Th. van Genuchten, 1996. Parameter estimation of soil hydraulic properties
from the tension disc infiltrometer experiment by numerical inversion. Water Resour. Res. 32[9], 2683-2696.
Šimůnek, J. & M. Th. van Genuchten, 1997. Estimating unsaturated soil hydraulic properties
from multiple tension disc infiltrometer data. Soil Science 162 [6], 383-398.
Šimůnek, J., R. Angulo-Jaramillo, M. Schaap, J.-P. Vandervaere & M. Th. van Genuchten,
1998a. Using an inverse method to estimate the hydraulic properties of crusted soils fromtension disc infiltrometer data. Geoderma 86 [1-2], 61-81.
Šimůnek, J., M. Šejna & M. Th. van Genuchten, 1998b. The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 2.0, IGWMC-TPS-70, International Ground Water ModelingCenter, Colorado School of Mines, Golden, CO, 186pp.

444

Šimůnek, J., D. Wang, P. J. Shouse & M. Th. van Genuchten, 1998c. Analysis of a field tension
disc infiltrometer experiment by parameter estimation. Intern. Agrophysics 12, 167-180.
Šimůnek, J., R. Kodešová, M. M. Gribb & M. Th. van Genuchten, 1999a. Estimating hysteresis
in the soil water retention function from modified cone penetrometer test. Water Resour. Res.
35 [5], 1329-1345.
Šimůnek, J., M. Šejna & M. Th. van Genuchten, 1999b. The HYDRUS-2D software package forsimulating the two-dimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 2.0, IGWMC-TPS-56, International Ground Water ModelingCenter, Colorado School of Mines, Golden, CO, 251pp.
Šimůnek, J., O. Wendroth & M. Th. van Genuchten, 1999c. A parameter estimation analysis of
the laboratory tension disc experiment for determining soil hydraulic properties, Water
Resour. Res., in press.
Šimůnek, J. & J. W. Hopmans, 2000. Parameter optimization and nonlinear fitting, In: J. H.
Dane & G. C. Topp (eds.), Methods of Soil Analysis, Part 1, Physical Methods, Third edition,
SSSA, Madison, WI, in press.
Šimůnek, J. & A. J. Valocchi, 2000. Geochemical Transport, In: J. H. Dane & G. C. Topp (eds.),
Methods of Soil Analysis, Part 1, Physical Methods, Third edition, SSSA, Madison, WI, inpress.
SNPI (Swedish Nuclear Power Inspectorate), 1986. INTRACOIN Final Report, Levels 2 and3, Model Validation and Uncertainty Analysis. SKI 86:2.
SNPI (Swedish Nuclear Power Inspectorate), 1987. INTRAVAL Project Proposal. SKI 87:3.
Timlin, D. & Y. Pachepsky, 1998. Measurement of unsaturated soil hydraulic conductivities
using a ceramic cup tensiometer. Soil Sci. 163 [8], 625-635.
U.S. Department of Energy, 1986. Environmental assessment: Yucca Mountain site, Nevada
Research and Development Area, Nevada. DOE/RW-0073, Vol. 2, U.S. Department ofEnergy, Office of Civilian Radioactive Waste Mangement, Washington, DC.
van Genuchten, M. Th., 1980. A closed-form equation for predicting the hydraulic conductivityof unsaturated soils. Soil Sci. Soc. Am. J. 44, 892-898.
Wooding, R. A., 1968. Steady infiltration from large shallow circular pond. Water Resour. Res.
4, 1259-1273.
Yeh, W. W-G. 1986. Review of parameter identification procedures in groundwater hydrology:
The inverse problem. Water Resour. Res. 22 [2], 95-108.

445