Advances in Water Resources 24 (2001) 827±841
www.elsevier.com/locate/advwatres

Estimating unsaturated soil hydraulic parameters
using ant colony optimization
K.C. Abbaspour a,*, R. Schulin b, M.Th. van Genuchten caSwiss Federal Institute for Environmental Science and Technology (EAWAG), Ueberlandstrasse 133, P.O. Box 611, 8600 Dubendorf, Switzerlandb
Swiss Federal Institute of Technology, Institute of Terrestrial Ecology, ETH Zurich, SwitzerlandcGeorge E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 West Big Spring Road, Riverside, CA 92507, USA
Received 1 February 2001; received in revised form 25 March 2001; accepted 12 April 2001

Abstract
Although models are now routinely used for addressing environmental problems, both in research and management applications,the problem of obtaining the required parameters remains a major challenge. An attractive procedure for obtaining model
parameters in recent years has been through inverse modeling. This approach involves obtaining easily measurable variables (model
output), and using this information to estimate a set of unknown model parameters. Inverse procedures usually require optimizationof an objective function. In this study we emulate the behavior of a colony of ants to achieve this optimization. The method uses thefact that ants are capable of ®nding the shortest path from a food source to their nest by depositing a trail of pheromone during theirwalk. Results obtained with the ant colony parameter optimization method are very promising; in eight dierent applications wewere able to estimate the `true' parameters to within a few percent. One such study is reported in this paper plus an application toestimating hydraulic parameters in a lysimeter experiment. Despite the encouraging results obtained thus far, further improvements
could still be made in the parameterization of the ant colony optimization for application to estimation of unsaturated ¯ow andtransport parameters. Ó 2001 Elsevier Science Ltd. All rights reserved.
Keywords: Inverse modeling; Unsaturated soil; Hydraulic parameter; Ant colony optimization; Parameter estimation

1. Introduction
During the past several decades, considerable advances have been made in the mathematical descriptionof water ¯ow and solute transport processes in variably
saturated soils and groundwater systems. A large number of models of dierent degrees of complexity anddimensionality are now available for predicting subsurface ¯ow and transport [21,23,31]. Still, eective application of such models to practical ®eld problems suers
from the lack of knowledge of model parameters andtheir uncertainty. Model parameters are often dicult, if
not impossible, to measure directly because of time,money, instrumentation, scale, and conceptual constraints. To deal with the issue of model parameter
identi®cation, inverse modeling (IM) has provided an
attractive alternative method to direct measurement.

*
Corresponding author. Tel.: +41-01-823-5359; fax: +41-01-8235375.
E-mail address: abbaspour@eawag.ch (K.C. Abbaspour).

IM constitutes estimation of input parameters fromreadily measurable model output variables. Input
parameters are obtained by minimizing an objectivefunction describing the dierence between measured andsimulated data. Most optimization algorithms used forparameter estimation in variably saturated ¯ow andtransport studies [10,13,22,25,32] are gradient-type
methods that have the disadvantage of being very sensitive to the initial guesses of the parameters, while also
being prone to convergence to local minima [14,19]. In
this paper we propose a global search algorithm forestimating model parameters based on optimization by acolony of ants.
Social insects, such as ants, bees, termites and wasps,
often exhibit a collective problem-solving ability [5]. In
particular, ants are capable of ®nding the shortest path
between their nest and a food source [4]. This optimization capability is apparently achieved by secreting
recruitment pheromone on their path when the ants arein a transport mode [20].
The process by which ants are able to establish theshorter of two paths between two points A and B is0309-1708/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved.
PII: S 0 3 0 9 - 1 7 0 8 ( 0 1 ) 0 0 0 1 8 - 5

828

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

town only once by keeping track of towns visited in a``tabu'' list.
Let sij t be the intensity of the trail on edge i; j attime t. Each ant at time t chooses the next town, where it
will be at time t  1. Let us de®ne a ``run'' as Z moves
required by the Z ants in the interval t; t  1. For Y
towns it hence will take Y runs of the model to completean ``iteration'' during which each ant has completed atour. The trail intensity sij is updated according to theequation
(a)

(b)

(c)

Fig. 1. Example of real ants ®nding the shortest path: (a) initially antschoose paths ACB and ADB with equal probability; (b) path ACB,
being shorter, receives more pheromone with time; (c) more pheromone attracts more ants; hence, eventually all ants will follow pathACB.

sij t  1  qsji t  Dsij ;where q < 1 is a coecient such that (1 q) represents
the evaporation (loss) of trail between times t and t  1,
and
Dsij 

ZX
k1

illustrated in Fig. 1. Initially, ants have an equal chanceof choosing between paths ACB and ADB. However,after some time path ACB, by virtue of being shorter,
receives more pheromone from ants going in both directions (ACB, BCA). The shorter path with extra
pheromone tends to attract more ants, which in turn willlead to more pheromone accumulation. The collective
behavior that emerges is a form of autocatalytic (positive feedback) process, which will soon result in all ants
following the shorter path.
Inspired by this ant-based optimization principle,
researchers have proposed several algorithms for solution of classical combinational optimization problems
such as the traveling salesman problem (TSP) [15], thequadratic assignment problem [11], or the job-shop
scheduling problem [7], all of which were reasonably
successful. This paper presents an adaptation and application of the ant colony optimization (ACO) technique to the inverse estimation of unknown parametersin unsaturated ¯ow models.

2. The ant colonization search algorithm
Since, to the best of our knowledge, ACO has notbeen used before in the area of soil and hydrologicalsciences, we ®rst brie¯y review its benchmark application to the TSP [9]. In the TSP, the objective is to ®ndthe shortest path that connects (visits) each of Y towns
once. Let zi t; i  1; . . . ; Y  be the number of ants intown i at time t, and let Z be the total number of agents(ants). Each arti®cial ant in town i is an agent who
places a substance called ``trail'' (pheromone for realants) on a visited path, and then chooses to visit the next
town j with a probability that is a function of the town
distance, dij , and of the amount of trail present on theconnecting path. Also, each agent is trained to visit a1Dskij ;

2

where Dskij is the quantity per unit of length of trail
substance (pheromone in the case of real ants) placed onedge i; j by the kth ant between times t and t  1. The
quantity Dskij is given by
8Q< Lk if kth ant uses edge i; jDskij 3in its tour between t and t  1;:0 otherwise;
where Q is a constant related to the quantity of trail laidby ants, and Lk is the tour length of the kth ant.
To satisfy the constraint that an ant must visit all Y
towns, a data structure called the tabu list is associatedwith each ant; this tabu list saves the towns alreadyvisited up to time t and forbids an ant to visit the same
town again before Y runs (a complete tour) have been
completed. When a tour is completed, the tabu list is
used to compute the ants' current solution (i.e., thedistance of the path covered by an ant). The tabu list isthen emptied and the ant is free to again choose. We
de®ne tabuk as the dynamically growing vector that
contains the tabu list of the kth ant.
Let us de®ne the ``visibility'' gij as 1=dij (i.e., the reciprocal of the distance between two towns). Thisquantity does not change during a run of the ACO, as
opposed to the trail s that changes according to Eq. (1).
We can then de®ne the transition probability from town
i to town j for the kth ant as
8ef< P sij t gij 
if j 2 allowedk ;ksik te gik f
pij t 4
: k2allowedk
0
otherwise;
where allowedk  fY tabuk g and where e and f areparameters that control the relative importance of trail
versus visibility. The transition probability is a trade-o
between the visibility (which states that close towns will
be selected with a higher probability, thus implementing
a greedy heuristic strategy), and the trail intensity attime t (which states that paths with a lot of trac must

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

be highly desirable, thus implementing the autocatalytic
process). A greedy heuristic strategy is a strategy thathelps ®nd acceptable solutions in the early stages of thesearch process.
After Y runs all ants have completed a tour, and theirtabu lists will be full; at this point for each ant k thevalue of the tour length Lk is computed and the values
Dskij are updated according to Eq. (3). Also, the shortest
path found by the ants (i.e., mink Lk , k  1; . . . ; Z) is
saved and all tabu lists are emptied. This process is iterated until the tour counter reaches the maximum (userde®ned) number of iterations.
Dorigo et al. [8] compared ACO with several otherpopular algorithms for solving the Oliver30 problem (a30-town problem) and concluded that with few exceptions, ACO always found the best-known solution.

829

Fig. 2. Illustration of parameter discretization, values of objectivefunction g for dierent parameter values, gmin ; gcr , and parameter
updating based on the criteria g 6 gmin , and g 6 T .

3. Application of ACO to inverse modeling
The ®rst step in any IM study is the de®nition of an
objective function. Our objective in using IM is to ®nd aparameter set b that minimizes the general functiong b 

l
X
i1

wi

qi
ri X
X
j1

k1


wijk vijk

v0ijk

2

;

5

where g is the objective function, b is the vector of input
parameters, l is the number of observed variables, vijk is
the measured value of the ith variable at the jth locationand the kth time, v0ijk is a corresponding simulated value,
wi is the weight of the ith variables, wijk is the weightassociated with individual measurements, qi is the total
number of measurements over time and ri is the total
number of measurements over space for the ith variable.
The algorithm ACO±IM described below is an adaptation of a previously developed program, Sequential
Uncertainty FItting (SUFI) [3], to ACO. The main difference between the two programs is in the updating of
the prior parameter uncertainty as described in the following paragraphs.
The main concept behind the SUFI algorithm can besummarized in the following steps: (1) depict each unknown parameter by an interval based on available information, i.e., [100±800] for the imaginary parameter p
in Fig. 2; (2) discretize each interval into a number ofstrata, let the middle of each stratum represent thatstratum; (3) run the simulation model of choice for all(exhaustive sampling), or a randomly selected subset of
all the possible parameter combinations (an iteration);(4) score each stratum on the bases of the smallness ofthe value of the objective function in such a way that
small values of objective function receive higher scores;(5) eliminate from the ends of the parameter intervals
the strata with small or no scores, thereby decreasing theinitial range of parameter uncertainty, and ®nally; (6)repeat the process until a desired stopping rule isreached. In SUFI algorithm, each stratum, initially set
to a score of zero, receives a score of 1 if the value of theobjective function, g, for that stratum satis®es the condition g 6 gcr , where gcr is a user-de®ned value. In theillustrative example in Fig. 2, for the conditiong 6 gcr  gmin , where gmin is the minimum value of theobjective function obtained in an iteration, only one
stratum receives score, and hence the updated parameter(p10 ), based on this criterion, would have a very narrow
range [500±600]. This case will provide a fast but lessthan optimum solution as the global minimum (at about
400 in Fig. 2) will be missed. For the criteriong 6 gcr < T , where T has a value greater than gmin , theupdated parameter would look like the second line in
Fig. 2 (p20 ) with a range of [200±700], which is smaller
than the initial range of [100±800] but still much larger
than the g 6 gcr  gmin  criterion. For this case, there isa larger probability of ®nding the global minimum but alarger number of iterations must be performed.
For parameter updating, the user in SUFI program ispresented by score distributions for several gcr conditions, and is asked to update each parameter bystudying the scoring patterns of the dierent strata [3].
In the following paragraphs we present a more formal
description of the above steps, and point out that thescoring, and hence the updating of the parameters in
ACO±IM is determined by a function that tries toachieve a balance between dierent interests, mainly, adesire to reach a fast solution modi®ed by a desire to
avoid the local minima.
3.1. Formulating ACO±IM
The ACO±IM problem can be described as follows.
Let g b : D ! R, as given by (5), be a continuous and830

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

bounded function, and b be a b-dimensional state vector. Also, let R be real numbers, and D be the search
space of b. The objective is to ®nd the state vector b thatminimizes g b. Without loss of generality, we can takeD as the hyperparallelepipedD  fbi jbi 6 bi 6 bi ; i  1; 2; . . . ; bg;
where bi and b
i denote, respectively, the lower and theupper bounds of the parameter bi . The upper and thelower bounds are based on prior information using existing knowledge.
Since ACO has been developed for ordered problems,
its application to a continuous problem requires ®rst
discretization of the search space, D. This is done by
dividing the interval bi ; bi  of each parameter bi into a
number, say mi , of strata. If each stratum be represented
by the value at the middle of the stratum, then there will
be M  m1 m2 ; . . . ; mb permutations or possible pathways
through the space of input parameters. We refer to M
as the exhaustive sampling scheme (ESS). Let
bij i  1; . . . ; b and j  1; . . . ; mi  be the stratum j ofparameter bi in the search space of D.
The above discretization scheme is schematicallyshown in the right-hand side of Fig. 3 for four parameters, with each parameter interval bi ; bi  divided
into ®ve strata. There are no hard rules for choosing thestrati®cation of the parameters, except that initially theyshould uniformly cover the complete parameter space,
unless additional information exists that would favor
one region of parameter space above another. While
having more strata would speed up convergence of the
ACO problem, this advantage should be balanced
against the speed of the simulation program since thenumber of calls to a simulation program rapidly increases as the number of strata increases. To limit the
number of function calls, one could select only a random
subset of M, or N, which we refer to as a random
sampling scheme (RSS), and/or invoke parallel processing. The ACO±IM algorithm easily lends itself toparallel processing techniques. Based on our experienceFig. 3. Graphical representation of discretized parameter space on theright and the pathway structure on the left.

so far, N P 0:1 M generally leads to optimum solutions,while and N < 0:1 M often leads to less than optimum
solutions.
3.2. Assigning tasks to ant agents
To satisfy the constraint that each ant must follow
only one pathway, we generate a pathway-structure list
as shown in the left side of Fig. 3. Using this list, eachant will follow one of M, or N pathways, i.e., thepathway (1,1,1,2), or (2,4,2,4) as shown in Fig. 3.
Therefore, at the most we need a total of M ants, witheach ant acting like an agent having the following tasks:· selecting a pathway from the pathway list and remembering the values of the parameter strata along
its pathway,· at the end of the pathway, passing the parameter values to a simulation program,· getting simulated values back from the simulation
program,· calculating the value of the objective function for its
pathway, and· on the basis of the value of the objective function,
placing a certain amount of trail (pheromone in thecase of real ants) on each stratum visited along its
pathway.
Once all ants complete the above tour, one `ant-cycle'
or one `iteration' is completed. Eventually the trail will
be used to score a pathway. The low scoring pathways
will vanish (in the case of real ants the pheromone
evaporates on pathways with less trac) and newpathways evolve around high scoring paths.
3.3. Calculating the trail and scoring parameter strata
Let su I be the intensity of trail on pathwayu u  1; . . . ; N  and iteration I. For the example in Fig.
3, ant number 209 places the same amount of trail on
stratum 2 of the ®rst parameter, stratum 4 of the second
parameters, stratum 2 of the third parameter, and stratum 4 of the fourth parameter. In the TSP problem thetrail intensity was de®ned as a function of the length of apath. In ACO±IM we use value of the objective functionof a pathway to assign the trail. The trail intensity oneach path is calculated in such a way that the larger
values of the objective function (i.e., longer paths forreal ants) receive smaller trails. This is accomplishedaccording to the following equation
(
 

u gcr
exp 4:6 ggmin; gu 6 gcr ;gcrsu I 6
0;gu > gcr ;where gu is the value of the objective function for the uthpathway, gcr is a critical value (to be de®ned later) above
which su  0 (see Fig. 2), and gmin is the minimum of all
gu s u  1; . . . ; N  for an iteration. Note that in an

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

ACO±IM problem the trail intensity is no longer afunction of time, as it was in the TSP, and is de®ned per
iteration. The above equation states that pathways withvalues of the objective function larger than gcr receive no
trail (for real ants the trail on longer paths evaporates
quickly), while below gcr , the intensity of trail is exponentially larger for smaller values of the objectivefunction (for real ants a shorter path from nest to foodsource receives increasingly more pheromone). At
gmin ; su I is set equal to an arbitrary value of 100 (thisvalue would be equivalent to Q in the TSP, which isincidentally also equal to 100).
Note that any single stratum bij in D may be thecrossroad of many ant pathways. The trail share (calculated with Eq. (6)) of each stratum bij from eachpathway can be summed to yield /ij as
X/ij I 
su I;
7
u2crossing pathways

where crossing pathways are all the pathways that crossstratum bij . In addition, let rij (I) be the standard deviation of all values of the objective function involving
bij . The value of rij (I) depicts the sensitivity of the bij
stratum with respect to other strata in D. Fig. 4 illustrates a pattern of both /ij and rij distributions over the
strata of the hypothetical four-parameter example inFig. 3.
To assess the evolution of new pathways throughparameter spaces we must formulate a trade-o between
the following three considerations: (1) the autocatalytic
process that accounts for rapid discovery of good solutions. This objective is achieved through Eq. (6) where
pathways with smaller values of the objective function
receive exponentially higher trails, (2) the greedy heuristics strategy that helps one to ®nd acceptable solutions in the early stages of the search process. Thisobjective is achieved by selecting gcr to be as close as
possible to gmin (see Fig. 2), and (3) the distributed
computation strategy that avoids premature conver-

(a)

831

gence. This objective is achieved by giving more weight
to strata with higher sensitivities rij .
In Fig. 4(a), regions of larger trail intensity (shownwith larger bars) are generated by smaller values of theobjective function, and hence are more likely to containthe desired optimum solution. This suggests that we may
eliminate regions with zero or small amounts of trail (amove that may constitute a greedy behavior). On theother hand, strata with large sensitivities (Fig. 4(b)) may
not have been represented well by their middle points
and should not be eliminated.
In order to decide how to update the parameters, wede®ne the term `score', Sij , for each bij calculated by thefollowing expression similar to the transition probabilityof the TSP problem
AN/ij rij8
Sij  P PAN ;gcr Trijij /ij
where we let the value of A  1, and calculate N and T
by usingrij
N  cN ;lij
9rg
T  gmin  cT ;lg
where lij is the mean of trail of the bij strata, rij is thestandard deviation of trail on the bij strata, lg and rg
are, respectively, the mean and the standard deviation of
the objective function for an iteration, and cN and cT areconstants that were determined to be about 0.3 and 0.5,
respectively (see Appendix A for more details). The ratioof standard deviation to mean is known as the variability and our objective in scoring each stratum
through Eq. (8) is to assign a larger value to T if thevariability of the objective function is large (refer toFig. 2) and to assign also a larger value to N if thesensitivity rij  of the stratum bij is large, hence satisfying the criteria for distributed computation.
The parameters A, N, and T in (8) control the relative importance of the trail intensity (which implements
the autocatalytic process) versus sensitivity (which addresses the demand for distributed computation). Although the greedy nature of arti®cial ants dictates thatparameters be updated on the basis of the smallest valueof T (i.e., T  gmin ) alone, the addition of rij and /ij
in¯uences an ant's decision, the eect of which is similar
to the eect of a mutation operator in a genetic algorithm. [24] We refer henceforth to the above scoring
strategy as the ANT scoring strategy. A de®nition of theparameters involved in ACO±IM appears in Table 1.

(b)3.4. Updating of parameters and stopping rules
Fig. 4. Typical trail intensity distribution (a) and sensitivity distribution (b) for each parameter stratum. The larger bars indicate larger
trail values and higher sensitivities.

Once the scoring is completed, parameters are updated by eliminating the strata with small scores from

832

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

Table 1De®nition of the parameters involved in ACO±IM
ParameterDe®nition

A, N, T
cN ; cT
gmin
gcr
Sij
bij
lg ; rgsu/ij
rijParameters controlling the eect of /ij ; rij , and gcr , respectively, on Sij
Constants associated with parameters N and T
Minimum value of the objective function in an iteration
Critical value of objective function above which a pathway receives no trail
Score of stratum j of parameter iStratum j of parameter bi
Mean and standard deviation of objective function in an iteration
Amount of trail laid on pathway u by a crossing ant
Total amount of trail accumulated on stratum j of parameter iStandard deviation of objective function for stratum j of parameter iboth ends of the intervals. This results in parameterswith narrower ranges (uncertainties). Note that if thehigh scoring strata fall on either end of the parameterintervals, the parameters could be extended in that direction. This makes the initial guess of the parameterintervals not very crucial, and at the same time allows adeleted parameter stratum to re-enter the calculation.
The system is subsequently reinitialized with updated
parameters and the process repeated.
The iterations continue until a stopping rule is satis®ed. This may occur when (1) a desired value of theobjective function is reached, or (2) no more changes inthe value of the objective function are obtained in consecutive iterations, or (3) the number of precision digits
used to signify a parameter becomes more than the accuracy needed to represent that parameter value, or (4)the 95% con®dence interval of the variable distributions,
resulting from the propagation of uncertainty in theparameters, contains 95% of the measurements (plustheir errors).
The above distribution is obtained by Monte-Carlothrough the sampling of the ant paths at the end of eachrun by simply looking at the distribution of variablevalues corresponding to the dierent paths. The 95%
con®dence interval for a variable (95PCI) is obtained by
calculating the cumulative distribution and marking thevariable values at 0.025 and 0.975.
For environmental projects, where there are usually
measurements errors and the models often do not capture all of the relevant process and boundary conditions,
we propose using option 4 above as a more appropriate
stopping rule. In Fig. 5 we illustrate the concept of`conditioned parameter distribution'. The ®rst iterationof the ACO±IM is based on a prior estimate of theuncertainty domain of the model parameters andtherefore, because of the often large initial uncertainties,
the 95PCI of any model variable is large. This is shownby the space within the example curves `a' in Fig. 5. As
iterations proceed, the uncertainty domain of theparameters become smaller and the parameter distributions are more and more conditioned on the measurements of the variable(s) used in the objectivefunction. Thus, conditioned parameter distribution in

the context of this study refers to parameter distributions such that, when propagated stochastically througha simulation program, the 95PCI of the simulated
variable would contain more than 95% of the measured
points (i.e., curves `b' in Fig. 5). In ACO±IM, the iterations can be continued until the upper and the lower
limits of the 95PCI coincide with a single curve. This
curve is produced by a set of single-valued parameters
generally referred to as the `best-®t parameters' (curve`c' in Fig. 5). However, as illustrated by curve `c', ®tted
parameters produce simulations that often miss most of
the measured points. In our opinion, ®tted parameters
are inadequate for analysis of environmental problems if
used without the uncertainty associated with them. In
gradient type optimization programs such as Levenberg±Marquardt (LM) [17], the covariance of theparameter matrix and hence the 95% con®dence interval
associated with each parameter are based on linear regression analysis. The problem, therefore, with the calculated parameter uncertainties are that they hold only

Fig. 5. Simulations based on the `conditional parameter distribution'
concept vs. the `best-®t parameters'. Best-®t parameters produce simulations that usually miss most of the measurements, that is, curve d,whereas conditioned parameter distributions produce simulations that
respect most of the measurements, that is, curves c.

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

approximately for the non-linear analysis [13]. Instead,to account for parameter uncertainties, we suggest obtaining conditioned parameter distributions, where mostof the data points are respected within the 95PCI.
Furthermore, the problem of algorithm convergence,as the result of correlation among parameters is perhaps less severe in obtaining conditioned parameter
distributions than ®tted parameters. This is because theformer considers ranges rather than point values.
However, as with any other algorithm, in case ofstrong parameter correlation in ACO±IM it would also
be necessary to constrain one of the dependent
parameters.

4. Testing ACO±IM4.1. Test case 1: application of ACO±IM to an unsaturated ¯ow problem
To test the ACO algorithm, we used Example 2 fromthe HYDRUS-1D and HYDRUS-2D Windows-based
software packages [26,27]. We chose this example,®rstly, because it considers a simple and yet non-trivial
one-dimensional water ¯ow in a two-layered ®eld pro®leof the Hupselse Beek watershed in the Netherlands, andsecondly, because we can assume that we know the `true'
parameters and test how closely ACO±IM can estimatethem.
The two-layered system being simulated is schematically described in Fig. 6. The equation being solved is a
one-dimensional isothermal Darcian ¯ow in a variably
saturated rigid porous medium, expressed by the following form of the Richards equation:


oh
o
oh

K K
S;
10
ot ox
ox
where h is the volumetric water content cm3 cm 3 , h is
the pressure head (cm), S is the sink term d 1 , K is theFig. 6. Flow system and the ®nite element mesh for the ®rst example.

833

hydraulic conductivity cm d 1 , z is a vertical coordinate (cm), and t is time (d). Atmospheric data and observed groundwater levels provided the required
boundary conditions for the numeric model. The soil
pro®le consisted of two layers: a 40-cm thick A-horizon,
and a B/C-horizon that extended to a depth of 300 cm(Fig. 4). The depth of the root zone was 30 cm. The soil
surface boundary conditions involved actual precipitation and potential transpiration rates for a grass cover.
The surface ¯uxes were incorporated by using average
daily rates distributed uniformly over each day. The
bottom boundary condition consisted of a prescribeddrainage ¯ux-groundwater level relationship given by[27]
q h 

a1 exp a2 j h

G j;

11

where a1 and a2 are empirical parameters, and G is thereference position of the groundwater level.
The initial moisture pro®le was taken to be in equilibrium with the initial groundwater level at 55 cm below
the soil surface. Soil hydraulic characteristics were described according to the van Genuchten model [29] forthe water retention curve,
Se h h h hr1m
hs hr
1  jahjn 

h<0

12

and
h h  hs

hP0

13

and the van Genuchten±Mualem model [18] for thehydraulic conductivity function(m  2
Ks Se0:5 1
1 Se1=m; h < 0;K h 
14
h P 0;
Ks ;where a and n are shape parameters, Se is eective water
saturation, m  1 1=n; hr and hs are the residual andsaturated water contents, respectively; and Ks is thesaturated hydraulic conductivity.
The water retention and ¯ow curves for both soils areshown in Fig. 7. For the inverse problem we assumed
that the following six parameters were unknown:
a1 ; n1 ; Ks1 a2 ; n2 ; and Ks2 . The information usually available for inverse analyses are some variables (e.g., pressure head, water content, concentration, or ¯ux from thebottom of a soil column or a lysimeter) measured atcertain depths and at certain times. In this test case wearbitrary assumed to have 20 pressure head datameasured on days 120, 160, 200, 240, 270, and at nodes
11, 25, 40, and 50; and 12 water content data measured
on days 130, 180, and 220, and at nodes 11, 21, 51, 63.
These `measurements' were taken from the output ®lesof Example 2 of HYDRUS_2D. Then, the objectivefunction was written as834

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

Fig. 7. Unsaturated hydraulic properties of the ®rst and second layers of the ®rst example.

g a1 ; n1 ; Ks1 ; a2 ; n2 ; Ks2  6:72  10 6

4 X
5 
X
hjk
j1

 40:48

k1

4 X
3 
X

hjk

j1

k1

h0jk

h0jk

2

2

;

15

where h is the measured pressure head, h0 is the simulated pressure head, h is the measured water content, h0
is the simulated water content, and the weights were
calculated according to the relationship 1=kr2  for eachvariable with r being the standard deviation of a variable and k being the number of samples [6].
Initially, we assumed a large range of uncertainties
for the unknown parameters, as listed in Table 2. We
estimated the unknown parameters using various strati®cation strategies and various number of searching ant
agents. All tests with the number of ant agents exceeding
10% of the exhaustive set M led to basically the same
results (2% dierences in the estimated parameters).
In the run reported here, initially we divided eachparameter into 5 strata, which produced a total ofM  15 625 exhaustive runs. But we executed only 20%of the exhaustive set and used N  3125 ant agents todo the search. We stopped after the sixth iteration withthe results summarized in Table 2. The entire run took
about two hours on a Pentium 600 MHz computer. The
strati®cation and the scoring for the greedy and the
ANT strategy are shown in Table 3 for the ®rst iteration. Note once again that the dierence between
greedy and ANT strategy is that for greedy strategy weonly use gmin as the only criteria for scoring the parameter strata. Hence, only one stratum per parameterreceives score, while for the ANT strategy more strata
combinations receive scores. The consequence of this isthat the updated parameters based on the greedy
strategy are more compact than a corresponding ANT
strategy. For convenience, in ANT strategy, we divided
the scores of all strata for each parameter by thesmallest score; hence, there is always one stratum with ascore of 1.
After the ®rst iteration, we followed separately thegreedy path and the ANT path for parameter updating.
Table 2 shows the results of both paths. For the greedy
path the third and the fourth iterations produced thesame value of the objective function, hence we considered the third iteration as the ®nal one. While for the
ANT strategy we stopped after the sixth iteration whenthe measured and expected simulated pressure head andTable 2
Initial parameter uncertainty range and ®nal estimates obtained using the ANT and the greedy strategies
Unknown parameter(true values)

Initial uncertainty
range

Final parameter and uncertainty
range based on ANT strategy

Final parameter and uncertainty
range based on the greedy strategy

a1 (0.0174)[0.0±0.05]

n1 (1.3757)[1.0±2.0]

Ks1 (29.75)[0.0±200.0]

a2 (0.0139)[0.01±0.03]

n2 (1.6024)[1.5±3.0]

Ks2 (45.34)[0.0±200.0]

Objective function

0.53

0.0173[0.0172±0.0175]
1.38[1.371±1.382]
28.85[28.74±29.82]
0.01386[0.01381±0.01393]
1.604[1.5983±1.6044]
44.6[44.43±45.76]
0.000019

0.02425[0.022±0.028]
1.29[1.28±1.36]
96.0[80.0±112.0]
0.01435[0.0112±0.0148]
1.65[1.56±1.8]32.8[16.0±40.0]
0.029

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841
Table 3
Summary of parameter strati®cation, strata scores, and updated
parameters based on greedy and ANI strategies for iteration 1 ofexample 1
Iteration 1
Parameter strati®cation
0.0
a1
n1
1.0
Ks1
0.0
0.01
a2
n2
1.5
Ks2
0.0
Greedy strategy:
a1
n1
Ks1
a2
n2
Ks2

0
0
0
1
1
1

0
1
0
0
0
0

0.01
1.2
40.0
0.016
1.8
40.0

0.02
1.4
80.0
0.022
2.1
80.0

Score
1 0 0
0 0 0
Updated parameters
1 0 0 !
0 0 0
0 0 0
0 0 0

0.03
1.6
120.0
0.028
2.4
120.0

0.04
1.8
160.0
0.034
2.7
160.0

0.05
2.0
200.0
0.04
3.0
200.0

a1  0:02±0:03
n1  1:2±1:4
Ks1  80±120
a1  0:01±0:016
n2  1:5±1:8
Ks2  0±40

ANT strategy:
a1
n1
Ks1
a2
n2
Ks2

0
0
1
8
4
1

3
7
1
5
2
1

Score
2 1 0
4 1 0
Updated parameters
1 1 1 !
3 2 1
2 1 1
1 1 1

a1  0:01±0:04
n1  1:2±1:8
Ks1  0±200
a1  0:01±0:034
n2  1:5±2:4
Ks2  0±200

water content curves were practically indistinguishable
as shown in Fig. 8. The sequential nature of the algorithm allows us to continue iterations until arriving at
the exact solution. This would not be possible with thegreedy strategy.
It appears that the ANT scoring strategy, while using
the greedy force through T to obtain a solution at theearly stages of the search, receives proper signals fromthe autocatalytic parameter A and the distributed computation parameter N to jump over local minima andarrive at an optimum solution.
Fig. 9 illustrates the 95PCI for the pressure head andwater content based on the values of the third, ®fth, and835

sixth iterations at a depth of 10 cm below the soil surface
as a function of time. Also shown in Fig. 9 are themeasured values and the expected values of the simulated variables for the sixth iteration. As iterations
continue, the uncertainty ranges become smaller, leading to narrower 95PCI. Note that for this ``well-de®ned''
problem (i.e., the measurements were model outputs
without any measurement errors, the boundary conditions as well as the processes were exactly de®ned) the
expected simulation (i.e., at the 50% con®dence interval)
along with the 95PCI for the sixth iteration all coincidewith the measured data. This normally would not happen in a real-world application.
4.2. Test case 2: comparison of ACO±IM with LM
methodIn this section we describe a comparison of the resultsof ACO±IM with the LM [17] method implemented in
HYDRUS-1D [26]. The example used is based on an
experiment carried out in a large lysimeter (3 m2 by 1.5m) planted with forest vegetation designed to study theeect of elevated CO2 and nitrogen on forest ecosystems[28]. We chose this speci®c example because the data has
all the noise that one would expect in a ®eld experiment
but with controlled boundary conditions. The important
processes included ¯ow, transport of bromide, and root
water and solute uptake. Another reason for using thisexample was to demonstrate that IM could be morethan parameter ®tting, and has the potential for being
used for system analyses.
The bromide CaBr2  tracer experiment began inJune 1996, when 30 g of bromide was diluted in 10 l ofwater, and distributed homogeneously below the canopy. The drained water was collected for a period ofover one year after the tracer application and the bromide concentration was measured on a regular basis.
Pressure head was measured with vertically installed
tensiometers at 25, 50, and 75 cm depth. Three replicate
instruments were installed at each depth and theirFig. 8. Measured (symbols) and simulated pressure head and water content obtained by the greedy and ANT scoring strategies.

836

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

Fig. 9. Graphs of the 95% con®dence intervals for pressure head and water content for the third, ®fth, and sixth iteration showing the progression toa solution.

averages used to represent the pressure head at that
depth. Water content was measured with time-domainre¯ectometry (TDR) probes of 25 cm length, verticallyinstalled with centers at depths 25, 50, and 75 cm. Again,
all installations were in triplicates, and the averages wereused to represent the water content at a depth. The
lysimeter was regularly irrigated and actual evapotranspiration calculated using the following water balance
equationETa  I

D

Dhz ;

16

where ETa is actual evapotranspiration, I is irrigation, Dis drainage, and Dhz is the dierence in water storage, allexpressed as amounts cumulated for the period of the
ETa measurement.
The governing ¯ow equation and the soil hydraulicproperties were as described for the previous test case.
Transport was described by a convection±dispersion
expression for a miscible, non-sorbing, and conservative
component as given by


ohc
o
oc
oqc

hD
SCs ;
17
ot
oz
oz
oz
where c is the solute concentration in solution g cm 3 ,
D is the dispersion coecient cm2 d 1 , q is Darcian¯ux cm d 1 , S is the root water uptake, and Cs is theconcentration of soil solution taken up by plant roots(g cm 3 ).
The ¯ow boundary conditions were expressed as:


oh
K
18
1  qT t; z  0;
oz


oh
19
K
1  qB t; z  L;oz
where qT is the net in®ltration cm d 1  at the top, andqB is the measured ¯ux cm d 1  at the bottom boundary. The initial and boundary conditions for the concentration were expressed as:c z; t  0; t  0;ochD  qc  qc0 t;ozoc 0; z  L;oz20z  0;

21
22

where c0 is the concentration of the in®ltrating water,and q0 is the Darcian ¯ux density evaluated at the soil
surface.
The objective function was formulated as
g a1 ; n1 ; Ks1 ; a2 ; n2 ; Ks2  1:13  10 6

3 X
41 
X
j1

 4:16

3 X
41 
X
j1

k1

h0jk

hjk

k1

hjk

h0jk

2

;

2

23

where the three space measurements were furnished withtensiometer and TDR data at 25, 50 and 75 cm depths.
Since the TDR data represented an average over 25 cm
vertical depth, the model outputs were also averaged
over 25 cm. The measurement times for both variables
were on days (13, 14, 16, 17, 20, 24, 26, 33, 40, 48, 54, 58,
65, 73, 79, 90, 97, 104, 111, 119, 126, 133, 139, 155, 162,
169, 174, 279, 284, 291, 300, 306, 312, 319, 327, 333, 340,
347, 354, 361, 368) where day 13 corresponded to June5, 1996.
Although the lysimeter was originally uniformly
packed, as the result of settling and root distribution amore compacted bottom layer was formed at about 40
cm depth. Therefore, we considered two layers and estimated six parameters as in the test example 1, mainly,
a1 ; n1 ; Ks1 ; a2 ; n2 ; and Ks2 . The strati®cation strategy fora; n; and Ks were 4, 4, and 10, respectively, and we
simulated only 10% of the exhaustive set and used onlyN  2560 ant agents.
Table 4 summarizes the initial and ®nal parameters
obtained by ACO±IM and the ®nal results obtained byK.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

837

Table 4
Initial parameter uncertainty range and ®nal estimates obtained using Levenberg±Marquardt and ACO±IM algorithms
Unknown parameter

Initial uncertainty distribution
used in ACO±IM

Expected parameter values and the®nal uncertainty distribution

Fitted parameter and the 95% CI
determined by Levenberg±Marquardt

a1[0.0±1.0]n1[0.1±3.0]

Ks1[300±100 000]a2[0.0±1.0]n2[1.1±3.0]

Ks2[100±100 000]

Objective function

0.53

0.937
[0.5±1.0]
1.288[1.1±1.6]
39 110[300±60 000]
0.4062[0.25±0.5]1.288[1.1±1.6]360[100±7000]
0.65

0.867
[0.478±1.256]
1.282[1.242±1.321]
31 627[1804±61 450]
0.124
[0.077±0.172]
1.518[1.389±1.648]
159[13±305]
0.63

LM algorithm along with the associated uncertainties.
Fig. 10 shows the comparison of LM and ACO±IM in
simulating pressure head and water content at the depthof 25 cm. For the other depths Figures are not shown
but the comparisons were similar with both models
having practically indistinguishable curves. Only slightdeviations of water content at 50 and 75 cm depths resulted in ACO±IM having a slightly larger objectivefunction. A review of the parameters in Table 4 reveals
that the identical simulation results in Fig. 10 were obtained with widely dierent parameters indicating thatthe set of hydraulic parameters considered was not very
important in this particular study. This is also evident
from the large uncertainties both algorithms determined
for the parameters.
We began the ACO±IM simulation with very large
initial parameter uncertainties so as to capture all themeasured data within the 95% con®dence band. But thiswas never achieved for the entire pro®le regardless of thedegree of parameter uncertainties. Figs. 11 and 12 show
the 95PCI for pressure head and water content for thesecond iteration at the depths of 25 and 75 cm, respectively. Although most of the data at the 25-cm depth

was captured within the 95PCI, but for the 50-cm depth(Figure not shown) and the 75-cm depth a large numberof pressure head data was missed. For this example weconsider further conditioning (iteration) of parameterson the measured data not to be justi®ed, and the parameters reported in Table 4 with the large range of
uncertainties associated with them are the best results
one could obtained under the given conditions.
To test the hydraulic parameters obtained in Table 4,
we simulated the breakthrough of bromide. Keeping thehydraulic parameters constant, we ®tted two transportparameters, dispersivity and maximum plant uptake.
Again, these parameters were not very sensitive, andirrespective of the initial uncertainty distribution, the95PCI could not contain the measured bromide concentration in leachate water. The measured values alongwith the simulations with LM and ACO±IM are shownin Fig. 13. It is evident that the measured Br breakthrough curve exhibits a more complex behavior thanwhat was simulated. One explanation for the lack ofsigni®cance in the parameters is that this particular
problem, contrary to the previous example, is processrather than parameter-controlled. In other words,

Fig. 10. Comparison of the ®nal simulation results obtained by LM and ACO±IM algorithms for the second example.

838

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

Fig. 11. The 95% con®dence interval for the pressure head and water content based on the second iteration at the 25-cm depth.

Fig. 12. The 95% con®dence interval for the pressure head and water content based on the second iteration at the 75-cm depth.

Fig. 13. Measured and simulated Br breakthrough curves.

variables water content, pressure head, and concentration are also partly governed by processes not considered by the simulation model. This revelation by inverse
analysis is quite useful and indicates that IM can be used
for process analysis as well as parameter estimation as
also shown by the authors elsewhere [1]. In this example
it would be meaningless to condition any further theparameters on the data and expect the parameters toaccount for responses that were caused by processes
other than those included in the simulation model.
The peak dominated Br breakthrough curve is indicative of preferential ¯ow in the soil system. This
could also be inferred from the very large estimates ofthe saturated hydraulic conductivity produced by bothinverse algorithms for the ®rst layer. If there are processes in the system not accounted for by the simulation
models, inverse models try to compensate for missing
processes by unrealistically over- or underestimatingsome of the parameters. This is a pitfall of IM that needs
more attention. Furthermore, the progressively
worsening of simulations as a function of increasing
depth could be a sign of inadequate modeling of root
distribution and, hence, water and passive solute uptake
by plants. Water content simulations are relatively better than pressure heads for all depths (Figs. 11 and 12)
because they are averages over 25 cm depths.
Finally, some of the advantages and disadvantages ofthe dierent IM algorithms can be summarized as below: (1) In ACO±IM degree of parameter conditioning
may be directly controlled through iterations. The moreiterations we perform the more conditioned theK.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

parameters being sought will become on the set ofmeasured data. In general, we do not recommend obtaining highly conditioned parameters, as they will
perform poorly when simulating other variables notincluded in the objective function [2]. In LM algorithm,
degree of conditioning can only be controlled indirectly
by restarting the program with several initial guesses,
and in general, a local minimum nearest to the initial
estimates is produced. (2) In LM algorithms parameter
estimates are given as single-valued parameters plus asummary statistics that expresses the parameter uncertainty usually in terms of the 95% con®dence interval.
The parameter uncertainties need yet to be propagatedbefore we can assess the degree of conditionality, while
in ACO±IM the uncertainties, which play a central role,
are already propagated. After each iteration, the 95PCI
is reviewed as a guide for further iteration. (3) The LM
method requires many trials with initial guesses of theunknown parameters and at the end one is never surewhether to perform another iteration of not. While in
ACO±IM this uncertainty does not exist. (4) In ACO±
IM, the parameters A, N, and T need yet a more robust
estimate, the process is more interactive than LM (which
may be considered a disadvantage by some, although wefavor this interaction), and subjective interference withthe parameter updating may be required, which calls forsome understanding of the problem and depends, morethan the LM method, on the experience of the user. The
subjective interference may be in the form of setting
upper and lower limits for the parameters, although thisoption also exists in the inverse option of the HYDRUS-1D program. Finally, in the course of the above
comparison, we realized that the two algorithms could
actually compliment each other in a way that an initial
run with ACO±IM could identify the general location of
the initial parameter estimates, which can then be fed to
LM to obtain a ®nal solution.
5. Conclusions
We introduced a new algorithm for inverse parameter
estimation using an optimization procedure inspired byant colony foraging behavior. The results presented inthis paper are very encouraging for ACO±IM method.
In the ®rst test, using a well-de®ned problem, the exact
parameters were identi®ed. While in the second test, a
comparison of the ACO±IM and the popular LM
algorithm produced practically identical simulation
results.
Based on the elementary behavior of an individual ant,
the greedy force combined with the autocatalytic process
would lead to sub-optimal solutions. However, when an
entire colony operates, complicated dynamics are exhibited through distributed computations where greed gives
the right signals to the autocatalytic process to converge839

to an optimum solution. Introduction of the parameters
A, N, and T achieved a balance between dierent dynamics that enables ants to jump over local minima.
Despite the encouraging results obtained in this ®rst
study, application of the ant optimization technique toenvironmental parameter estimation problems is still at
its infancy. The ACO±IM system could bene®t from
improved computational eciency, as well as a morerobust ANT scoring strategy. This is due to the fact thatthe values of N and T through cN and cT coecients
seemed to vary slightly for dierent cases studied. Future application of the ACO may need to address theseissues in more detail.
Appendix A
Parameters A, N, and T in Eq. (8) control the relative
importance of the trail, sensitivity, and critical value ofthe objective function. For convenience and in order tohave less parameters to determine, A was set equal to 1,
N was expressed in terms of the variability of the trail
received by each stratum, and T was expressed in termsof the variability of the objective function as follows:rij
N  cN ;lij
T  gmin  cT

rg
;
lg

where lij is the mean of trails placed on the bij strata, rij
is the standard deviation of trail on the bij strata, lg andrg are, respectively, the mean and the standard deviationof the objective function for an iteration, and cN and cT
are constants. Note that gmin is the smallest value that T
can obtain. For large variability in the objective function, the value of T should be larger in order to have abetter chance of avoiding local minima (refer to Fig. 2).
To determine the values of cN and cT we used threedierent models and tested three dierent problems. Fulldetails of the procedure can be found in Lerch [16]. The
three models used were HYDRUS_2D [27], HYDRUS1D [26], and MACRO [12]. The ®rst test was similar to
the example 1 in this manuscript with the objectivefunction written in terms of the water content andpressure head. In this test we assumed the parameterssaturated hydraulic conductivity Ks , and van Genuchten shape factors a, and n of the ®rst layer were
unknown. For the second test we used the program
HYDRUS-1D and simulated the breakthrough of aboron H3 BO4  pulse through Glendale clay loam [30].
In this example we treated the parameters saturated
hydraulic conductivity Ks , adsorption coecient ks ,and fraction of exchange sites in equilibrium withsolution concentration ( f ) as unknown and minimized
an objective function based on the breakthrough840

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841

Fig. 14. Normalized dierence between actual and modelled variables as a function of cT and cN . Typical run showing cT  0:5 and cN  0:3 to be
the best values.

concentration of boron. The third example, simulated
by MACRO [12], dealt with the break through of bromide from an unsaturated soil column. In this example
we estimated boundary hydraulic conductivity (KSM),
pore size distribution index (ZLAMB), and tortuosity
factor (ZN) and minimized an objective function basedon the breakthrough concentration of bromide.
In each example, we ®rst ran the programs with a setof known parameters, noted the desired simulated
variable as `observation', assumed an uncertain range
for the parameters we wished to treat as `unknown', anddetermined the sensitivity of cT and cN in a classical way
by keeping one of them constant and varying the other.
After some preliminary runs, we tested values:cT 2 f0:4; 0:45; 0:5; 0:55; 0:6; 0:65; 0:7g, and cN 2
f0:1; 0:2; 0:3; 0:4; 0:5; 0:7g. For each case we determined the `unknown' parameters by running ACO±IM
algorithm. Our objective was to obtain a set of cT and cN
parameters that would in each case lead to parameters
closest to the real parameters.
Fig. 14 illustrates an example of the results obtained
for the ®rst test. This Figure shows the normalized differences between the estimated and actual parameters by
ACO±IM versus dierent values of cT and cN . With
small variations, in each case the best values of cT and cN
were centered around 0.5 and 0.3, respectively. The
above values were used for examples one and two in the
manuscript. However, we feel because of the small deviations in cT and cN in the three examples considered,more work may be needed to identify a better formulation for A, N, and T and a more robust estimation ofcT and cN parameters.
References[1] Abbaspour KC, Kasteel R, Schulin R. Inverse parameter estimation in a layered unsaturated ®eld soil. Soil Sci 2000;165(2):109±5123.
[2] Abbaspour KC, Sonnleitner MA, Schulin R. Uncertainty in estimation of soil hydraulic parameters by inverse modeling: example
lysimeter experiments. Soil Sci Soc Am J 1999;63(3):501±9.

[3] Abbaspour KC, van Genuchten MTh, Schulin R, Schlappi E. A
sequential uncertainty domain inverse procedure for estimating
subsurface ¯ow and transport parameters. Water Resour Res
1997;33(8):1879±92.
[4] Beckers R, Deneubourg JL, Gross S. Trails and U-turns in theselection of the shortest path by the ant Lasius niger. J Theor Biol
1992;159:397±415.
[5] Bonabeau E, Theraulaz G, Deneubourg JL, Aron S, Camazine S.
Self-organization in social insects. Trends Ecol Evol 1997;12:
188±93.
[6] Clausnitzer V, Hopmans JW. Non-linear parameter estimation:
LM_OPT. General-purpose optimisation code based on Levenberg±Marquardt algorithm. Land Air and Water Resources Paper
No. 100032, University of California Davis, CA, USA, 1995.
[7] Colorni A, Doringo M, Maniezzo V, Trubian M. Ant system forjob-shop scheduling. Belg J Oper Res, Stat Comput Sci
1993;34:39±53.
[8] Dorigo M, Maniezzo V, Colorni A. The ant system: optimization
by a colony of cooperating agents. IEEE Trans Syst Man Cybern± Part B 1996;26(1):29±41.
[9] Dorigo M. Optimization learning and natural algorithms. PhD
Thesis, Dip. Elettronica e Informazione, Politecnico di Milano,
Italy, 1992.
[10] Eching SO, Hopmans JW. Optimization of hydraulic functions
from transient out¯ow and soil water pressure head data. Soil Sci
Soc Am J 1993;57:1167±75.
[11] Gambardella LM, Doringo M. Ant-Q: A reinforcement learning
approach to the traveling salesman problem. In: Prieditis A,
Russell S, editors. Proceedings of ML-95, 12th International
Conference on Machine Learning, 1995; Morgan Kaufmann,
Sanfransisco, CA. p. 252±60.
[12] Jarvis NJ. The MACRO model: Technical description and sample
simulations. Reports and Dissertation 19, Department of Soil
Science, SLU, Uppsala, Sweden, 1994.
[13] Kool JB, Parker JC. Analysis of the inverse problem for transient
unsaturated ¯ow. Water Resour Res 1988;24(6):817±30.
[14] Kool JB, Parker JC, van Genuchten MTh. Parameter estimation
for unsaturated ¯ow and transport models ± a review. J Hydrol
1987;91:255±93.
[15] Lawler EL, Lenstra JK, Rinnooy-Kan AHG, Shmoys DB,
editors. The traveling salesman problem. New York: Wiley; 1985.
[16] Lerch M. Determination and sensitivity analyses of ant colony
optimization ± inverse modeling parameters. Diploma Thesis,
Swiss Federal Institute of Technology, Zurich, 2001.
[17] Marquardt DW. An algorithm for least-squares estimation of
non-linear parameters. SIAM J Appl Math 1963;11:431±41.
[18] Mualem Y. A new model for predicting the hydraulic conductivityof unsaturated porous media. Water Resour Res 1976;12(3):513±22.

K.C. Abbaspour et al. / Advances in Water Resources 24 (2001) 827±841[19] Murtagh BA, Saunders MA. MINOS/AUGMENTED user's
manual. Tech. Rep. 80-4, System Optimization Laboratory,
Department of Operational Research, Stanford University, Stanford CA, 1980.
[20] Nakamura M, Kurumatani K. A mathematical model for theforaging of an ant colony and pattern formation of pheromone
trail. In: Fundamental theories of deterministic and stochastic
models in mathematical biology, 1995; Institute of Statistical
Mathematics. p. 120±31.
[21] NRC, Groundwater models, scienti®c and regulatory applications. Committee on Groundwater Modeling Assessment, Water
Science and Technology Board, National Research Council,
National Academy Press, Washington, DC, 1990.
[22] Romano N, Santini A. Determining soil hydraulic functions from
evaporation experiments by a parameter estimation approach:experimental veri®cations and numerical studies. Water Resour
Res 1999;35(11):3343±59.
[23] Segol G. Classic groundwater simulations; proving and improving
numerical models. Englewood Clis, NJ: Prentice-Hall; 1994.
[24] Sheble GB, Brittig K. Re®ned genetic algorithm ± economic
dispatch example. IEEE Trans Power Systems 1995;10(1):117±23.
[25] Simunek J, van Genuchten MTh. Parameter estimation of soil
hydraulic properties from multiple tension disc in®ltrometer data.
Soil Sci 1997;162:383±98.
[26] Simunek J, Huang K, Sejna M, van Genuchten MTh. The
HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably
saturated media, version 1.0. IGWMC ± TPS ± 70. International

841

Groundwater Modeling Center, Colorado School of Mines,
Golden, Colorado, 1998.
[27] Simunek J, Sejna M, van Genuchten MTh. The HYDRUS-2D
software package for simulating two-dimensional movement ofwater, heat, and multiple solutes in variably saturated media,
version 2.0. IGWMC ± TPS ± 53. International Groundwater
Modeling Center, Colorado School of Mines, Golden, Colorado,1999.
[28] Sonnleitner M, Gunthardt-Goerg M, Bucher-Wallin I, Attinger
W, Reis S, Schulin R. In¯uence of soil type on the eect ofelevated atmospheric CO2 and N deposition on the water balanceand growth of a young spruce and beech forest. Water, Air, Soil
Pollut 2001;126:271±90.
[29] van Genuchten MT. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci Soc Am J
1980;44:892±8.
[30] van Genuchten MT. Non-equilibrium transport parameters frommiscible displacement experiments. Research Report No. 119, US
Salinity Laboratory, Riverside, CA, 1981.
[31] van Genuchten MT, Sudicky EA. Recent advances in vadose zone
¯ow and transport modeling. In: Parlange MB, Hopmans JW,
editors. Vadose zone hydrology cutting across disciplines. New
York: Oxford University Press; 1999. p. 155±93.
[32] van Genuchten MTh, Leij FJ, Yates SR. The RETC code forquantifying the hydraulic function of unsaturated soils. EPA/600/
2-91-065, US EPA, Oce of Research and Development, Washington, DC, 1991.