VZJ sign up for citetrack
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 24 August 2006
Published in Vadose Zone J 5:917-933 (2006)
DOI: 10.2136/vzj2005.0117
© 2006 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)
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 Similar articles in ISI Web of Science
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 Mertens, J.
Right arrow Articles by Barkle, G. F.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Mertens, J.
Right arrow Articles by Barkle, G. F.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Mertens, J.
Right arrow Articles by Barkle, G. F.
Related Collections
Right arrow Soil Physics
Right arrow Lysimeter/Rhizosphere Studies
Right arrow Water Flow Models

SPECIAL SECTION: PARAMETER IDENTIFICATION AND UNCERTAINTY ASSESSMENT IN THE UNSATURATED ZONE

Multiobjective Inverse Modeling for Soil Parameter Estimation and Model Verification

J. Mertensa,*, R. Stengera and G. F. Barkleb

a Lincoln Environmental Research, Ruakura Research Centre, Private Bag 3062, Hamilton, New Zealand. J. Mertens, current address: Division Soil and Water Management, Katholieke Universiteit Leuven, Kasteelpark Arenberg 20, 3001 Leuven, Belgium
b Aqualinc Research Limited, P.O. Box 14-041, Enderley, Hamilton, New Zealand
* Corresponding author (mertensja{at}yahoo.co.nz)

Received 21 September 2005.



   ABSTRACT
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Due to the exponential increase in computational power and increasing awareness of problems associated with vadose zone parameter estimations based on laboratory and in situ measurements during the last decades, the process of automatic model calibration against laboratory or field data is being increasingly used. This is often referred to as inverse modeling and has as a major limitation the inability to identify a single optimal parameter set. A coupled HYDRUS1D-SCE (shuffled complex evolution) simulation global optimization technique was developed and its suitability for multiobjective inverse modeling evaluated. In particular, the trade-off between goodness of fit against leachate volume and soil moisture content in unirrigated and irrigated lysimeters was evaluated. After identification of the most sensitive model parameters using a Monte Carlo based sensitivity analysis, the technique was capable of finding effective Pareto optimal parameter sets that well simulated leachate volume and soil moisture content in both unirrigated and irrigated lysimeters. The parameter variation along the Pareto fronts was large and differences existed between soil hydraulic parameter distributions along the Pareto fronts of the irrigated and unirrigated treatments. The multiobjective optimization technique was then adopted for the verification of the suitability of the conceptual model of equal parameter sets for both treatments. The technique was able to objectively reject the hypothesis of equal parameter sets for both treatments. This is probably due to (i) physical parameter changes with time due to the effect that long-term irrigation has on soil structure, and (ii) differences in acting transport mechanisms between the unirrigated and irrigated lysimeters.

Abbreviations: SA, sensitivity analysis • SCE, shuffled complex evolution • WCR, water content reflectometer


   INTRODUCTION
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
THE SUCCESS of model applications describing the movement of water and chemicals through the unsaturated zone relies heavily on the suitability of the conceptual model as well as the quality of the model parameters, especially the unsaturated hydraulic properties (Si and Kachanoski, 2000). Laboratory and in situ measurements of soil hydraulic properties have been common practice during the last decades (Klute, 1986); however, most methods still remain tedious, expensive, time consuming, and may include large errors (Abbaspour et al., 1999). Mallants et al. (1997) showed that hydraulic conductivity values measured in the laboratory on open-ended columns decrease with increasing column size due to the effect of the greater pore continuity at smaller scales. Similarly, as shown by Reynolds et al. (2000) for saturated conductivities, different experimental techniques may yield different parameter estimates. Jacques et al. (2002) found that alternative methods for deriving hydraulic properties from single-disk tension infiltrometer measurements resulted in different unsaturated conductivity estimates. Additionally, recent studies (Mertens et al., 2005; Oliver and Smettem, 2005) have shown that the measured soil hydraulic parameters estimated from laboratory or in situ measurements may be different from the parameters needed by a field-scale model to adequately predict the observed moisture content, pressure head, or leachate flux. For these reasons and because they are objective and reproducible, inverse modeling procedures are increasingly used to identify model parameters (Roulier and Jarvis, 2003).

Inverse modeling or automatic model calibration is defined as the process of estimating model parameters by matching a mathematical model to measured data representing the system response at discrete points in space and time (Finsterle, 2004). A large number of studies have been conducted that compare different inverse modeling algorithms (e.g., Duan et al., 1992; Gan and Biftu, 1996; Cooper et al., 1997; Kuczera, 1997; Franchini et al., 1998; Thyer et al., 1999; Madsen et al., 2002). The main conclusion from these studies is that the global population–evolution-based algorithms are more effective than multistart local search procedures, which in turn perform better than pure local search methods. The reason that local search methods have not generally given satisfactory results is that most practical problems involving calibration of nonlinear hydrological models have response surfaces that are multimodal, i.e., there are several locations in parameter space where the value of the function is a "local" minimum. In such cases, the points where a local search algorithm terminates will depend on the location where it is started (Sooroshian and Gupta, 1995). Therefore, local search optimization strategies are inappropriate for calibration of highly nonlinear hydrological models.

Traditionally, inverse modeling techniques have aimed at finding an optimal set of parameter values within some particular model structure. The limitations of the optimal parameter set concept have been discussed by Beven and Binley (1992) and Beven (1993), and Beven and Freer (2001) introduced the notion of "equifinality" to define the fact that there may be many parameter sets within a model structure that are equally acceptable in simulating the system. These may often come from very different regions in the parameter space. Given the observations available, there may be no rigorous basis for differentiating between these parameter sets. The most important implication of the "equifinality" problem is the non-uniqueness of the parameters found by the inverse modeling or calibration process. Therefore Beven (1993) proposed the Generalized Likelihood Uncertainty Estimation or GLUE methodology. The GLUE accepts that it may not be possible to distinguish between different parameter sets and aims at ranking the parameter sets based on some likelihood scale. Gupta et al. (1998) stated that the GLUE approach still has weaknesses that need to be addressed (e.g., the selection of prior parameter distributions, the likelihood criterion, and the cut-off thresholds). Gupta et al. (1998) introduced an automated procedure for multiobjective calibration. The ideas presented in Gupta et al. (1998) have similarities and notable differences from the GLUE approach. The major difference is the focus on the multiobjective nature of the model calibration problem. Application of automatic procedures has mainly been based on a single overall objective measure of comparison (e.g., root mean squared error between observed and simulated values). The use of a single measure is often inadequate to properly measure all the important characteristics of the system that are reflected in the observations used by the hydrologist to evaluate the quality of the calibrated model for the specific model application being considered (Madsen, 2000b).

In a multiobjective context, parameter combinations exist that are "equally good" due to trade-offs between the different objectives. Optimizing M objectives yields an (M – 1) dimensional surface when all of the objectives conflict with one another. It is the shape of this surface that allows the modeler to analyze tradeoffs between conflicting objectives. In the simplest case of two objectives, points on the Pareto curve have the characteristic that no other point has both a smaller value of the first objective and a smaller value of the other objective. Moving along the Pareto front results in smaller values of one objective but at the expense of the value of the second objective. It is the shape of this Pareto curve or front between two objectives that is very informative. In the simplest case of two objectives, the Pareto front is a curve on an x–y graph. A "sharp" Pareto front (Fig. 1a ) indicates that the conceptual model can simultaneously meet both objectives adequately. On the other hand, a rather "linear" Pareto front (Fig. 1b) indicates that the model is not capable of simulating both objectives simultaneously. For example, very good simulations of soil moisture content result in poor simulations of soil tension and vice versa. If this is the case, the reasons for the unsharp Pareto front must be investigated. If, in this example, it is impossible to fit soil moisture content and soil tension simultaneously, it means there is something wrong with either the measurement techniques of both objectives (e.g., measurement of soil moisture content throughout a larger volume as opposed to the small unknown volume of soil tension measurements) or the conceptual model (e.g., the assumed shape of the soil hydraulic relation between tension and moisture content) of the system. Therefore, multiobjective optimization as presented here can be used as a technique for identification and verification of the assumed conceptual model structure.


Figure 1
View larger version (5K):
[in this window]
[in a new window]
  		Fig. 1. (a) "Sharp" and (b) rather "linear" Pareto fronts. F1 and F2 are the objective functions under consideration.

 
This rationale for multiobjective "equivalence" of several parameter sets is different from the rationale of what Beven and Binley (1992) call "equifinality" of parameter sets. The arguments used here are based on the multiple ways in which the best fit of a model to the data can be defined and not on the probabilistic representation of parameter uncertainty as in the GLUE methodology. Pareto-optimal parameter sets may overlap with the behavioral parameter sets defined in a GLUE procedure but will, in general, not be equivalent. The GLUE concept assesses the parameter uncertainty from a statistical point of view, according to some likelihood measure that is, in general, based on a single aggregate objective measure.

Global optimization algorithms have been successfully applied in rainfall runoff and groundwater modeling for more than a decade (e.g., Duan et al., 1992; Gan and Biftu, 1996; Cooper et al., 1997; Kuczera, 1997; Franchini et al., 1998; Thyer et al., 1999; Madsen et al., 2002) while in vadose zone modeling studies, more traditional gradient-based methods have been common practice (e.g., Kool et al., 1987; Kool and Parker, 1988; Simunek and van Genuchten, 1997; Romano and Santini, 1999; Si and Kachanoski, 2000; Meadows et al., 2005; Kelleners et al., 2005). Only recently have vadose zone studies been published introducing global optimization algorithms in unsaturated zone modeling studies (Abbaspour et al., 2001; Lambot et al., 2002; Hupet et al., 2003; Vrugt et al., 2003a; Mertens et al., 2004, 2005; Abbaspour et al., 2004; Ritter et al., 2004, 2005). This study used the HYDRUS-1D software to simulate water transport in lysimeters with four different soil types and two treatments (irrigated or unirrigated lysimeters). Sensitivity of model parameters with respect to soil moisture content and leachate volumes was evaluated using a Monte Carlo based sensitivity analysis. Sensitive model parameters were incorporated in the inverse modeling and the model was autocalibrated on the basis of daily soil moisture content measurements and leachate volumes measured biweekly. A local Levenberg–Marquardt search procedure is readily available in the HYDRUS-1D software; however, for reasons discussed above, this study did not use the readily available local search procedure but used the SCE algorithm (Duan et al., 1992), which is a population-evolution-based global optimization method. The algorithm was implemented in MATLAB and the HYDRUS-1D software was executed from within the MATLAB environment. The Pareto front was estimated on the basis of repeated single-objective optimizations with a weighted sum of each of the objectives as the objective function. This is in contrast with the recently developed PROSCE (Pareto-based multiobjective shuffled complex evolution; Madsen and Khu, 2002) and MOSCEM (multiobjective shuffled complex evolution metropolis; Vrugt et al., 2003b) algorithms in which the multiobjective optimization is solved using a single optimization run.

The objectives of this study were to: (i) develop the joint HYDRUS1D-SCE simulation-optimization technique, (ii) use the multiobjective inverse modeling technique to investigate the trade-off in goodness of fit between leachate volumes and soil moisture content; (iii) evaluate the distributions of the Pareto-optimal parameters along the Pareto front for different soil types and two treatments (irrigated and unirrigated); and (iv) investigate the capability of the multiobjective inverse technique to verify the suitability of the conceptual model in simulating both the unirrigated and irrigated treatments simultaneously.


   MATERIALS AND METHODS
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Lysimeter Experiment Setup
The leachate volume and soil moisture content data used in this study were measured in a companion lysimeter experiment, described in detail by Barton et al. (2005) and Sparling et al. (2006). Briefly, water dynamics of grassed intact soil cores (500-mm diameter by 700-mm depth) of four soil types of vastly differing characteristics were monitored during a 4-yr period in a lysimeter facility. Half of the lysimeters were irrigated weekly with secondary-treated domestic effluent (10 mm h–1 for 5 h; ~2650 mm yr–1), while the other half received rainfall only (~1260 mm yr–1). The soils were: a well-drained soil derived from sandy rhyolitic tephra (Pumice soil, Hewitt, 1998; Typic Udivitrand, Soil Survey Staff, 1996); a poorly drained soil formed in clayey estuarine alluvium (Gley soil, Hewitt, 1998; Typic Endoaquept, Soil Survey Staff, 1996); a well-drained soil with a high allophane content (Allophanic soil, Hewitt, 1998; Typic Hapludand, Soil Survey Staff, 1996); and a well-drained soil formed in a sand dune (Recent soil, Hewitt, 1998; Typic Udipsamment, Soil Survey Staff, 1996). Figure 2 shows soil profiles of the four soil types at the location of the lysimeters. Clay contents were a maximum of 5% in all horizons of the Pumice and Recent soils, ~30% in the Allophanic soil, and ~50 and 70% in the topsoil and subsoil, respectively, of the Gley soil. Saturated hydraulic conductivity (Ks) of each horizon was measured on 12 replicate A horizon cores (76 mm in height, 98 mm internal diameter) taken at the same locations where the lysimeters were excavated (Barton et al., 2005). Two and 4 yr after the beginning of the lysimeter study, four lysimeters of each soil and treatment combination were destructively sampled and Ks of the A horizons measured again on 16 cores (Barton, personal communication, 2003).


Figure 2
View larger version (103K):
[in this window]
[in a new window]
  		Fig. 2. Profiles of the four different soil types at the locations where the lysimeters were installed.

 
Volumetric moisture content was recorded in the top three horizons of one lysimeter per soil and treatment using CS615 WCRs (water content reflectometers) from Campbell Scientific (Campbell Scientific, Inc., 1996). The raw readings were corrected for temperature following Western et al. (2001) and subsequently transformed into volumetric moisture contents using horizon-specific calibration equations (Stenger et al., 2005). Daily averages of the WCR data recorded in 15-min resolution were used for the modeling. The leachate volumes were measured manually in typically 2-wk intervals. The data used here are averages of five lysimeters.

Model Setup
The HYDRUS-1D model (Simunek et al., 1998) solves the Richards equation (Richards, 1931) for one-dimensional vertical flow using the van Genuchten (1980)–Mualem (1976) soil hydraulic functions. The simulation was set up for a 2-yr period starting 1 Jan. 2001 to 1 Jan. 2003. Locally measured daily rainfall and estimated daily Penman–Monteith potential evapotranspiration (Allen et al., 1998) time series were used as atmospheric upper boundary conditions. If the lysimeter was irrigated, the applied irrigation volume was added to the daily rainfall value. The depth of the model was 0.7 m, which corresponds to the depth of the lysimeters. For each of the soil types, the model comprised three soil horizons (A, B, and C horizons). The depths of the horizons in the model corresponded with the measured data at the field sites where the lysimeters were installed. The bottom boundary condition of the model was set to be a seepage face, mimicking the way the bottom of the lysimeters were set up. This requires saturation to occur at the bottom boundary before water can leach out of the lysimeter. The calculation mesh consisted of 140 nodes. The average distance between nodes was 5 mm, the nodes being closer together near the top boundary condition. The initial conditions were set to a pressure head of –500 mm uniformly across the model domain.

The root water uptake was simulated according to Feddes et al. (1978), using the parameters for grass from the HYDRUS-1D database. Since no data were available on the rooting depth of the grass, rooting depth (d) was considered a parameter in the optimization algorithm. An exponential root distribution (R) varying between 0 at the bottom of the root zone and 1 at the soil surface was incorporated. Soil depth (z) ranged between 0 at the soil surface and d (de Rosnay and Polcher, 1998; Feddes et al., 2001):

Formula 1[1]
where c is a fitting constant depending on the biome considered and was assumed to be 4 after de Rosnay and Polcher (1998) for grassland.

Sensitivity Analysis
The classical statistical approach for assessing parameter uncertainty is based on a multinormal approximation of the probability density function of the model parameters around the estimated optimum; however, as shown by Duan et al. (1992) and more recently by Vrugt et al. (2002), this multinormal approximation around the optimum is often insufficient and misleading. Therefore, a sampling-based global sensitivity analysis (Spear and Hornberger, 1980; Mertens et al., 2005) was adapted in this study. The method is based on Monte Carlo sampling, where a number of parameter sets are randomly generated (uniform distribution) within their respective feasible regions and the corresponding model output is evaluated (Spear and Hornberger, 1980; Madsen, 2000a). In this study, 20 000 parameter sets were randomly generated and the corresponding model output was evaluated. The parameter sets were subsequently sorted with respect to the objective function value (e.g., goodness-of-fit measure) and two subsets were formed consisting of the best and worst parameters. The distributions of the two samples (X1 and X2) were compared using the two-sample Kolgomorov–Smirnov test. The test statistic (t) is given by

Formula 2[2]
where F(X1) and F(X2) are the sample cumulative distribution functions of X1 and X2, respectively. The larger the test statistic, the higher the probability that the two samples have different distributions. The test statistic can be associated with a particular probability level (p, significance level) for a proper assessment of the significance of the sensitivity of the parameter (Conover, 1980).

For the SA (sensitivity analysis), the Pumice soil's moisture content and leachate volumes were used as observations. The modeled soil profile consisted of three soil horizons: the A horizon between 0 and 23 cm, the B horizon between 23 and 40 cm, and the C horizon between 40 and 70 cm. For each of the three soil horizons, the sensitivity of the model to five soil hydraulic parameters was evaluated: Ks, {theta}s (saturated volumetric water content), {theta}r (residual volumetric water content), {alpha}, and n (both shape parameters). Additionally, d was also considered a parameter. In total, the sensitivity of 16 parameters was evaluated. In the SA as well as in the optimization procedure described below, the same upper and lower parameter bounds were used for all soil horizons. The parameter bounds were therefore selected to be very wide so that they comprised soil hydraulic parameter combinations that might be representative for all soil horizons found in the four different sol types. The minima and maxima are presented in Table 1.


View this table:
[in this window]
[in a new window]
  		Table 1. Lower and upper parameter bounds used in the sensitivity analysis and in the shuffled complex evolution parameter optimization algorithm.

 
The 20 000 parameter sets were randomly generated for each parameter between the bounds presented in Table 1. Since Ks, {theta}r, and {alpha} ranged over several magnitudes, a random number was generated between the natural logarithm of the bounds and thereafter back-calculated for use in the model. For a particular parameter set, the daily soil moisture contents at three observation nodes (7.5, 18, and 35 cm) were extracted from the simulation results as well as the leaching volume. The RMSE was calculated against the observed soil moisture contents at the three respective depths and the leachate volumes:

Formula 3[3]
where n is the number of observations (daily soil moisture content or biweekly leachate volume), Obsi is the observed moisture content or leachate volume, and Simi is the simulated moisture content or leachate volume. The average biweekly leachate volume of all five lysimeters was used in the objective function. Since only one of the five lysimeters of each soil type was equipped with WCR probes, the soil moisture content data from this lysimeter were used in the objective function. For the calculation of the RMSE, the first 3 mo of simulation were not included since this period was considered initialization.

Multiobjective Inverse Modeling for Soil Parameter Estimation
As shown in Mertens et al. (2005), the SCE algorithm (Duan et al., 1992) has proven to be a successful automatic calibration method for the estimation of effective parameter values of unsaturated zone parameters. For a detailed description of the SCE algorithm, see Duan et al. (1992). The most important algorithmic parameter is the number of complexes p. Sensitivity tests have shown that the number of parameters comprised in the optimization procedure is the primary factor determining the proper choice of p. In general, the larger the value of p, the higher the probability of converging into the global optimum but at the expense of a larger number of model simulations. Kuczera (1997) suggested setting p equal to the number of parameters included in the optimization; in this study, p was set to 8 (Mertens et al., 2005). Other SCE control parameters were set to values recommended by Duan et al. (1994), which corresponds to a population size of 216. At each location, the same upper and lower bounds as the ones used in the sensitivity analysis were used in the SCE algorithm (Table 1). On the basis of test optimizations and visual inspection of the convergence of objective function and "effective" parameters, the stopping criterion for the optimization routine was set to 5000 model evaluations. In all cases discussed below, the SCE algorithm converged (change in objective function <0.1% within three consecutive loops) within 5000 model evaluations.

The RMSE was chosen as the objective function for the optimization. In this study, two objectives (F1 and F2) were considered. Objective F1 evaluated the goodness of fit of the model against soil moisture contents and is defined as the sum of the RMSE values at the three different depths, whereas objective F2 evaluated the goodness of fit of the model against leachate volume. The result of a multiobjective calibration is not a single unique set of parameters but consists of the "Pareto set" of solutions (Deb, 2001; Gupta et al., 1998; Madsen, 2003; Mertens et al., 2004). In the case of two objectives, points on the Pareto front have the characteristic that no other points in the parameter space have smaller values for both F1 and F2. The Pareto front corresponds to the optimal solutions with a trade-off between different objectives, i.e., a point on the Pareto front with a better fit to the observed soil moisture content time series than any other point on the Pareto front will necessarily have a worse fit with respect to the measured leachate volume.

The multiobjective calibration problem was solved by aggregation of the two objectives into one measure. To compensate for differences in magnitude between both objective terms, the two objectives were transformed to a common distance scale (Madsen, 2003; Mertens et al., 2004). The transformation functions were deduced from the initial population in the SCE algorithm. Both objective functions were first normalized and then transformed to obtain the same distance from the origin to the optimum of the initial sample. The transformation function (g) is given by

Formula 4[4]
where Fi is the objective function of the ith objective (i = 1, soil moisture content; i = 2, leachate volume), {sigma}i is the standard deviation of the ith objective function of the initial sample, and {varepsilon}i is a transformation constant given by

Formula 5[5]
The aggregated objective function is given by

Formula 6[6]
where {omega} is a weighting factor. If {omega} = 0.5 both objectives are given the same weight and the SCE optimization provides a balanced optimum (Madsen, 2003; Mertens et al., 2004). In this study, the Pareto front was calculated by performing several individual SCE optimizations using different weights (0 ≤ {omega} ≤ 1). In the two-dimensional multiobjective framework as it was applied in this study, three "extreme" points on the Pareto front can easily be determined. The first extreme point is the result of the SCE optimization without incorporating leachate volume ({omega} = 1) in the objective function as expressed in Eq. [6]. The second extreme point corresponds with the balanced optimum and is the result of the SCE optimization with both objectives having equal weights ({omega} = 0.5). The third extreme point is the result of the SCE optimization incorporating only leachate volume ({omega} = 0). To limit computational time, only these three "extreme" weights (15000 simulations) were considered for the estimation of the Pareto front. This necessarily means that, in between, an approximation of the "true" Pareto curve was calculated and that possibly nonconvex parts of the Pareto front might not be found; however, to achieve our objectives it was not necessary to describe the Pareto front exactly. This study used the overall shape of the Pareto front for the evaluation of the trade-offs between different objectives and hence we believe this approximation suffices to address the issues discussed here.

Multiobjective Inverse Modeling as Model Verification Technique
A question of interest is whether the conceptual model of equal soil hydraulic parameters for both treatments is suitable for the modeling of leachate volume and soil moisture content for both treatments simultaneously. In other words, can soil moisture content and leachate volumes from irrigated lysimeters be modeled using the parameters estimated from unirrigated lysimeter data? The most straightforward approach to answer this question is probably to perform a model run using the soil hydraulic parameters corresponding to the balanced optimum of the unirrigated treatment and evaluate its goodness of fit against leachate volumes and soil moisture contents for the irrigated treatment. To evaluate this goodness of fit, however, a judgment call would have to be made based on some goodness-of-fit measure (e.g., RMSE). There is no objective way to evaluate whether the parameters estimated from the unirrigated treatment perform "well enough" in the irrigated treatment. An alternative approach to answer the question is to use both unirrigated and irrigated data in a single objective function and evaluate whether an acceptable fit for both treatments is obtained. Again, some kind of judgment call based on a performance measure would have to be made. In this study, a multiobjective framework was formulated to answer this question in a more robust and objective manner.

A new objective function was formulated that evaluates the goodness of fit of a parameter set with respect to the balanced objective function of both the irrigated and unirrigated lysimeters:

Formula 7[7]

Formula 8[8]

Formula 9[9]
Objective Funirr evaluates the goodness of fit (RMSE) of the model against the balanced objective function for the unirrigated lysimeter, whereas objective Firr evaluates the goodness of fit (RMSE) of the model against the balanced objective function for the irrigated lysimeter. The aggregated objective function now evaluates the goodness of fit of the model to the soil moisture content and leachate flux for both the irrigated and unirrigated lysimeters simultaneously. This means that during the optimization, the model was run twice for each parameter set: once using the boundary conditions of the unirrigated lysimeter and evaluated against the corresponding observed time series and subsequently an identical analysis for the irrigated lysimeter treatment.


   RESULTS AND DISCUSSION
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
SUMMARY AND CONCLUSIONS
 REFERENCES
 
Sensitivity Analysis
The Monte Carlo analysis was performed and model water balance errors checked for each parameter combination. If water balance errors ended up being >2%, the parameter combination was discarded (7000 out of the 20000). Subsequently, the Kolmogorov–Smirnov test was performed for the two different treatments, i.e., irrigated and unirrigated. For both treatments, the best (X1) and worst (X2) parameter sets were selected for the RMSE against (i) soil moisture contents at three depths (equal to the sum of 7.5-, 18-, and 35-cm depth RMSE), (ii) leachate volumes, and (iii) a balanced RMSE where {omega} = 0.5 in Eq. [6]. Earlier work (Mertens et al., 2005) showed that the sample size of X1 and X2 is not crucial and therefore a sample length of 100 parameters was chosen. Consequently, the parameters were ranked according to their probability level (P) in the Kolgomorov–Smirnov test from most sensitive to least sensitive parameters with respect to the three different objectives. The result of the multiobjective SA for both treatments is presented in Table 2.


View this table:
[in this window]
[in a new window]
  		Table 2. Parameters{dagger} ranked according to their sensitivity level (high to low sensitivity) or significance level (P) with respect to soil moisture, leachate volume, and the balance between both as a result of the multiobjective Monte Carlo sensitivity analysis for both the unirrigated and irrigated treatments.

 
The SA revealed that parameter sensitivity differed between objectives and between the irrigated and unirrigated treatment. Retention curve parameters (in particular saturated moisture content) were the most sensitive parameters when only optimizing against soil moisture contents while conductivity curve parameters (in particular saturated hydraulic conductivity) were the most sensitive parameters when optimizing solely on leachate volume. The balanced approach resulted in a combination of both parameter types being sensitive. This suggests that when optimizing both the retention and the hydraulic conductivity curve, it is desirable to have both objectives incorporated in the objective function. In case of the irrigated treatment, fewer retention curve parameters were found to be sensitive with respect to the soil moisture content than the unirrigated treatment. This is a reflection of differences in information content present in both treatments. Research into data requirements for inverse modeling has led to the understanding that the information content of the data is far more important than the amount of data used (Vrugt et al., 2001). As shown below, the soil moisture content as well as the leachate volume time series in the unirrigated lysimeters varied across a much wider range than in the irrigated lysimeters where the soil remained close to saturation all year and leachate remained constantly high. Hence, the unirrigated time series contained more information about the shape of the retention curve than the irrigated treatment. Since {theta}r for all three layers was the only parameter found insensitive to the balanced cases of both the irrigated and unirrigated treatments, it was not incorporated in the further inverse modeling and fixed to a value of 0.01.

Multiobjective Inverse Modeling for Soil Parameter Estimation
Unirrigated Lysimeters
For each soil type, the SCE algorithm was run three times using {omega} = 0, 1, and 0.5 in Eq. [6]. Figure 3 shows the results of the inverse modeling for the Pumice soil type. It shows that single-objective optimization against either leachate volume ({omega} = 0, Fig. 3 a) or solely soil moisture content ({omega} = 1, Fig. 3 b) resulted in an acceptable fit with respect to either considered objective, but performed badly with respect to the objective not incorporated. In the case of the balanced optimum ({omega} = 0.5, Fig. 3 c), the SCE algorithm succeeded in finding an effective parameter set that yielded acceptable fits to both objectives simultaneously. Figure 4 shows the Pareto front corresponding to the multiobjective inverse modeling results presented in Fig. 3. The three extreme Pareto points are indicated on Fig. 4 as triangles. The other Pareto points shown in Fig. 4 are points that have the property that no other points have both lower F1 and F2 when considering all the evaluated points in the three optimization runs. The three "extreme" points are important because they determine the shape of the Pareto front. The Pareto front in Fig. 4 indicates a significant trade-off between goodness of fit to both objectives, as expected from Fig. 3. It is the shape of the Pareto, however, that is important for the further evaluation of the trade-off. Figure 4 shows a "sharp" Pareto front. Thus, despite the significant trade-off, the "sharp" structure of the front reveals that including a second objective to a single-objective optimization ({omega} = 0 or 1) results in only a minor decrease in the goodness of fit to the original objective while providing an acceptable fit to the newly incorporated objective. If a rather linear Pareto front was observed, the "cost" associated with including a second objective in terms of goodness of fit to the original objective would be high. Very similar fits to the soil moisture content and leachate volume time series were obtained for the remaining three soil types (Table 3). Additionally, all Pareto fronts for the four soil types were found to be very sharp. Therefore, it can be concluded that the balanced optimum ({omega} = 0.5) provides a good balance between the goodness of fit of both soil moisture content and leachate volume for all soil types.


Figure 3
Figure 3
View larger version (70K):
[in this window]
[in a new window]
  		Fig. 3. Result of the inverse modeling for the unirrigated Pumice soil type (a) {omega} (weighting factor) = 0, (b) {omega} = 1, and (c) {omega} = 0.5. Simulated and observed soil moisture contents at three depths as well as simulated and observed biweekly leachate and cumulative leachate volumes are plotted.

 

Figure 4
View larger version (96K):
[in this window]
[in a new window]
  		Fig. 4. Pareto front for the unirrigated Pumice soil type indicating the trade-off between goodness of fit to the soil moisture content [g1(F1)] and leachate volume [g2(F2)]; 11 000 out of a total of 15 000 simulations shown.

 

View this table:
[in this window]
[in a new window]
  		Table 3. Optimized transformed soil moisture content [g1(F1)] and leachate volume [g2(F2)] RMSE values for the four different soil types and the three considered weight factors ({omega}).

 
Figure 5 presents the normalized range of the Pareto optimal parameter values for the unirrigated Pumice soil type presented in Fig. 4. The parameter values are normalized according to their parameter bounds used in the optimization (Table 1). Moves along the Pareto front, a large variation is observed for most of the parameters. Similar large variations in parameters along the Pareto front were found by Madsen (2003) in the case of distributed hydrological modeling with multiple objectives. It must be stressed that the variability of the Pareto optimal parameter sets is based solely on the concept of multiobjective equivalence of parameter sets and not a statistical measure for parameter uncertainty. As stated by Madsen (2000b), for a proper assessment of the confidence of a model simulation, other methods for the assessment of the parameter uncertainty (e.g., Markov chain Monte Carlo) should be performed complementary to this approach.


Figure 5
View larger version (44K):
[in this window]
[in a new window]
  		Fig. 5. Normalized range of the Pareto optimal parameter values for the unirrigated Pumice soil type shown in Fig. 2 (thetaS = saturated volumetric water content; alpha and n are shape parameters; Ks = saturated hydraulic conductivity; d = rooting depth; 1, 2, and 3 = A, B, and C horizons). Extreme optimal Pareto parameter sets are indicated by bold and dashed lines.

 
Figure 6 compares the cumulative distributions of the measured and Pareto optimal top-layer Ks values for the four soil types. Note that the SA showed that the balanced model output was most sensitive to the Ks value of the top layer. The Ks measurements within the same soil type and performed at the same time range by factors of 10 up to 1000 (e.g., the Gley soil in 1999) between the lowest and highest measurements; however, it is commonly acknowledged that Ks measurements can show this kind of variability. The measurements in 1999 were performed on undisturbed soil sampled at the field site where the lysimeters were installed, while for all other years, the soil cores were taken during the destructive sampling of four of the unirrigated lysimeters of each soil type. The figure additionally shows that for the Pumice and Allophanic soil types, the cumulative distribution of the Pareto optimal Ks does not overlap with any of its measured equivalents. Without performing any further statistical analysis, it can be concluded that Pareto-optimal parameters may or may not differ from directly measured parameters or from estimated parameters based on laboratory measurement. A discussion of the different reasons why Pareto optimal parameters and measured parameters may differ significantly is given in Mertens et al. (2005): (i) the difference between the scale on which Ks is measured (100-cm3 rings) and the model scale (one-dimensional column of 70-cm depth), (ii) the measurement and the consequent interpretation technique used, (iii) the global optimum of the objective function has been reached by the autocalibration algorithm, and (iv) the hydrological modeling has uncertainty (input uncertainty, parameter uncertainty, and conceptual model uncertainty). Although it is hard to quantify, we believe the dominant reason that caused the wide gap between measured and effective hydraulic conductivities in this study is related to the scale effect of the measurements. Mallants et al. (1997) showed that the measured hydraulic conductivity as well as its coefficient of variation decreases with increasing column size. A possible conceptual model error in this study could be the absence of preferential flow in the model; however, it is unlikely that in the unirrigated lysimeters' preferential flow played a dominant role since preferential or macropore flow only becomes significant under saturated conditions, which hardly occurred in the unirrigated lysimeters. Additionally, McLeod et al. (2001) concluded, on the basis of Br– breakthrough curves measured on the same soil type, that bypass flow is not dominant.


Figure 6
View larger version (36K):
[in this window]
[in a new window]
  		Fig. 6. Cumulative distributions of measured (1999, 2001, and 2003 or 1999, 2002, and 2004) and Pareto optimal top-layer saturated conductivities (Ks) for all unirrigated soil types.

 
Irrigated Lysimeters
As a comparison of results from the SCE applied to unirrigated lysimeters, the multiobjective inverse modeling approach was applied to the data for the irrigated Pumice soil type. As shown in Fig. 7 , the SCE algorithm again succeeded in finding parameter sets that fit both leachate volume and soil moisture content simultaneously. The Pareto front is not presented here but looks almost identical to the unirrigated Pareto front (Fig. 4). As in the unirrigated treatment, despite the significant trade-off, the "sharp" structure of the front indicates that the balanced optimum provides a good balance between the goodness of fit of both soil moisture content and leachate volume. Figure 8 shows the cumulative distributions of the soil hydraulic parameters along the Pareto front for both the unirrigated and irrigated Pumice soil type. The figure shows that for most soil hydraulic parameters, the range of parameter variation along the Pareto front is larger for the irrigated treatment than the unirrigated treatment. A possible reason for this is the fact that the irrigated time series contains less information about the shape of the soil retention curves, as is also suggested by the results of the SA. Without performing any statistical analysis, it is obvious that significant differences exist between parameter distributions along the Pareto fronts of both treatments. Most obvious are the higher {theta}s and Ks values found along the Pareto front of the irrigated treatment than the unirrigated treatment.


Figure 7
Figure 7
View larger version (74K):
[in this window]
[in a new window]
  		Fig. 7. Result of the inverse modeling for the irrigated Pumice soil type (a) {omega} (weighting factor) = 0, (b) {omega} = 1, and (c) {omega} = 0.5. Simulated and observed soil moisture contents at three depths as well as simulated and observed biweekly leachate and cumulative leachate volumes are plotted.

 

Figure 8
View larger version (38K):
[in this window]
[in a new window]
  		Fig. 8. Comparison between cumulative soil hydraulic parameter distributions along the Pareto front for the irrigated (Irr) and unirrigated (Non-Irr) Pumice soil type (thetaS = saturated volumetric water content; alpha and n are shape parameters; K = hydraulic conductivity; subscripts 1, 2, and 3 = A, B, and C horizons).

 
Multiobjective Inverse Modeling as Model Verification Technique
The trade-off between goodness of fit to the unirrigated and irrigated treatments for the Pumice soil type is presented in Fig. 9 . Similar to the Pareto fronts observed for leachate volume and soil moisture content in both the irrigated and unirrigated treatments, this Pareto front indicates a significant trade-off between both objectives. Unlike the Pareto fronts between goodness of fit to leachate volume and soil moisture content, however, this front is not sharp.


Figure 9
View larger version (88K):
[in this window]
[in a new window]
  		Fig. 9. Pareto front for the Pumice soil type indicating the trade-off between goodness of fit to the unirrigated (objective function Fnon-irr) and irrigated (objective function Firr) treatments.

 
The balanced optimum is found close to the {omega} = 0 extreme, in which only the irrigated lysimeter is considered. The Pareto front shows that the SCE algorithm does not succeed in finding a parameter set that yields acceptable fits to both the unirrigated and irrigated lysimeters simultaneously. Since the lysimeters data used for the two treatments were randomly collected at the same field site, we have no reason to believe that there were systematic differences between the lysimeters used for each treatment before the start of the effluent irrigation; however, the study shows that the conceptual model of equal soil hydraulic parameters for both treatments is not valid. Different soil hydraulic parameters are needed to model leachate volume and soil moisture contents of the irrigated and unirrigated treatments. Even though this conclusion is solely based on a visual inspection of the Pareto front, it is more robust and objective than a judgment call based on a single goodness-of-fit measure in case one of the two approaches described above were used. Again it must be acknowledged that the presented Pareto front is only an approximation of the true Pareto front. It can be expected that as the complexity of the front increases, the aggregated weighting technique used in this study will fail to fully capture the Pareto front; however, we believe that the quantified solutions as well as the presented approximation of the Pareto front do demonstrate the existing trade-offs and support our conclusions.

The only difference between both treatments is the application of secondary-treated effluent. Therefore, the reasons for the different effective parameter sets between irrigated and unirrigated lysimeters must be related to the long-term physical effect irrigation has on the soil structure and possible differences in water transport processes. The only available data to investigate the physical effect of irrigation on soil hydraulic parameters are the measured top-layer Ks measured at different times. Figure 10 shows the measured top-layer Ks data for the Pumice soil type with time and between treatments. The original Ks measurements were made in 1999. In 2001 and 2003, four irrigated and four unirrigated lysimeters were destructively sampled and Ks measured again. It is apparent from Fig. 10 that the Ks values of the irrigated treatment were higher in 2001 (median of 6827 mm d–1) and 2003 (median of 4154 mm d–1) than the initial measurements in 1999 (median of 1350 mm d–1). On the contrary, the Ks of the unirrigated treatment appears to have increased only slightly between 1999, 2001 (median of 1589 mm d–1), and 2003 (median of 2513 mm d–1). The range of measured Ks values in both treatments increased markedly with time. Figure 10 suggests that irrigation increases Ks; however, no statistically sound conclusions can be drawn since only 16 measurements were taken at each time. Effluent irrigation has been found in most studies to decrease the hydraulic conductivity of the soil (Balks et al., 1997; Magesan et al., 1999); however, increases have been reported as well (Mathan, 1994). Unfortunately, this analysis cannot be statistically elaborated since the only measured data available are a limited number of top-layer Ks values. No measured retention curve data with time are available to validate the change in soil hydraulic parameters due to irrigation in a more robust way. It is clear from Fig. 8, however, that higher Ks values are required to fit leachate and soil moisture content for the irrigated treatment than for the unirrigated treatment.


Figure 10
View larger version (28K):
[in this window]
[in a new window]
  		Fig. 10. Cumulative distributions of the measured top-layer saturated conductivities (Ks) for the unirrigated (non-irr) and irrigated (irr) Pumice soil type in 1999, 2001, and 2003.

 
Next to the increased Ks values, a difference in water transport processes between the unirrigated and irrigated lysimeters might exist. As argued above, it is unlikely that preferential flow played a significant role in the unirrigated lysimeters; however, preferential flow in the irrigated lysimeters cannot be excluded. Even though McLeod et al. (2001) concluded that bypass flow was not dominant, they observed an early arrival of the Br– peak in the pumiceous soil type, which they attributed to "physical exclusion". Additionally, the irrigation rate used in this study was double that of McLeod et al. (2001). Therefore, higher soil moisture contents and saturation occurring more frequently increased the chances of preferential flow. Since this process is not present in the conceptual model, the model might be "in error". A model in error means that effective parameter values become dependent on the boundary conditions, which was observed in this study (i.e., effective parameter values are a function of irrigation or no irrigation). The error in the conceptual model related to the absence of preferential flow would result in increased effective Ks values compared with the effective Ks values estimated on the unirrigated lysimeters, as is obvious from Fig. 8. Both a physical increase in Ks as well as macropore flow in the irrigated lysimeters might have resulted in the fact that effective soil hydraulic parameters differed between treatments.


   SUMMARY AND CONCLUSIONS
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
REFERENCES
 
The SA shows that when optimizing both the retention and the hydraulic conductivity curve, it is desirable to have both objectives incorporated in the objective function. Additionally, sensitivity ranking of parameters differed between the irrigated and unirrigated treatment but both show that model results are not sensitive to the value of {theta}r. The SCE algorithm succeeded in finding "effective" parameter sets that fit both data types well from effluent-irrigated and unirrigated treatments. The Pareto fronts were "sharp", meaning that parameter sets exist that fit both leachate volume and soil moisture content simultaneously. Comparing a limited number of measured saturated conductivities for the top layer with the saturated conductivities along the Pareto front indicates that the measured and effective parameter sets do not necessarily coincide. It has been shown that measured and "effective" parameters may differ for various reasons (Mertens et al., 2005; Oliver and Smettem, 2005). This work showed that Pareto optimal parameter sets differ between the unirrigated and the irrigated treatments. The question therefore arose whether the conceptual model of equal soil hydraulic parameters for both treatments was valid. To answer this, the multiobjective inverse modeling technique was set up using the leachate volume and soil moisture content of both the unirrigated and irrigated Pumice soil type. The corresponding Pareto front was not sharp and indicated that no effective parameter set could be found that could adequately describe leachate volumes and soil moisture contents for both treatments simultaneously. That result confirms the earlier finding from the single optimization analyses that had resulted in very different parameter sets for the two treatments. As the soil lysimeters used for both treatments were installed at the same site, the only difference between the treatments was the top boundary condition.

We believe that the reasons for the difference in the Pareto optimal parameter sets and the unsharp Pareto front must be related to the physical effect that effluent irrigation has on the soil hydraulic parameters and possible differences in water transport mechanisms. Comparing measured top-layer Ks with time suggests that Ks values for the Pumice soil type increased due to long-term effluent irrigation. Unfortunately, no measured soil water retention data with time are available to validate the change in soil hydraulic parameters due to irrigation in a more robust way. Although there is no physical evidence to support the existence of macropore flow in the irrigated lysimeters, it cannot be excluded. This would result in the model being "in error" and would therefore make the effective parameters a function of the boundary conditions.

This study confirms that the SCE optimization framework is capable of estimating soil hydraulic parameters in a multiobjective context and emphasizes the importance of the corresponding shape of the Pareto front. Pareto optimal parameters vary substantially along the Pareto front according to the weight given to each of the objectives. Analysis of the Pareto front allows identification of possible "errors" in the conceptual model and is therefore useful for conceptual model verification.


   ACKNOWLEDGMENTS
 
We would like to thank the New Zealand Foundation for Research Science and Technology (FRST) for funding this work as part of Lincoln Environmental Research's Groundwater Quality Protection Programme. Appreciation is also given to Landcare Research Ltd for providing the leachate volume data and access to the lysimeter facility. Thanks also to Malcolm McLeod from Landcare Research Ltd for supplying the soil profile photos. We are grateful to Diederik Jacques, Jan Diels, and Catherine Moore for their useful comments and critical review.


   REFERENCES
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES