SSSAJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 25 January 2008
Published in Soil Sci Soc Am J 72:305-319 (2008)
DOI: 10.2136/sssaj2007.0176
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
Related Collections
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Water Flow Models
Right arrow Volcanic Soils

SOIL PHYSICS

Comparison of Three Multiobjective Optimization Algorithms for Inverse Modeling of Vadose Zone Hydraulic Properties

Thomas Wöhlinga,*, Jasper A. Vrugtb and Gregory F. Barklec

a Lincoln Environmental Research, Lincoln Ventures Ltd., Ruakura Research Centre, Hamilton, New Zealand
b Center for Nonlinear Studies (CNLS), Los Alamos National Lab., Los Alamos, NM 87545
c Aqualinc Research Ltd., P.O. Box 14-041, Enderley, Hamilton, New Zealand

* Corresponding author (woehling{at}lvlham.lincoln.ac.nz).


   ABSTRACT
TOP
 NOTES
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Inverse modeling has become increasingly popular for estimating effective hydraulic properties across a range of spatial scales. In recent years, many different algorithms have been developed to solve complex multiobjective optimization problems. In this study, we compared the efficiency of the Nondominated Sorting Genetic Algorithm (NSGA-II), the Multiobjective Shuffled Complex Evolution Metropolis algorithm (MOSCEM-UA), and AMALGAM, a multialgorithm genetically adaptive search method for multiobjective estimation of soil hydraulic parameters. In our analyses, we implemented the HYDRUS-1D model and used observed pressure head data at three different depths from the Spydia experimental field site in New Zealand. Our optimization problem was posed in a multiobjective context by simultaneously using three complementary RMSE criteria at each depth. We analyzed the trade-off between these criteria and the adherent Pareto uncertainty. The results demonstrate that all three algorithms were able to find a good approximation of the Pareto set of solutions, but differed in the rate of convergence to this distribution. Small differences in performance of the various algorithms were observed because of the relative high dimension of the optimization problem in combination with the presence of multiple local optimal solutions within the three-objective search space. The Pareto parameter sets yielded satisfactory results when simulating the transient tensiometric pressure at predetermined observation points in the investigated vadose zone profile. The overall best parameter set was found by AMALGAM with RMSE values of 0.14, 0.11, and 0.17 m at the 0.4-, 1.0-, and 2.6-m depths, respectively. In contrast, the fit errors were substantially higher at these respective depths, with RMSE values ranging from 0.87 to 1.49 m, when using soil hydraulic parameters derived from laboratory analysis of small vadose zone cores.

Abbreviations: AMALGAM, a multialgorithm genetically adaptive search method • MOSCEM-UA, Multiobjective Shuffled Complex Evolution Metropolis algorithm • MVG, Mualem–van Genuchten • NSGA-II, Nondominated Sorting Genetic Algorithm • OI, Oruanui ignimbrite • P1 and P2, Palaeosols • SCE-UA, Shuffled Complex Evolution algorithm • TI, Taupo ignimbrite


   INTRODUCTION
TOP
 NOTES
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Water transport and the fate of solutes in the vadose zone are among the important processes soil scientists have been trying to quantify for a long time to address emerging environmental problems. When investigating these processes with conceptual or deterministic models, we have to provide information about the hydraulic properties of the vadose zone under investigation. One approach to obtain this information is the sampling of undisturbed vadose zone cores that are subsequently analyzed in the laboratory. Unfortunately, it has been a major challenge to integrate these small-scale measurements of the soil hydraulic properties into hydrologic models that apply across a range of spatial scales. Because of the high nonlinearity of soil hydraulic functions, their application across spatial scales is inherently problematic. Specifically, the averaging of processes determined from discrete small-scale samples may not be representative of the key hydrologic processes of the larger spatial domain. In addition, the dominant hydrologic flow processes may vary between spatial scales, so that potentially different models need to be used to describe water flow at the soil pedon, field, or watershed scale (Vrugt et al., 2004).

A promising approach to obtain "effective" hydraulic properties is the inversion of the problem, i.e., the search of effective parameters that yield the best attainable fit between model predictions and observations at the scale of interest (e.g., Inoue et al., 1998; Simunek et al., 1998; Zurmühl and Durner, 1998; Si and Kachanoski, 2000; Vrugt et al., 2001; Abbasi et al., 2003; Hupet et al., 2003; Roulier and Jarvis, 2003; Abbaspour et al., 2004; Ritter et al., 2004, 2005; Zhang et al., 2004; Kelleners et al., 2005). The application of inverse modeling to estimate vadose zone hydraulic properties across spatial scales is very promising, yielding parameters that represent effective conceptual representations of spatially and temporally heterogeneous watershed properties at the scale of interest. Physically based flow and transport models with complex multidimensional governing equations typically require small grid resolutions to minimize numerical errors, however, thereby requiring significant computational resources for simulation, and efficient optimization algorithms for model calibration.

Because of the subjectivity and time-consuming nature of manual trial-and-error parameter estimation, there has been a great deal of research into the development of automatic methods for model calibration. Automatic methods seek to take advantage of the speed and power of computers, while being objective and easier to implement than manual methods. These algorithms may be classified as local search methodologies when seeking for systematic improvement of the objective function using an iterative search starting from a single arbitrary initial point in the parameter space, or as global search methods in which multiple concurrent searches from different starting points are conducted within the parameter space.

In the field of vadose zone hydrology, nonlinear gradient-based search algorithms such as the Levenberg–Marquardt algorithm (Marquardt, 1963) have found widespread use and implementation to solve parameter estimation problems. Such derivative-based search methods will evolve toward the global optimum in the search space in situations where the objective function exhibits a convex response surface in the entire parameter domain. Unfortunately, numerous contributions to the hydrologic literature have demonstrated that the response surface seldom satisfies these conditions, but instead exhibits multiple optima in the parameter space with both small and large domains of attraction, discontinuous first derivatives, and curved multidimensional ridges (e.g., Schwefel, 1995). Local gradient-based search algorithms are not designed to handle these peculiarities, and therefore they often prematurely terminate their search, with their final solution essentially being dependent on the starting point in the parameter space. Another problem is that many of the hydraulic parameters typically demonstrate significant interaction because of limited information content of the observational data, further lowering the chance of finding a single unique solution with local search methodologies. These considerations inspired researchers in many fields of study to develop efficient global optimization methods that use multiple concurrent searches from different starting points to reduce the chance of getting stuck in a single basin of attraction. Global optimization methods that have been used for the estimation of the unsaturated soil hydraulic properties include the annealing simplex method (Pan and Wu, 1998), genetic algorithms (Takeshita and Kohno, 1999; Vrugt et al., 2001), multilevel grid sampling strategies (Abbaspour et al., 2001), ant-colony optimization (Abbaspour et al., 2001), and shuffled complex methods (Vrugt and Bouten, 2002; Mertens et al., 2005, 2006).

