Estimating Soil Hydraulic Parameters Using A Metaheuristic Algorithm Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Abstract: This paper propose a complete mechanism for assessing the movement of water in the unsaturated media. In the direct problem, the flow equation for unsaturated media based on Richards equation is numerically modeled with finite differences weight schemes. In the inverse problem the soil hydraulic parameters given by Mualem - van Genuchten model can be estimated using metaheuristic optimization algorithms. Metaheuristic algorithms have the advantage that the function chosen for optimization does not need to be continuous or differentiable because we are not using the gradient or Hessian matrix from classic optimization methods. From metaheuristic algorithms, we choose ant colony optimization because it is a new approach that has performed much better than other optimization algorithms. In the field of soil, hydrogeology and geostatistics Ant colony optimization applications based on the behavior of ants seeking an optimal path between their colony and the source of food is just in the beginnings. For this purpose the ant colony optimization is modified to estimate Mualem-van Genuchten unsaturated hydraulic soil parameters.

Abstract: Aceasta lucrare propune un mecanism complet de evaluare a curgerii apei în mediul nesaturat. In problema directa, ecuatia de curgere a lui Richards pentru mediul nesaturat este modelata numeric cu diferenÅ£e finite ponderate. In problema inversa parametrii hidraulici ai solului dati de modelul Mualem - van Genuchten pot fi estimati cu algoritmi metaeuristici. Algoritmii metaeuristici au avantajul că funcÅ£ia aleasa pentru optimizare nu trebuie să fie continua sau diferenÅ£iabilă pentru că nu folosim gradientul sau matricea hessiana din metodele clasice de optimizare. Dintre algoritmii metaeuristici s-a ales algoritmul musuroiului de furnici (ant colony optimization) deoarece cu aceasta noua abordare s-au obtinut rezultate mai bune comparativ cu alti algoritmi de optimizare. În domeniul solului, hidrogeologiei si geostatisticii aplicarea algoritmului musuroiului de furnici bazat pe comportamentul furnicilor care caută o cale optimă între colonia lor ÅŸi sursa de hrana este abia la inceput. In acest scop s-a adaptat algoritmul musuroiului de furnici pentru a estima parametrii hidraulici ai solului nesaturat dati de modelul Mualem - van Genuchten.

Keywords: modeling unsaturated zone, parameter estimation, ant colony optimization, Richards equation, finite differences

1. Introduction

The vadose zone or unsaturated zone is situated above the water table and regulates the transfer of water from the land surface to the groundwater and vice versa. The water flow in unsaturated media is a part of the hydrological cycle essentially controlling interrelationships between precipitation, infiltration, surface runoff, evapo-transpiration and groundwater recharge. During the past several decades, considerable advances have been made in the conceptual understanding and mathematical description of water flow process in variably saturated media. A large number of models of different degrees of complexity are now available for predicting subsurface flow. Still, effective application of these models to practical field situations suffers from two general problems. One issue concerns the inability of models, due to simultaneous hydrologic processes operative in the subsurface. The second issue is that even if we have a model developed that for all pertinent processes, model users would still have difficulty collecting enough meaningful field data to effectively run the model. The Richards equation is used to mathematically describe the unsaturated water flow. For the numerical solution of Richards equation we are using finite differences method (FDM). FDM approximates partial differential equation by replacing the derivatives from Richards equation with finite differences weight representations. In the vadose zone, the presence of three phases i.e., the solid phase, the liquid phase and the gaseous phase (mostly air at atmospheric pressure), make the relations between water content, matric potential and the hydraulic conductivity to be highly nonlinear. The unsaturated soil hydraulic properties, the water retention and hydraulic conductivity are taken as nonlinear function of pressure head. The soil hydraulic properties are assumed to be represented by Mualem - van Genuchten (MVG) expressions. For the determination of the unknown soil hydraulic parameters, we are using predictive methods. The use of measured primary outputs of hydrologic simulation models to estimate the model input parameters, has received special attention in the last years. Much literature is devoted to this alternative process referred as inverse modeling, model calibration, parameter fitting or parameter estimation. Reviews of the subject are given by……. Unknown parameters of the model are identified by using optimization techniques. The procedure generally involves minimization of the square difference function of some measured and simulated flow variable. Many studies investigate the best way to obtain accurate hydraulic information (for example hydraulic conductivity) of the subsurface from measurements (such as hydraulic head). Due to the limited amount of information and their associated uncertainty, these problems are often ill posed and non-unique. To minimize these shortcomings, different methods have been applied. Heuristic or approximate algorithms aim to find a good solution to a problem in a reasonable amount of computation time. A metaheuristic is a heuristic method, a high-level algorithm that combining lower-level techniques for exploration and exploitation of the search space and can be applied for solving a very general class of problems so we can apply it in the field of hydrogeology. Metaheuristics generally produce higher quality results than simple heuristics but they take a lot longer as they have to generate and evaluate many solutions rather than just one and have good result until now in difficult optimization problems. This is a new approach for modeling flow and transport in unsaturated zone. The ACO method is a metaheuristic inspired by the fact that ants are able to find the shortest route between their nest and a food source. This is accomplished by using pheromone trails as a form of indirect communication. Usage of this method in modeling unsaturated flow is just starting to develop. For the future papers the aim will be to estimate parameters in inverse problem for groundwater hydrology using more metaheuristic algorithms and to compare these algorithms and their results


The forward problem (FP) refers to mathematical model and numerical model to solve resulting partial differential equations. The direct flow problem can be formulated for any set of initial and boundary conditions, and can be solved with an analytical or numerical method.

2.1. The mathematical model of water movement in unsaturated media

The general equation of water movement in unsaturated media combined with boundary and initial conditions gives the mathematical model. Unsaturated soils distinguish from saturated soils by negative pore water pressures and saturated hydraulic conductivity, which is a constant, replacing it with unsaturated hydraulic conductivity, a nonlinear function that depends on fluid pressure or moisture content.

Combining Buckingham-Darcy's law for the water flux with the continuity equation that expresses the fundamental principle of mass conservation yields a three-dimensional, general equation of water movement in unsaturated media:


where X =(x,y,z) is the spatial variable; t is time [T]; q=q(X,t) is the specific discharge [LT-1]; H=H(X,t) is the total hydraulic head [L]; H = h - z when the origin is placed at ground level and z is directed downward; h=h(X,t) is pressure head [L]; θ=θ(X,t) is the volumetric water content [L3L-3]; z is the elevation head [L]; K=K(θ) is the unsaturated hydraulic conductivity [LT-1]; S=S(X) is the volumetric sink term representing sources of water [L3L-3T-1].

In order to obtain the general equation expressed only in function of the matric potential h, K(θ) will be replaced by K(h). Taking into account that the θ= θ(h(X,t) we obtain a general flow equation known as the Richards equation (2):


where C(h) is soil water capacity, representing the slope of the soil water retention curve.

Mualem-van Genuchten model for unsaturated soil function

Richards equation is highly nonlinear in the vadose zone because θ, C, and K vary as function of water pressure. For the parameters above, in this study, the Mualem -van Genuchten (MVG) constitutive relations are used. MVG (Mualem, 1976; van Genuchten, 1980) describe the shape of water retention curves and unsaturated hydraulic conductivity, as a function of pressure head h,


where ; θ=θ(h) - the water retention curve; θs - is saturated water content; θr -is residual water content (defined as the water content at which ); α, n, m = curve shape empirical parameters; Ks is saturated hydraulic conductivity.

1.3. Numerical solutions of unsaturated flow. The finite-difference method

The solution of unsaturated flow problems requires the use of indirect methods of analysis, based on approximations or numerical techniques. In this paper, we present various discretization schemes that use a weighted average of implicit and explicit method, for modelling unsaturated flow.

The governing equation, for unsaturated one-dimensional vertical transient flow through the rigid porous medium, is described by Richards equation:


Any initial conditions in terms of pressure head or water content can be invoked. Dirichlet or Neumann boundary conditions at the top or the bottom of the profile are associated to the partial differential equation (4):

or with z=0 or z=-L (5)

where hD(t) and qN(t) are respectively the prescribed pressure head (Dirichlet type) and the net flux (Neumann type) and L is the length of the domain.

We discretize the soil profile, the spatial domain into finite increments of Δz, and the time domain into finite increments of Δt. Denote the pressure head value at i-th depth increment and at the j-th time step by: .

The hydraulic conductivity, which is a function of the matric head, relates transfer between two layers and thus must be averaged to represent both layers:


The finite difference (FD) method approximates the derivatives in the governing flow equations by difference equations applied over small intervals. (Kumar, 2002; Taheri, Ataie, 2009)


We take


where .

Equation (9) gives the next few important cases:


If b=0 then the scheme is a fully explicit forward differences (it uses only known information from time step j to predict the unknown at time step j+1).

If b=1 then the scheme is a fully implicit backward differences (the spatial derivatives in this case are evaluated at time step j+1)

If b=1/2, the scheme is known as the Crank-Nicolson method (is the average of forward differences at time j and backward differences at j + 1)

The fully implicit method and Crank-Nicolson are unconditionally stable for any time step size.

The Richards equation after discretization, linearization and simplification can be expressed as:


The equation (10) is written for each node in which the pressure head must be calculated. The resulting linear system of N linear equations with N unknowns is a tridiagonal system and can be solved efficiently by the Thomas algorithm, a simplified form of Gaussian elimination. Input data are functions of time like measured soil-water contents, measured pressure-head profiles, under known boundary conditions and constitutive functions for the hydraulic properties are defined by equations (3). Providing adequate boundary and initial conditions, forward problem (FP) or the direct problem, solve this linear system when the model parameters are known and the state variables (e.g. pressure head at each node) can be obtained for different scenarios.


Inverse problem (IP) for determining soil hydraulic properties combine numerical solutions of nonlinear flow problems with optimization routines. IP estimate model parameters obtained from the observed data. Determining hydraulic properties of soil encompasses direct measurements and indirect estimation methods. In the field of vadose zone hydrology, inverse modeling usually involves the estimation of the soil water retention and unsaturated soil hydraulic conductivity. Solving the inverse problem consists of minimization of an objective function (minimizing the difference between observed and predicted model): the sum of squared differences between observed and computed pressure heads.

The major steps of the inverse modeling of flow parameters can be summarized as follows:

The conceptual model describing soil-water flow under a controlled transient flow experiment: the geometry of the flow model, the boundary and initial conditions are prescribed and various flow variables are measured.

The finite differences numerical method is applied to solve the Richards equation. Flow regime is simulated under test conditions, using initial estimates of the parametric soil hydraulic functions.

Selecting the parameters that should be estimated (parameters selection)

The parametric model that represents the water retention and hydraulic conductivity curves. Generally, in one-dimensional unsaturated flow problem we are using MVG model that contains five unknown parameters. These parameters are sought by the numerical inversion but depending of problem requirements, some of those parameters can be measured independently. If we assign to each N cell q parameters, vcalc will depend on qN parameters. If the number of observations is n, we define the vectors of length n, vcalc(p), (containing the flow parameters to be estimated) and vobs (containing measured pressure head and/or volumetric water content). This is a trade off in the inverse flow modeling. If too many parameters are defined to be estimated, then the flow model becomes large and unstable. On the other hand, if a small number of parameters are defined to be estimated, then the model may not be a representative one. (Farmani, 2007)

We define an objective function that describes the target variable and our goal is to make best match between the experimental data and simulated response. Objective function is defined by the difference between model predicted and measured pressure heads and/or water contents. The observation vector includes data that can be of different types, magnitudes and accuracy so each residual should be weighted according to the relationship:


where: p is the trial vector of unknown parameter values; is the measured value for v (pressure head and /or water content); is the corresponding predicted model for v (h and/or θ) using the soil hydraulic parameters of the optimized trial vector, p; n, nt are the number for all measurements of v in space and time; , where σ is the standard deviation and k is the number of samples.

Calibrating the model: The optimization procedure is a nonlinear least-squares algorithm, which can be solved by an iterative process, starting from an initial set of parameters p0. During this process, the parameters are iteratively adjusted so that the objective function is minimized and the model approximates as closely to the observed response of the system under study. Any nonlinear optimization requires the computation of a correction such as:


Calculation of the model (forward modelling); A simulation is performed using the current value of the parameter vector p to obtain the calculated system response vcalc(p). This simulation is repeated during the inverse flow modeling with updated parameters suggested by the optimization algorithm.

Optimization algorithm: Update of the parameters is made to obtain a better fit between calculated and observed data. To determine the correction , due to variability of physical parameters, we are using global search methods. Global optimization algorithms involve the evaluation of the function usually at a random sample of points in the feasible parameter space, followed by subsequent manipulations of this sample, using a combination of deterministic and probabilistic rules. Metaheuristic algorithms produce better results than heuristic algorithms combining lower-level techniques for exploration and exploitation of the search space using multiple concurrent searches from different starting points. They have the advantage of being able to be used for a wider class of problems and guarantee asymptotic convergence to the global optimum.

The inverse problem aim to find an optimal vector p* that minimizes the objective function, g(p), so that .


Heuristic or approximate algorithms aim to find a good solution (a local optimum) to a problem in a reasonable amount of computation time, with no guarantee of a global optimum solution. Metaheuristics are stochastic optimization algorithms, which employ some degree of randomness to find an optimal solution to hard problems. Metaheuristics generally produce higher quality results than simple heuristics, but they take a bigger amount of time to generate and evaluate many solutions rather than just one. Examples of metaheuristic algorithms are genetic algorithms (GA), simulated annealing (SA) , tabu search (TS), artificial neural network (ANN), evolutionary computation (EC), particle swarm optimization(PSO), ant colony optimization(ACO) and others. ACO algorithms, for most problems, are among the currently top-performing algorithms, because can run continuously, adapt to changes in real time, convergence still beeing faster than in the case of other algorithms.

3.1. General introduction to Ant Colony Optimization (ACO)

Ant Colony Optimization (ACO) is a member in swarm intelligence methods. Swarm intelligence (SI) method describes the collective HYPERLINK ""behavior of decentralized, self-organized systems, natural or artificial. ACO studies artificial systems that take inspiration from the behavior of real ant colonies and which are used to solve optimization problems. The first ACO system was introduced by Marco Dorigo in his Ph.D. thesis (1992) and was named Ant System (AS) being based on idea that the foraging ants can find the shortest path between their nest and a food source by marking the paths they follow with pheromone. This indirect communication via interaction with environment by depositing pheromones on the ground is called stigmergy and explains self-organizing capabilities of social insects. They have an autocatalytic behavior: the probability to choose a path increases with the number of times the same path had been chosen before. Each ant chooses an action, based on variable probability: first ants are forced to decide whether they should go, the choice that is made is a random decision. At a certain time, probability of choosing a branch of a path depends on the total amount of pheromone on the branch: accumulation it is faster on the shorter path and it is proportional to the number of ants that have used the branches. The process is thus characterized by a positive feedback loop that leads all the ants on a single path. Ants do not know the global structure of the problem; they discover the network systematically. (Dorigo, Socha, 2006)

3.2. Ant colony optimization algorithm

ACO is inspired by the foraging ant behavior and uses artificial ants. The real ants that want to find the food have inspired the definition of artificial ants that cooperate in finding good solutions to a given optimization problems. An artificial ant in ACO is a stochastic constructive procedure that incrementally builds a solution by adding opportunely defined solution components to a partial solution under construction. ACO is a combinatorial optimization problem of finding the best solution, with the least-cost, from all feasible solutions, that satisfy the constraints of the optimization problem. The goal is to find a global optimal feasible solution. (Dorigo, Stützle, 2004)

In ACO, artificial ants build a solution to a combinatorial optimization problem by traversing the so-called construction graph. The fully connected construction graph consists of a set of vertexes and a set of edges. The ants move from vertex to vertex along the edges of the graph, incrementally building a partial solution. Additionally, the ants deposit a certain amount of pheromone on the components that may depend on the quality of the solution found. Subsequent, the ants utilize the pheromone information as a guide towards more promising regions of the search space.

ACO algorithm:

a). Set parameters, initialize pheromone trails