Despite this progress made, the current generation of automatic frameworks for estimation of the unsaturated soil hydraulic parameters typically implements a single objective function for parameter estimation. This approach is based on classical statistical theory, and relies on the assumption that the underlying simulation model is a correct description of reality, and that the input and output data are observed without error, effectively neglecting uncertainty. These assumptions can be contested in subsurface flow and transport modeling. One response to this problem, highlighted and used in this study, is to pose the optimization problem in a multiobjective context (Gupta et al., 1998; Boyle et al., 2000; Madsen, 2000, 2003; Deb et al., 2002; Vrugt et al., 2003a). By simultaneously using a number of complementary criteria and analyzing the trade-offs between the fitting of these criteria, the modeler is able to better understand the limitations of model structures and gain insights into possible model improvements (Gupta et al., 1998; Vrugt et al., 2003a; Tang et al., 2006).

In this study, we set up a multiobjective optimization framework to inversely estimate the hydraulic parameters of a volcanic vadose zone using observed tensiometer data from three different depths at the Spydia field site in New Zealand. We used HYDRUS-1D (Simunek et al., 2005) for flow modeling, and defined a RMSE value at three depths to separately measure the ability of the model to simulate the observed tensiometric data at each of these locations. The resulting optimization problem was solved with three different multiobjective algorithms. Specifically, we compared the performance and efficiency of the Nondominated Sorting Genetic Algorithm (NSGA-II, Deb et al., 2002), the Multiobjective Shuffled Complex Evolution Metropolis algorithm (MOSCEM-UA, Vrugt et al., 2003b), and AMALGAM (Vrugt and Robinson, 2007), a recently developed multialgorithm genetically adaptive search method for estimating the Pareto set of solutions. We analyzed the trade-off between the three different RMSE objective functions, and studied the resulting parameter uncertainty by interpreting the Pareto distribution as well as the associated HYDRUS-1D model prediction uncertainty. To illustrate the advantages of our multiobjective framework, we also include a comparison against laboratory-measured values for the hydraulic properties.


   MATERIALS AND METHODS