while termination conditions not met do

b). Construct Ant Solutions

c). Daemon activities (Apply Local Search) {optional}

d). Update Pheromones

end while

Ant System (AS) was the first ACO algorithm proposed in the literature. Its main characteristic is that all the ants, after they have completed the tour, update the pheromone values.

AS algorithm:

a). Set parameters, initialize pheromone trails; Starting point depends on problem constraints; Initial pheromone is positive, but very small.

b). Construct Ant Solutions using pheromone trail; Each ant constructs a complete solution to the problem according to a probabilistic state transition rule; Each ant maintains a "tabu list" of vertex that it has already visited: it does not revisit any vertex.

Probabilistic selection of the edge (ij):



Ï„i,j is the amount of pheromone on a given edge (i,j) by k-th ant on arcs it has visited: ants which have detected a shorter path can deposit pheromone earlier than ants traveling on a longer path

trail visibility is the inverse distance between vertex i and j

is the neighborhood of ant k when is in node i; the neighborhood contains all the nodes directly connected to node i in the graph, except for the predecessor.

a is the importance of the pheromone intensity in the probabilistic transition

b is the importance of the visibility of the trail segment

c). Daemon activities {optional}

A local search algorithm is employed here to the constructed solutions: the locally optimized solutions are then used to decide which pheromones to update

d). Update Pheromones

Global updating rule: once all ants have built their complete tours, the pheromone is updated on all edges. The aim of the pheromone update is to increase the pheromone's values associated with good or promising solutions, and to decrease those that are associated with bad one and escape from a local minima. ACO implements a useful form of forgetting, favoring the exploration of new areas in the search space.

Using Lk, the tour length for the k-th ant and Q, a constant related to the quantity of trail laid by ants, the deposited amount (quantity of pheromone added to each edge) is :


The resulting total pheromone update function therefore is:


where ρ is the rate of pheromone evaporation ; m is the number of ants

Evaporation is applied uniformly to all edges with a simple decay coefficient ρ.

Hoos and Stützle developed MAX-MIN Ant System (MMAS) in 1996 as an improvement over the original Ant System idea: the pheromone levels are bond between a minimum and a maximum value and only the best ant can update the pheromone trails. MMAS introduced a local pheromone update after each construction step, called offline pheromone update, in addition to the pheromone update performed at the end of the construction process in AS.


Inverse methods have been applied successfully in saturated media, but estimating the parameters in vadose zone is quite difficult due to nonlinearity flow equation. Application of ACO in inverse modelling (ACO-IM) in vadose zone is an adaptation of Sequential Uncertainty FItting (SUFI), to ACO. SUFI is an iterative fitting procedure with associate uncertainty domains to each parameter, conditioning the unknown parameter estimates on an array of observed values. The stopping rule is provided by a critical value of a goal function. (Abbaspour et al, 1997, 2001)

The ACO-IM problem can be described as follows:

Define the objective function

Let continuous and bounded function define by:


Our objective in using ACO-IM is to find a parameter set, contains MVG parameters that minimize g(p).

Depict each unknown parameter by an interval based on available information; ACO has been developed for ordered problems. Its application to a continuous problem requires first discretization of the search space, where ai, bi are the lower and the upper bounds, based on a prior information using existing knowledge.

Dividing each interval [ai, bi] into a number of mi of subintervals by equal length

If each stratum is represented by the value at the middle of the stratum, then there will be permutations or possible pathways through the space of input parameters.

Run the simulation model of choice for all pathways (exhaustive sampling M), or for a randomly sampling scheme, a selected subset S of M that satisfies condition S≥0.1M so that an optimal solution will be found. To satisfy the constraint that each ant must follow only one pathway, we generate a pathway-structure list. Using this list, each ant will follow one of M, or S pathways.

Score each stratum on the bases of the smallness value of the objective function in such a way that small values of objective function receive higher scores. In ACO-IM, we use the value of the objective function of a pathway to assign the trail. Let pij , be the stratum j of parameter pi in the search space of D and be the intensity of trail on pathway u,and iteration I. The trail intensity on each path is calculated in such a way that the larger values of the objective function (i.e., longer paths for real ants) receive smaller trails. This objective is achieved through equation (16) where pathways with smaller values of the objective function receive exponentially higher trails. This is accomplished according to the following equation:


where gu is the value of the objective function for the u-th pathway, gcr is a critical value, gmin is the minimum of all gu's for an iteration

In SUFI algorithm, each stratum, initially set to a score of zero, receives a score of 1 if the value of the objective function, g, for that stratum satisfies the condition , where gcr is a user-defined value. The stratum with the largest number of hits is most likely to contain the desired estimates.

On the basis of the value of the objective function, placing a certain amount of trail (pheromone in the case of real ants) on each stratum visited along its pathway. Note that any single stratum pij in D may be the crossroad of many ant pathways. The trail share of each stratum pij from each pathway can be summed to yield Ï„ij(I) as:


Let σij (I) be the standard deviation of all values of the objective function involving pij. We are giving more weight to stratum with higher sensitivities σij, so regions of larger trail intensity are generated by smaller values of the objective function and are more likely to contain the desired optimum solution.

In order to decide how to update the parameters, will be define the term 'score', Sij , for each pij calculated by the following expression similar to the transition probability

( (19)

where μij, σij are the mean and the standard deviation of trail of the pij stratum, μg, σg are, respectively, the mean and the standard deviation of the objective function for an iteration and c1, c2 are constants experimentally found.

Eliminate from the ends of the parameter intervals the strata with small or no scores. The stratum with large sensitivities does not eliminate, maybe because they have not been represented well by their middle points.

The system is subsequently, reinitialized with updated parameters and the process repeated. The low scoring pathways will vanish (in the case of real ants the pheromone evaporates on pathways with less traffic) and new pathways evolve around high scoring paths. Repeat the process until a desired stopping rule is reached. This may occur when: a desired value of the objective function is reached, or we do not have more changes in the value of the objective function in consecutive iterations.


In this paper, we aimed to create a complex pattern of flow in unsaturated zone. It contains the numerical solution of the direct problem, using weighted finite difference method for modeling of Richards equation. FD method includes a combined version of classical discretization schemes such as Crank-Nicholson scheme or fully implicit. For a proper description of unsaturated flow, a correct description of the hydraulic functions, K(h) and θ(h), is important and necessary. MGV model was chosen for these functions because they are used in software programs that model the unsaturated media, such as Hydrus-1D, 2D. In the analysis of the unsaturated zone, one of the most challenging problems is to use inverse theory in the search for an optimal parameterization of the porous media, because direct measurements of hydraulic properties of soils are tedious and time consuming. Inverse procedures require optimization of an objective function. The objective of calibration (parameter optimization) of any model of a physical system is to identify the values of some parameters in the model, which are not known a priori. This is achieved by feeding the model with input data and comparing the computed output variables to the variable values measured in the physical system. The optimization approach is an automatic calibration of the model based on the minimization of objective function. In this study, we have used the behavior of a colony of ants to achieve this optimization. The method uses the fact that ants are capable of finding the shortest path from a food source to their nest by depositing a trail of pheromone during their walk. ACO is a parallel probabilistic algorithm, adaptable and collaborative system, based on feedback that gives high quality solutions emerging via global collaboration. We are using a modified method for vadose zone using a sequential fitting. Future work will use ACO variants as MMAS or hybrid method. We can extend the idea to the modeling of the transport of pollutants through unsaturated zones. Our target is a coupled saturated-unsaturated model of the flow and transport where we can use in parameters estimation metaheuristic algorithms. Future work will deal with the extension to 2D or 3D model and with the application on real test cases.