TOP
 NOTES
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Field Experiments
The Spydia site is located in the Tutaeuaua subcatchment (Landcorp's Waihora Station, E 175.79977, S 38.61423) north of Lake Taupo, New Zealand, on a sheep and beef farm under pastoral land use. The vadose zone materials at Spydia encompass a young volcanic soil (0–1.6-m depth), which belongs to the Oruanui loamy sand series (ashy/pumiceous, mesic Andic Haplorthod or Podzolic Orthic Pumice soil) developed on the underlying unwelded Taupo ignimbrite (TI, 1.6–4.4 m) from the ~1800 yr BP Taupo eruption (Rijkse, 2005). Two older buried soils (Paleosols, P1 and P2), presumably late Pleistocene to Holocene in age, follow at about 4.5- to 5.8-m depth. Their parent material was inferred to be weathered tephra or tephric loess. Below, Oruanui ignimbrite (OI) material is found, which derived from an eruption about 26,500 yr BP (Wilson, 2004).

Altogether, 100 undisturbed soil cores of about 350 cm3 were taken at 10 nominal depths down to 7.0 m below the ground surface (five vertical and five horizontal replicates at each depth). The samples were analyzed in the laboratory using standard soil physical experiments (Blakemore et al., 1987; Klute and Dirksen, 1986) to obtain the saturated and near-saturated hydraulic conductivity and the water retention function. The physical properties and textural data of the materials from the 10 sampled depths are summarized in Table 1 . The gravel fraction, i.e., the portion with particle size >0.002 m, is quite large in the TI materials (up to 37%), whereas this fraction is negligible in the Paleosols. The Paleosols have significantly different hydraulic properties compared with the other materials in the Spydia vadose zone, exhibiting a particular large total porosity, {eta}, and a much larger water holding capacity (Table 1). In this study, we focused on the unsaturated zone (0–4.2 m below ground surface), where the water moves predominantly in the vertical direction. In all our calculations, we did not include the vadose zone below 4.2 m because field observations of water level and tensiometric pressure below these depths suggested that this domain is influenced by lateral groundwater flow during at least some parts of the year.


View this table:
[in this window]
[in a new window]

  		Table 1. Average physical characteristics and textural data for the Spydia vadose zone materials as derived from laboratory analysis of core samples. The gravel fraction includes all particles >0.002 m; BD denotes dry bulk density, {eta} is the total porosity, {theta} is the volumetric water content, Ks,median is the median saturated hydraulic conductivity, and K–0.4,median is the median unsaturated hydraulic conductivity at –0.4 kPa pressure head. Variations in sampling depth for the lower Taupo ignimbrite and the Paleosol P1 are related to variations of the stratum boundaries at these depths.

 
The saturated hydraulic conductivity, Ks, was determined from the vadose zone samples for both horizontally and vertically extracted cores at each depth. The Ks values derived from the horizontal samples did not differ systematically from values derived for the vertical samples at these individual depths; however, the Ks values varied significantly with 8.5 to 71, 31 to 77, 19 to 42, 2.6 to 80, and 1.9 to 20 µm s–1 at the various sampling depths relevant for the model application considered in this study. Table 1 summarizes the median saturated conductivity, Ks,median, at the various sampling depths. Notice that the unsaturated conductivity at a tension of 0.4 kPa, K–0.4,median, determined from the vadose zone samples is about 0.1 to 2.1 orders of magnitude smaller than the Ks,median for the various depths. The hydraulic conductivity dropped only 0.1 to 0.2 orders of magnitude in the cores from the upper soil (0.0–0.6 m). On the contrary, the reduction in hydraulic conductivity was largest in the upper TI materials (1.8- and 3.5-m depths) and in the Paleosols (Table 1).

Soil water retention data was derived from a subsample of five cores per depth at tensions of 5, 10, 20, 40, 100, 300, and 1500 kPa using the hanging water column method and pressure plates. The water contents obtained at the individual tensions varied by 5 and 15% (v/v) between the replicates. Water retention curves were fitted to the data at each depth using the program SHYPFIT (Durner, 1995), which is based on a nonlinear optimization method called the modified "golden search" algorithm (Press et al., 1992). We used the van Genuchten model (van Genuchten, 1980) to portray the water retention functions. The saturated water content, {theta}s, is often approximated by the (calculated) total pore space. It is typical, however, that only a portion of the total pore space contributes to significant water movement because some pores are not connected or are blocked. Therefore the "effective" {theta}s is typically smaller than the total porosity. For this reason, both the saturated and the residual water contents were optimized (using the SHYPFIT feature) during the fitting procedure. The resulting Mualem–van Genuchten (MVG) parameters for the depths down to 4.2 m below ground surface are listed in Table 2 .


View this table:
[in this window]
[in a new window]

  		Table 2. Mualem–van Genuchten parameters of the hydraulic properties as derived from laboratory analysis of vadose zone samples: {theta}r and {theta}s denote the residual and saturated water content, respectively; {alpha} and n are shape parameters of the water retention function.

 
To measure pressure head in the vadose zone, 15 tensiometer probes (Type T4e, UMS, Munich, Germany, accuracy ±0.5 kPa) were installed outward from a central caisson with 2.3-m diameter at five depths (0.4, 1.0, 2.6, 4.2, and 5.1 m, with three replicates at each depth). The pressure head was continuously recorded at 15-min intervals since March 2006 using a compact FieldPoint controller (cFP2010, National Instruments Corp., Austin, TX) programmed for daily remote data transfer.

Model Setup
For our inverse modeling exercise, we used the HYDRUS-1D model (Simunek et al., 2005). This model simulates water flow in variably saturated porous media and utilizes the Galerkin finite element method based on the mass conservative iterative scheme proposed by Celia et al. (1990). The model solves the one-dimensional Richards equation:

Formula 1[1]
where {theta} is the volumetric water content [L3 L–3], t represents time [T], z is the vertical coordinate (positive upward) [L], h denotes the pressure head [L], K is the unsaturated hydraulic conductivity function [L T–1], and S is a sink term representing processes such as plant water uptake [L3 L–3 T–1]. The soil hydraulic properties are described by the modified MVG model (Vogel et al., 2001):

Formula 2[2]

Formula 3[3]
where Se is the effective water content, {theta}r and {theta}s denote the residual and saturated water content, respectively [L3 L–3], {alpha} [L–1] and n (unitless) are parameters that define the shape of the water retention function, Ks represents the saturated hydraulic conductivity [L T–1], l is the pore-connectivity parameter of Mualem (1976), and hs < 0 is the maximum pressure head allowed at the soil surface, i.e., the non-zero minimum capillary height (Vogel et al., 2001). In this study we assume that m = 1 – 1/n and n > 1.

The initial and boundary conditions used to solve Eq. [1] are

Formula 4[4]

Formula 5[5]
and

Formula 6[6]
where hi(z) is the initial pressure head derived from linear interpolation of observed tensions at the 0.4-, 1.0-, 2.6-, and 4.2-m depths, hL(t) denotes the prescribed (observed) pressure head at the bottom boundary L = –4.2 m (depth of the model is 4.2 m), q0(t) is the net infiltration rate (i.e., precipitation minus evaporation), and hA signifies the minimum pressure head allowed at the soil surface. Equation [6] describes the atmospheric boundary condition at the soil–air interface (Simunek et al., 1996), which switches between a prescribed flux condition and a prescribed head condition, depending on the prevailing transient pressure head conditions near the surface. Because our study considered a relatively coarse-textured soil with large infiltration capacity, we neglected infiltration-excess overland flow and used the limits hA = –200 m and hs = –0.02 m.

Simulations were set up for a period of 282 d (11 Apr. 2006–18 Jan. 2007). The initial pressure heads at 11 Apr. 2006 were –0.41, –1.38, –1.18, and –0.85 m at the 0.4-, 1.0-, 2.6-, and 4.2-m depths, respectively. Daily values of potential evaporation were calculated by the Penman–Monteith equation (Allen et al., 1998) using data from the nearby Waihora meteorological station (500-m distance). Precipitation was recorded event based on site using a 0.2-mm bucket gauge and upscaled to hourly values for use in our calculations. The plant water uptake, S in Eq. [1], was simulated by the model of Feddes et al. (1978) using HYDRUS-1D default parameters for grass and a rooting depth of 0.35 m.

The HYDRUS-1D model was set up initially for five horizons corresponding to the different materials identified in the upper vadose zone profile (0–4.2-m depths). To reduce the number of parameters in our multiobjective parameter estimation framework, we approximated the vadose zone with three different horizons. This was achieved by combining the first three horizons of the more recent materials (0–0.69-m depths), and the other two simulation layers being the disturbed Taupo ignimbrite (0.69–1.6 m) and the in situ Taupo ignimbrite (1.6–4.2 m), respectively. Before our optimization runs, we conducted a sensitivity analysis of the HYDRUS-1D spatial discretization using various values for the nodal distance. Our numerical experiments demonstrated that a uniform distance of {Delta}x = 0.02 m (with a total of 211 nodes in the HYDRUS-1D model) was numerically stable for physically realistic parameter combinations, provided accurate simulations of the flow regime and pressure heads, and resulted in a considerable time savings for forward simulations compared with much denser spatial discretizations.

Fitting Criteria
Three different criteria were used to measure the agreement between observed and simulated tensiometric data: the RMSE (e.g., Hall, 2001), the coefficient of determination (R2), and the coefficient of efficiency of Nash–Sutcliffe (American Society of Civil Engineers, 1993):

Formula 7[7]
where hp and ho denote the simulated and observed pressure head and N is the number of data points used in the evaluation. The Ce is a widely used fitting criterion and analagous to the coefficient of determination as used in the analysis of variance. The Ce may resolve negative numbers if the mean square error exceeds the variance of the observations (Hall, 2001). Model predictions are considered satisfactory if the values of both R2 and Ce are close to unity.

Inverse Modeling Procedure
The aim of the inverse modeling procedure is to find values for the model parameters that provide the best attainable fit between model predictions and corresponding observations. In this study, we considered a multiobjective framework with three different criteria:

Formula 8[8]
where F1 to F3 are defined by the RMSE of the fit between the simulated and observed pressure heads at the depths of 0.4, 1.0, and 2.6 m in the vadose zone profile, and u is a vector of k model parameters to be optimized. In this instance, they are the MVG parameters for each of the layers in the HYDRUS-1D model. Previous studies have shown that of the six MVG model parameters, {theta}r is typically the least sensitive to the calibration data (Inoue et al., 1998; Zurmühl and Durner, 1998; Simunek et al., 1998; Kelleners et al., 2005; Mertens et al., 2006). To reduce the number of parameters to be optimized and thus reduce the computational requirements of the optimization, {theta}r was set to zero. Hence, five MVG model parameters were estimated in each of the three layers, namely {theta}s, Ks, {alpha}, n, and l, resulting in a 15-dimensional optimization problem.

To solve the multiobjective framework expressed in Eq. [8], we used three multiobjective optimization algorithms, namely the NSGA-II, the MOSCEM-UA, and the AMALGAM algorithm. Multiobjective optimization does not result in a single unique set of parameters but consists of a Pareto set of solutions (Gupta et al., 1998; Madsen, 2003; Vrugt et al., 2003a). Pareto efficient or Pareto optimal parameter sets represent trade-offs among the different objectives having the property that moving from one solution to another results in the improvement of one objective while causing deterioration in one or more others (Gupta et al., 1998; Deb, 2001; Vrugt et al., 2003a). Ideally, the multiobjective optimization algorithm should find all Pareto optimal solutions. The trade-off between the different objectives is analyzed by mapping the approximated Pareto set from the parameter (decision) space to the objective space.

To contrast the results of the different multiobjective search algorithms, we isolated four different Pareto points from the analyses that we believe are essential and most informative about the efficiency of the individual algorithms. The first three Pareto points are the best found individual solutions with respect to each of the individual objectives (subsequently referred to as Pareto extremes RMSE1, RMSE2, and RMSE3). The fourth solution is the minimum overall RMSE, i.e., where the average RMSE of all three objectives is at its minimum (RMSE0). This point was subjectively chosen and subsequently referred to as the compromise solution.

NSGA-II
The Nondominated Sorting Genetic Algorithm NSGA-II (Deb et al., 2002) was derived from its parent algorithm NSGA (Srinivas and Deb, 1995). The method uses a genetic algorithm for population evolution in combination with a fast nondominated sorting approach to classify solutions according to level of nondomination and a crowding distance operator to preserve solution diversity. The NSGA-II algorithm has been successfully used for multiobjective optimization in a number of fields, such as mechanical and chemical engineering and bioinformatics (e.g., Deb and Reddy, 2003; Salazar et al., 2006; Zhang et al., 2006; Gupta et al., 2007; Khosla et al., 2007). The basic steps of the algorithm can be summarized as follows: An initial population (a number of MVG parameter sets) of size s is generated, and this population is sorted based on the nondomination level of the individuals. Each individual (single parameter set) is then assigned a rank (fitness) value and the crowding distance, a measure of how close an individual is to its neighbors, is computed. Individuals are selected from the population by using a binary tournament selection, which is based on the rank and the crowding distance. The selected population generates children (offspring) using traditional evolutionary operators (crossover and mutation). The resulting population comprising the current population and offspring is sorted again based on nondomination. Then the best s individuals constitute the new population. The process then repeats to generate subsequent generations. An extensive description of the NSGA-II algorithm, including a detailed explanation of tournament selection, crossover, mutation, nondomination, and crowding distance appears in Deb et al. (2002) and Seshadri (2006), and so is not repeated here.

An advantage of NSGA-II is that the method requires minimum user input. Apart from the number of objectives (M = 3) and the number of parameters (k = 15), only two algorithmic parameters have to be defined by the user: s and the number of generations (gen). In this study, we used s = 200 and gen = 100 to obtain a total of 20,000 model evaluations. The feasible parameter space was constrained to the bounds given in Table 3 .


View this table:
[in this window]
[in a new window]

  		Table 3. Lower and upper parameter bounds used in the NSGA-II, MOSCEM-UA, and AMALGAM multiobjective optimization algorithms.

 
MOSCEM-UA
The Multiobjective Shuffled Complex Evolution Metropolis algorithm (MOSCEM-UA; Vrugt et al., 2003a) was developed for solving the multiobjective optimization problem of hydrologic models. Various case studies have demonstrated that the algorithm is capable of generating a fairly uniform approximation of the "true" Pareto front within a single optimization run. It combines a Markov Chain Monte Carlo sampler with the Shuffled Complex Evolutionary algorithm (SCE-UA; Duan et al., 1992) and the probabilistic covariance-annealing process of the Shuffled Complex Evolution Metropolis algorithm (SCEM-UA; Vrugt et al., 2003b). MOSCEM-UA uses the concept of Pareto dominance to evolve the initial population of points toward a set of solutions stemming from a stable distribution. The algorithm can be summarized in three steps: First, the multiobjective vector F(u) is calculated for each point using a uniformly distributed initial population (of parameter sets) of size s. After sorting the population by fitness value, it is divided into a predefined number of complexes and within each complex a Markov chain (sequence) is initialized. Subsequently, in the second step, these sequences are evolved using their current location, and the covariance induced between the various individuals in the complex. The third step comprises the reshuffling of the complexes. Finally, the members of the complexes are combined in the "new" population and sorted in order of decreasing fitness value. If the convergence criterion is satisfied, the algorithm stops, otherwise it returns to the second step. For more details on the algorithm, refer to Vrugt et al. (2003a). Based on the findings in Vrugt et al. (2003a) and Tang et al. (2006), the MOSCEM-UA performance is most sensitive to three algorithmic parameters: the population size (s), the maximum (total) number of function evaluations (ndraw), and the number of complexes or sequences (q). For application in this study, we used s = 200, ndraw = 20,000, and for computational efficiency we decided to run MOSCEM-UA with a relatively low number of complexes q = 4.

AMALGAM
AMALGAM has recently been presented in Vrugt and Robinson (2007), and is a novel search algorithm that combines two concepts, simultaneous multimethod search and self-adaptive offspring creation, to ensure a reliable and computationally efficient solution to multiobjective optimization problems. The algorithm implements a population-based elitism search procedure to find a well-distributed set of Pareto solutions within a single optimization run. AMALGAM is initiated using a random initial population of size s, generated using Latin hypercube sampling. Then, each parent from the initial population is assigned a rank using the fast nondominated sorting algorithm (Deb et al., 2002) as used in the NSGA-II algorithm. The population of offspring (with size s) is subsequently created using the multimethod concept: instead of utilizing a single biological or physical model for reproduction, a diverse set of operators is run simultaneously. The individual operators are favored adaptively accordingly to their reproductive success during the search. The offspring and parent population is ranked together and members for the next population are chosen from subsequent nondominated fronts based on their rank and crowding distance (Deb et al., 2002). The new population is then used to generate offspring and the aforementioned algorithmic steps are repeated until convergence is achieved. A detailed description of AMALGAM appears in Vrugt and Robinson (2007), and so is not repeated here. Note that the adaptive search properties of AMALGAM reduce the need for tuning of the algorithmic parameters. In fact, the main algorithmic parameter of AMALGAM is its population size. Various synthetic multiobjective benchmark studies have demonstrated that the method works well with a population size of about 50. To enable a fair comparison between the different algorithms, however, we implemented AMALGAM in this study with s = 100.

Determination of Single-Objective Solutions
We used the SCE-UA global optimization algorithm developed by Duan et al. (1992) to calculate the best attainable parameter values ui,opt for each of the objectives Fi. The SCE-UA algorithm has been reported to be consistent, effective, and efficient for watershed model calibration (e.g., Kuczera, 1997; Boyle et al., 2000; Lin and Radcliffe, 2006) but also for parameter estimation of other models including soil hydrologic models (Scott et al., 2000; Mertens et al., 2004, 2006; Yang et al., 2005). The method is based on a synthesis of the best features of several different optimization algorithms, including controlled random search, competitive evolution, genetic mutation, and complex shuffling. The effectiveness and the efficiency of the SCE-UA algorithm are sensitive to the choice of the algorithmic parameters (Duan et al., 1994). The most important parameter is the number of complexes, p. Its optimal value is strongly dependent on the nature and dimensionality of the problem under investigation: the higher the degree of difficulty (rough response surfaces or spaces with many local minima), the larger the number of complexes required to locate the global optimum. We initially analyzed the performance of the SCE-UA method for our problem with different values of p and found it to be most efficient using p = 17. All other algorithmic parameters were chosen accordingly to the recommendations of Duan et al. (1994). The parameter space was constrained by the bounds specified in Table 3.

Due to the global nature of the search, it is possible that many physically unrealistic parameter combinations are generated, which typically result in nonconvergent numerical solutions of the HYDRUS-1D model. To prevent the optimization algorithms from further exploring solutions in this space, we assigned a large RMSE value to these parameter combinations in all our optimization runs. In each optimization trial, a total of 20,000 HYDRUS-1D model evaluations were conducted, after which the results presented here were reported. In addition, the first 10 d of the simulation were not used in our analysis and the computation of our objective functions to reduce sensitivity to state value initialization.

Initial Sampling Distribution
Typically, in subsurface flow and transport modeling, we have very limited information about the location of the Pareto distribution in the parameter space before any experimental data is collected. A common approach is therefore to assign upper and lower bounds for each of the parameters, and to use uniform sampling in this space to create the initial population of points to be iteratively improved with the optimization algorithm. This sampling strategy is referred to as initial sampling without prior knowledge and was used in this study. An alternative implementation, tested by us in preliminarily to this study, was introduced by Vrugt et al. (2003a). Prior information from the single-criterion solutions (i.e., the Pareto extremes) of the individual objectives was used for the construction of the initial sample. The preliminary study was conducted for both the NSGA-II and the MOSCEM-UA algorithms. We found that using prior information for generating the initial population of parameter sets appeared to result in quicker convergence of the algorithms. The overhead in generating this information, however, required an unacceptably large number of HYDRUS-1D model evaluations (60,000 runs) due to the slow convergence of the single-objective SCE-UA algorithm optimizations. In light of these considerations, we therefore only report the outcomes of our optimization runs using uniform prior distributions for each of the individual parameters.


   RESULTS AND DISCUSSION
TOP
 NOTES
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Laboratory Parameter Estimation
The parameters of the five retention functions (Table 2) and an assumed pore-connectivity parameter of l = 0.5 were used to parameterize the five-layer HYDRUS-1D model. Three simulation runs were conducted using minimum, median, and maximum Ks values from the vadose zone sample analysis. Figure 1 shows the simulated pressure head at 0.4-, 1.0-, and 2.6-m depths in comparison to the observed data. It is evident from Fig. 1 that simulated tensions do not match the observations for any of the runs or the three depths. The poor model fit with median Ks values is confirmed by relatively large RMSE values of 1.31, 1.49, and 0.87 m for the 0.4-, 1.0-, and 2.6-m depths, respectively. Negative Ce values (c.f. Table 4 ) indicate a strong negative bias between observed and simulated pressure heads, with the discrepancy between model predictions and observations increasing with time (Fig. 1). All three simulations show the water in the vadose zone profile draining at a larger rate than observed. Since l was the only assumed MVG model parameter in these simulations, we also conducted a simulation with median Ks values and MOSCEM-UA optimized l values (l1 = 0.66, l2 = 0.13, l3 = 1.0, l4 = 0.99, and l5 = 1.0) instead of using l = 0.5 in the simulations. Notwithstanding, the RMSE between simulated and observed pressure head data at each depth remained relatively large, with values of RMSE = 1.06, 1.39, and 0.83 m, respectively (Fig. 1). Out of the four simulation runs, the fit to the observed pressure heads was best using minimum Ks values, but the coefficients of efficiency were still negative (Table 4).


Figure 1
View larger version (37K):
[in this window]
[in a new window]

  		Fig. 1. Observed and simulated pressure head at (a) 0.4-m depth, (b) 1.0-m depth, and (c) 2.6-m depth using minimum, median, and maximum values of the saturated hydraulic conductivity, Ks, derived from laboratory analysis of vadose zone samples. A simulation with median Ks and optimized pore-connectivity parameter l values is also shown.

 

View this table:
[in this window]
[in a new window]

  		Table 4. Root mean square error, coefficient of determination (R2) and coefficient of efficiency (Ce) for the fit between observed and simulated pressure heads using parameters derived from laboratory analysis and optimized parameter sets. Values are shown for the 0.4-m (L1), the 1.0-m (L2), and the 2.6-m (L3) depths.

 
It is obvious that the laboratory-derived hydraulic parameters do not yield an accurate prediction of the pressure heads at the various observation depths. These results confirm the commonly observed inability of small-scale measurements to provide meaningful estimates of the soil hydraulic properties and hydrologic flow regime across larger spatial scales. This necessitates the use of inverse modeling to back out the appropriate values of the hydraulic parameters, preferably including an estimate of their uncertainty.

Single-Criterion Solutions
We identified the single-criterion solutions of the individual objectives using the SCE-UA optimization algorithm. Convergence of this algorithm was achieved within 20,000 HYDRUS-1D model evaluations in the cases of objectives F1 and F2. The SCE-UA was unable to converge for objective F3 within a reasonable number of model evaluations, resulting in considerable uncertainty in the optimal values that minimize this particular objective. The best attainable parameter values with respect to the three individual objectives resulted in RMSE values of 0.09, 0.08, and 0.14 m for the 0.4-, 1.0-, and 2.6-m depths, respectively. Large R2 and Ce values for objectives F1 and F2 (Table 4) suggest that a good model fit of the simulated and observed pressure head can be obtained for the depths 0.4 and 1.0 m. Although the observed pressure heads are least variable at the 2.6-m depth, the model fit for the corresponding objective F3 was worse than for F1 and F2 (Table 4). We assume that this is a result of the inability of the SCE-UA algorithm to identify defined areas of attraction in the parameter space for this particular case. The SCE-UA derived parameter values for each of the three objectives are listed in Table 5 .


View this table:
[in this window]
[in a new window]

  		Table 5. Optimized parameter sets as derived from the SCE-UA algorithm for single objectives (minimum RMSE for each of the objectives F1–F3), as well as for NSGA-II, MOSCEM-UA, and AMALGAM compromise solutions (smallest overall RMSE); {theta}s denotes the saturated water content, {alpha} and n are shape parameters of the water retention function, Ks represents the saturated hydraulic conductivity, and l is a pore-connectivity parameter.

 
Multiobjective Inverse Parameter Estimation
The results of the multiobjective optimization with the NSGA-II algorithm, the MOSCEM-UA algorithm, and AMALGAM are presented in Fig. 2 to 5GoGoGo and discussed below. The parameter values of the four Pareto points analyzed and the fitting criteria measuring the agreement between observed and simulated tensiometric data are listed in Tables 4 and 5. Figure 2 depicts the bi-criterion F1–F2, F1–F3, and F2–F3 fronts of the full three-dimensional Pareto trade-off space (objective space), with the final Rank 1 Pareto points (for Pareto ranks, cf. Vrugt et al., 2003a) of the last 2000 HYDRUS-1D evaluations indicated with black circles. It seems that considerable trade-off exists between objectives F1 and F2 (Fig. 2a, 2d, and 2g). The trade-off is largest for the NSGA-II Pareto solutions (Fig. 2a) and smallest for AMALGAM (Fig. 2g), indicating that this latter algorithm finds better solutions for objective F2, exhibiting less trade-off. The NSGA-II algorithm has sampled the resulting Pareto front most densely and uniformly, whereas MOSCEM-UA produced an indistinct front exhibiting large gaps. To compare the results, it should be considered that the number of MOSCEM-UA Rank 1 Pareto points was 55. In contrast, all 200 points of the last NSGA-II generation were Pareto Rank 1. This is caused by a difference in search strategy. The AMALGAM Rank 1 Pareto points are more evenly distributed along the Pareto front than the ones of the MOSCEM-UA algorithm but slightly less than the NSGA-II points.


Figure 2
View larger version (117K):
[in this window]
[in a new window]

  		Fig. 2. Pareto optimal solutions (solid circles) of the three-dimensional Pareto trade-off space as determined by the (a–c) NSGA-II algorithm, (d–f) MOSCEM-UA, and (g–i) AMALGAM: (a, d, g) F1–F2 plane, (b, e, h) F1–F3 plane, and (c, f, i) F2–F3 plane of the objective space. The compromise solutions are separately indicated in each panel by the + symbol.

 

Figure 3
View larger version (30K):
[in this window]
[in a new window]

  		Fig. 3. Normalized range of Pareto optimal parameter values for the three vadose zone materials in the Spydia HYDRUS-1D model including the compromise solution (solid lines) and the Pareto extremes (dashed and dashed-dotted lines), i.e., the best-fit solutions for each of the individual objectives F1–F3: (a) using the NSGA-II algorithm, (b) using the MOSCEM-UA algorithm, and (c) using AMALGAM.

 

Figure 4
View larger version (45K):
[in this window]
[in a new window]

  		Fig. 4. Water retention and hydraulic conductivity functions for (a–b) the 0.4-m depth, (c–d) the 1.0-m depth, and (e–f) the 2.6-m depth shown for the compromise solutions of NSGA-II, MOSCEM-UA, and AMALGAM as well as for vadose zone cores analyzed in the laboratory.

 

Figure 5
View larger version (55K):
[in this window]
[in a new window]

  		Fig. 5. Observed and simulated pressure head using the compromise solution parameter set and the best-fit parameter sets for the objectives F1, F2, and F3 (Pareto extremes) determined by AMALGAM at: (a) the 0.4-m depth, (b) the 1.0-m depth, and (c) the 2.6-m depth.

 
In contrast to the front between objectives F1 and F2, the F1–F3 and F2–F3 fronts exhibit a strong rectangular trade-off pattern with unevenly sampled points (Fig. 2b, 2c, 2e, 2f, 2h, and 2i), demonstrating that the HYDRUS-1D model is able to minimize these objectives simultaneously using only a single combination of values for the hydraulic parameters. The reason for the curved F1–F2 Pareto front is probably found in the model structure. In our HYDRUS-1D application, we assumed that water transport occurs solely through the matrix, and neglected the presence and dynamics of preferential flow paths, which, as we believe, at least temporally develop close to the soil surface and hence may have an impact on measured pressure heads at a depth of 0.4 m (measured in objective F1). Other potential structural inadequacies in the model include the zonation of the strongly stratified vadose zone at the experimental site into only three different layers and the assumption of uniform water uptake throughout the root zone. Of course, we could have included root water uptake in our model formulation, but decided to focus in this particular study on the values of the hydraulic parameters in the various layers.

Figure 3 shows normalized parameter ranges corresponding to the Pareto solutions derived with the NSGA-II, MOSCEM-UA, and AMALGAM algorithms. The 15 HYDRUS-1D model parameters are listed along the abscissa (indices correspond to the three layers in the model), while the ordinate corresponds to the parameter values scaled according to their prior uncertainty ranges (as given in Table 3) to yield normalized ranges. Each line in the graph represents one Pareto solution. It is interesting to note that the NSGA-II Pareto solutions (Fig. 3a) are grouped closely around the Pareto extremes, i.e., the best-fit solutions for each of the individual objectives, and that they express little uncertainty around these optimal values, suggesting that they are well identified by calibration against the observed tensiometer data. Also, notice that for at least some of the parameters, significant bifurcations are apparent. The dense clustering of Pareto solutions around each of the individual objectives is not as apparent for the MOSCEM-UA algorithm (Fig. 3b) or AMALGAM (Fig. 3c). Most of the HYDRUS-1D model parameters are poorly identifiable, with Rank 1 solutions found across a large portion of the prior defined parameter space. Only n1 and n2 appear to be identifiable, and well constrained with Pareto optimal values that cover a small portion interior to their initial range.

The Pareto extremes derived with the NSGA-II, MOSCEM-UA, and AMALGAM algorithms resulted in similar predictive performance and optimized parameter values as their counterparts derived separately with the single-objective SCE-UA optimization runs. This is true in the case of objectives F1 and F2 (Table 4). In the case of objective F3, all three multiobjective algorithms performed slightly better than the SCE-UA method, resulting in RMSE3 values of about 0.11 m (Table 4). The compromise solutions determined by the algorithms are indicated in Fig. 2 using the + symbol, and corresponding VGM model parameters are listed in Table 5. The best compromise solution was found with AMALGAM (RMSE0 = 0.14), followed by the MOSCEM-UA (0.16), and NSGA-II (0.17) algorithms. These differences in predictive performance appear rather small in the objective function space, but are caused by some notable differences in location of the Pareto solution set in the parameter space derived with the three algorithms. This can be seen by careful interpretation and comparison of the Pareto parameter panels, previously presented in Fig. 3.

The fact that different initial populations might result in slightly different optimized parameters was caused by the complexity of the optimization problem that we were solving. Detailed investigation of the response surface (not further illustrated here) showed the presence of a large number of local optimal solutions, probably causing premature convergence of the SCE-UA method and the MOSCEM-UA and NSGA-II methods. The presence of these local solutions is caused by the dimensionality of the problem, and perhaps more importantly the optimization of multiple objectives simultaneously that demonstrate to have significant trade-off. Another factor that might come into play is the presence of numerical errors in the transient pressure head simulations used in the computation of the objective function. These errors, albeit small, are nonstationary, and heavily depend on the parameters in the forward model. Some of these local optimal solutions might therefore be smoothed and disappear if lower tolerance values were used for the numerical errors in the Picard solution scheme in HYDRUS-1D; however, this would have caused a significant increase in the computational time needed to solve the inverse problem with the individual algorithms.

It is likely that the performance of the individual optimization algorithms would have been more similar in terms of parameter estimates and final objective function values if the MOSCEM-UA and NSGA-II would have been used in conjunction with a larger population size. Running the MOSCEM-UA and NSGA-II algorithms with a larger population size, however, would have significantly deteriorated their efficiency. Moreover, in our study we used common values of the population. Despite these findings, the strength of AMALGAM is its good and reliable performance with a relatively small population, inspiring confidence in the efficiency of this multimethod search algorithm.

Figure 4 compares the water retention and unsaturated hydraulic conductivity functions of the compromise solutions derived with the various optimization algorithms with those obtained from laboratory analysis using Ks,median. The water retention functions corresponding to the uppermost vadose zone layer (d = 0.0–0.69 m) are significantly different for each of the three compromise solutions (Fig. 4a) and are all quite different from the laboratory-measured curve. For the other two depths considered in this study and represented by F2 and F3 (Fig. 4c and 4e), the differences between the compromise solutions derived with the three different optimization algorithms appear smaller, but again show poor correspondence with independently derived retention functions from small soil cores. With the exception of the 0.4-m depth (AMALGAM), the compromise solutions for each depth suggest a much smaller water holding capacity than the laboratory analysis throughout the entire modeled vadose zone profile (Fig. 4a, 4c, and 4e).

The unsaturated hydraulic conductivity functions of the MOSCEM-UA and AMALGAM compromise solutions are similar at each depth (Fig. 4b, 4d, and 4f). The conductivity function of the NSGA-II solution again has a different shape at the 0.4-m depth but also at the 2.6-m depth. With the exception of the NSGA-II solution of objective F1, the optimized conductivity values corresponding to the compromise solutions are generally about one order of magnitude larger for all depths than their (median) laboratory estimates within the range of observed pressure heads in the field. On the contrary, optimized K values are about two orders of magnitude smaller at a pressure head of –4 m.

The optimized retention and hydraulic conductivity functions for each of the individual optimizations are significantly different from those obtained using laboratory analysis. It is likely that the discrepancies are caused by differences in spatial support of both methods and perhaps more importantly by structural inadequacies in the model. The water transport Eq. [1] used in this study accounts for flow in the vadose zone matrix only. The optimized parameter values for {theta}s and K, however, suggest the presence of faster flow paths at the Spydia site. With the optimized parameters, the model mimics preferential flow by lesser flow paths (smaller {theta}s values) and larger hydraulic conductivity than corresponding laboratory data. Inspections of the soil surface at the Spydia site did not reveal physical soil cracking but preferential flow can occur in the present vadose zone materials in the form of gravity-driven unstable flow processes (fingered flow).

To further analyze the performance of the compromise solutions, Fig. 5a to 5c depict time series plots of HYDRUS-1D model simulations of the pressure head at 0.4-, 1.0-, and 2.6-m depths. Each gray line going from left to right across the time series plot denotes the prediction of a single AMALGAM Pareto solution, whereas the solid black line represents the observed values. The pictures are very similar for the MOSCEM-UA and NSGA-II Pareto solutions and therefore are not shown here. Notice, that the Pareto uncertainty generally brackets the tensiometric observations but seems quite large, especially in the deeper vadose zone. Also notice that the best solutions for the individual depths (dashed-dotted lines) match well with the observations (solid lines). The fit is similarly good for the compromise solution parameter set (dashed lines), as previously reported RMSE values are within close range of the optimized values for each individual objective. In fact, the compromise simulations result in somewhat lower fits to each individual objective (Table 4), but are still much better than the HYDRUS-1D model simulations obtained using the laboratory-derived parameters from the small vadose zone cores. Both the AMALGAM Pareto extreme and the compromise solution at the 0.4-m depth is found in the center of the Pareto set of solutions (Fig. 5a), whereas for the other two depths these solutions are either found at the lower end of the Pareto prediction uncertainty intervals (1.0-m depth; Fig. 5b) or at the upper bound (2.6-m depth; Fig. 5c). The uncertainty bounds are slightly smaller for the MOSCEM-UA algorithm and are more centered on the optimal solution at the 2.6-m depth (not pictured). Both findings can probably be explained by the smaller number of Rank 1 solutions obtained with the MOSCEM-UA algorithm. The simulation with the AMALGAM compromise solution parameter set resulted in a better match to the observed pressure heads at the 0.4-m depth than the simulations with the corresponding NSGA-II and MOSCEM-UA parameter sets (Table 4). This is especially true after prolonged dry periods (Fig. 5a).

Convergence of AMALGAM vs. MOSCEM and NSGA-II
For optimization runs with computationally demanding numerical models such as HYDRUS-1D, the ideal global optimization algorithm would find the Pareto set of solutions quickly, i.e., with the smallest possible number of HYDRUS-1D model evaluations. For instance, a single optimization run with 20,000 HYDRUS-1D model evaluations required between 4 and 5 d time in sequential mode under the Matlab R2007a (64 bit)–Microsoft Windows XP Professional (64 bit) modeling environment on a Dell Precision 390 workstation with a Quad-Core Intel Core2 Extreme processor QX6700 (2.67 GHz) and 2GB of RAM.

We analyzed the convergence behavior of the three algorithms in this study using the four reference points of the Pareto set. Figure 6a shows the evolution of the compromise solution, RMSE0, with an increasing number of HYDRUS-1D evaluations. AMALGAM converges fastest to a function value RMSE0 < 0.16 m at about 3000 model evaluations. It is followed by the MOSCEM-UA algorithm at about 6500 evaluations. The slowest convergence was observed for the NSGA-II algorithm, which seemed to have difficulty finding preferred solutions within 20,000 model evaluations. We obtained similar results when we used a smaller population size of s = 100, as recommended by Deb et al. (2002) in the NSGA-II optimization (not shown). In terms of the compromise solution, AMALGAM also performed best with an overall lowest RMSE0 of about 0.14 m after about 12,500 HYDRUS-1D model runs (Fig. 6a). AMALGAM also proved to be very quick in locating the various Pareto extremes (Fig. 6b–6d). This makes AMALGAM the computationally most efficient algorithm and superior to the other two algorithms in the present study. Second in the rank of convergence performance was the MOSCEM-UA algorithm, which showed quicker convergence to smaller function values than the NSGA-II algorithm.


Figure 6
View larger version (25K):
[in this window]
[in a new window]

  		Fig. 6. Evolution of the best (minimum) root mean square error (RMSE) values as a function of the number of HYDRUS-1D model evaluations with the NSGA-II and MOSCEM-UA algorithms and AMALGAM: (a) compromise solution RMSE0, (b) best fit for 0.4-m depth RMSE1, (c) best fit for 1.0-m depth RMSE2, and (d) best fit for 2.6-m depth RMSE3. Single-objective SCE-UA runs are also separately included.

 
The single-objective SCE-UA runs are also depicted in Fig. 6b to 6d. All algorithms converged to about the same RMSE value for the objectives F1 and F2 within 20,000 model evaluations; however, NSGA-II, MOSCEM-UA, and AMALGAM converged significantly faster to smaller RMSE values for objective F3 than the SCE-UA algorithm. This is an interesting finding. Hence, it might prove more efficient to run a multiobjective algorithm to find single-objective solutions. This important result deserves more attention in future work.


   CONCLUSIONS
TOP
 NOTES
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
In this study, we compared the efficiency of three different multiobjective optimization algorithms for estimation of the vadose zone hydraulic properties in the HYDRUS-1D model using observed tensiometric data from the Spydia field site in New Zealand. The three algorithms were compared with respect to computational efficiency as well as parameter and Pareto uncertainty. In our multiobjective framework, we used three different RMSE objective functions, each measuring the discrepancy between observed and simulated tensiometric pressure data at a particular depth. The algorithms investigated in this study, namely the NSGA-II, the MOSCEM-UA, and the AMALGAM algorithms, found a similar good approximation of the Pareto set of solutions. AMALGAM and NSGA-II generated a similarly uniform distribution of Pareto points along the Pareto fronts, caused by their common method of nondominated sorting. The Pareto fronts showed a relatively large trade-off between Objectives 1 and 2 (RMSE at 0.4- and 1.0-m depths) compared with the trade-off between other objectives. The sharpest Pareto front between these two objectives, i.e., the smallest trade-off, was found by the AMALGAM algorithm. Consequently, AMALGAM performed best in finding the overall best attainable parameter set, which resulted in an averaged RMSE of 0.14 m. In addition, AMALGAM was also the fastest algorithm in finding the Pareto optimal parameter sets by using a relatively small population size of 100. The second best in terms of Pareto convergence was the MOSCEM-UA algorithm, followed by NSGA-II.

Differences between the compromise solutions of the various algorithms were rather small in the objective function space, but appeared to be considerable in the parameter space for at least a few hydraulic parameters. These differences were caused by the presence of numerous local optimal solutions in the response surface, which were probably created by the relative high dimensionality of the optimization problem in combination with multiple objectives that showed significant trade-off and the presence of nonstationary, numerical errors in the transient pressure head simulations used in the computation of the objective function.

Simulations of pressure head at the observation locations using either the Pareto extremes for the various objectives or the compromise solution parameter combination generally resulted in a very good match to the measured pressure head values. Perhaps more importantly, the Pareto uncertainty generally encapsulated the observed data, but was quite large, especially at the 2.6-m depth. Nevertheless, model predictions at the various depths obtained with the Pareto solutions were considerably better than HYDRUS-1D predictions obtained with laboratory-derived values for the hydraulic parameters from small vadose zone samples. Moreover, access to the entire Pareto distribution increases insight into model structural inadequacies, and helps select a single value for the various hydraulic parameters that provide acceptable trade-off between the various objectives, if needed for predictive simulation.


   ACKNOWLEDGMENTS
 
Thomas Wöhling and Gregory F. Barkle thank the New Zealand Foundation for Research, Science and Technology (FRST) for funding this work as part of LVL's Groundwater Quality Protection Programme (8137-ASXS-LVL). Jasper A. Vrugt is supported by a J. Robert Oppenheimer Fellowship from the LANL postdoctoral program.


   NOTES
TOP
 NOTES
ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher. Permission for printing and for reprinting the material contained herein has been obtained by the publisher.

Received for publication May 15, 2007.


   REFERENCES
TOP
 NOTES
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
T. Wohling, N. Schutze, B. Heinrich, J. Simunek, and G. F. Barkle
Three-Dimensional Modeling of Multiple Automated Equilibrium Tension Lysimeters to Measure Vadose Zone Fluxes
Vadose Zone J., November 17, 2009; 8(4): 1051 - 1063.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
J. A. Vrugt, P. H. Stauffer, Th. Wohling, B. A. Robinson, and V. V. Vesselinov
Inverse Modeling of Subsurface Flow and Transport Properties: A Review with New Developments
Vadose Zone J., May 27, 2008; 7(2): 843 - 864.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Wöhling, T.
Right arrow Articles by Barkle, G. F.
Related Collections
Right arrow Inverse Procedures/Parameter Estimation
Right arrow Water Flow Models
Right arrow Volcanic Soils

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Vadose Zone Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome