HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (15) Citing Articles via Google Scholar Google Scholar Articles by Vrugt, J. A. Articles by Bouten, W. Search for Related Content PubMed Articles by Vrugt, J. A. Articles by Bouten, W. Agricola Articles by Vrugt, J. A. Articles by Bouten, W. Related Collections Soil Physics Statistics Soil Models

DIVISION S-1—SOIL PHYSICS

Validity of First-Order Approximations to Describe Parameter Uncertainty in Soil Hydrologic Models

Institute for Biodiversity and Ecosystem Dynamics, Section Physical Geography, Univ. of Amsterdam, The Netherlands, Nieuwe Achtergracht 166, Amsterdam, 1018 WV

* Corresponding author (j.vrugt{at}science.uva.nl)

	ABSTRACT

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

 
Model nonlinearity and parameter interdependence violate the use of a first-order approximation to obtain exact confidence intervals of parameters in soil hydrologic models. In this study, the posterior distribution of parameters in soil water retention and hydraulic conductivity functions is examined using observed water retention data and a laboratory transient multistep outflow experiment. Parameter uncertainties obtained with traditional first-order approximations and uniform grid sampling strategies were compared with those obtained using the Metropolis algorithm, a Markov Chain Monte Carlo (MCMC) sampler. A diagnostic measure, based on multiple sequences generated in parallel, was used to test whether convergence of the Metropolis sampler to the posterior distribution had been achieved. Most significantly, as the Metropolis algorithm can cope with rough response surfaces generated by the objective function used, it not only successfully infers the multivariate posterior probability distribution of the model parameters, but also provides valuable insights in parameter interdependence in the full parameter space.

Abbreviations: CV, coefficient of variation • MCMC, Markov Chain Monte Carlo • RMSE, root mean squared error • SCE, shuffled complex evolution • VG, van Genuchten

	INTRODUCTION

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

 
NUMERICAL SIMULATION MODELS of water flow and solute transport in the vadose zone are important tools in environmental research. The accuracy of predictions with these models heavily relies on accurate estimates of the unsaturated soil hydraulic parameters. Because direct assessment of the hydraulic parameters is generally impossible, estimation is usually based on fitting a numerical solution of the Richard's equation to observations collected during an experiment. In recent years, much effort has been directed to finding the most likely parameter set from experimental laboratory experiments, such as Multi-step outflow (Eching and Hopmans, 1993; Van Dam et al., 1994), evaporation experiments (imnek et al., 1998a; Romano and Santini, 1999) and from field data (Abbaspour et al., 1999; Musters and Bouten, 2000; Vrugt et al., 2001b, c; among many others). Surprisingly, relatively little attention has been given to a realistic assessment of parameter uncertainty under the different types of experimental conditions. We share the recent opinion by Durner et al. (1997) that improved interpretation of parameter uncertainty can yield valuable information to enable a better judgment of the limits of our theoretical understanding of unsaturated water flow in soils.

Parameter estimates obtained from calibrated soil hydrologic models are generally error-prone, because the data used for calibration contain measurement errors and because the model never perfectly represents the system under study. Confidence intervals on the calibrated model parameters can be used to express the degree of uncertainty in these quantities, but calculation of confidence intervals is far from straightforward because the solution of Richards' equation is generally a highly nonlinear function of the model parameters (Christensen and Cooley, 1999). Without a realistic assessment of parameter uncertainty, it is not wise to undertake with any confidence tasks such as evaluating confidence limits on future hydrological responses or assessing the relationships between model parameters and easily measurable basic soil properties.

A frequently employed approach, which is particularly popular in the area of vadose zone hydrology, is to obtain confidence intervals of parameters by utilizing a first-order approximation to the model function near its minimum (see Carrera and Neuman [1986] in groundwater hydrology; Kool and Parker [1988] in unsaturated soil water flow; Kuczera and Parent [1988] in rainfall-runoff modeling). As this classical first-order approximation does not account for correlations between the parameter estimates, computed standard errors can appear too favorable (Hollenbeck and Jensen, 1998; Christensen and Cooley, 1999). Indeed, case studies have shown that these linearly calculated confidence intervals are generally smaller than their corresponding nonlinear counterparts (Cooley, 1993; Christensen and Cooley, 1999). In general, the first-order approximation yields reasonable results provided that the estimated parameter uncertainty does not extend beyond the range for which a first-order approximation of the model equation applies. In soil hydrologic models, however, parameter interdependence and model nonlinearity violate the use of this first-order approximation to obtain exact confidence intervals of the parameter.

Alternatively, a robust but computationally intensive method for the calculation of confidence intervals, contrasting the classical first-order approach, is the generation of contour plots (Toorman et al., 1992; Gribb, 1996; Romano and Santini, 1999; Vrugt et al., 2001a). This kind of exhaustive uniform grid sampling requires discretizing the parameter space and computing the objective function for each grid point, which is a rather primitive method. Recently, Abbaspour et al. (1997)(1999) used a similar kind of sampling strategy, entitled the Sequentual Uncertainty Fitting algorithm (SUFI) for estimating subsurface flow and transport parameters. However, as simple this approach might be, it requires massive computing resources for highly dimensioned parameter spaces. For example, if we wish to sample, on average, a parameter within a hypercube with resolution equal to one tenth of the parameter range in a six-parameter space, we need 10⁶ model runs. Moreover, failure to maintain an adequate sampling density may result in undersampling or too course sampling in certain regions of the parameter space. Consequently, this may lead to errors in the computed parameter confidence intervals.

The objective of the present study is two-fold. First we present an effective and efficient MCMC sampler for assessing parameter uncertainty in soil hydrologic models. The Metropolis algorithm has received considerable attention in the last decade in the Bayesian statistics literature and is intimately related to a probabilistic optimization technique known as simulated annealing (Kirkpatrick et al., 1983). Additionally, we present a diagnostic measure to test whether convergence of the MCMC sampler to a stationary target distribution has been achieved. A second objective is to compare the size of the confidence intervals computed using the Metropolis sampler with those obtained using traditional first-order approximations and uniform grid sampling strategies. Two case studies of increasing complexity serve to illustrate the potentials of the presented algorithm. The first case considers estimation of soil water retention parameters and their posterior distribution in the commonly used soil water retention model of van Genuchten (1980). The second case involves estimation of the posterior distribution of the unsaturated soil hydraulic parameters using a laboratory multistep outflow experiment. In this study, we are especially concerned with the value of additional soil water pressure head measurements within the soil core for the identification of the soil hydraulic parameters.

	MATERIALS AND METHODS

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

 
Experiments often aim at assessing the relationship between a response or output variable , subject to error, and input variables. Consider the linear regression model,

[1]

where = (₁, ₂,..., _m) is a m x 1 vector of model outputs, X = (X₁₁,..., X_mp) is an m x p matrix of input values, ß = (ß₁, ß₂,..., ß_p) is a vector of p unknown parameters and e is a vector of statistically independent errors with zero expectation and constant variance ².

We assume that the mathematical structure of the model is known to be correct and fixed and that a uniform prior distribution between the realistic upper and lower bounds on each of the model parameters is specified a priori (thereby defining the feasible parameter space). Taking a Bayesian perspective, the objective of model calibration is to infer the posterior probability distribution, p(ß|D), which describes what is known about the model parameters ß given the Data D and specified prior information.

Traditional First-Order Approximation
Let's assume that little is known about the model parmeters ß a priori relative to what the experimental data will tell us, that is, p₀(ß|D) is uniform on the parameter space. For the linear regression model, the posterior distribution of ß is exact and then expressed as (Box and Tiao, 1973),

[2]

where, ß_opt is the vector of most probable parameters,

[3]

Applying the profile likelihood ratio concept (Kalbfleish and Sprott, 1970) yields the confidence intervals (1 - ) of the parameters,

[4]

Note that the profile likelihood ratio concept is similar to confidence intervals generated using the asymptotic ²-distribution approach (e.g., Hollenbeck and Jensen, 1998). All relevant inferences about ß can be made from the knowledge that the posterior distribution of ß is the multivariate normal distribution N_p . Hence, the marginal posterior distribution of ß_i is the multivariate normal distribution, N_i , where _ii is the ith diagonal element of .

This posterior probability region is exact for linear models, but only approximate for nonlinear models. Model nonlinearity and parameter interdependence undermine the use of this first-order approximation to obtain exact confidence intervals of model parameters in nonlinear models. The first-order approximation only yields reasonable results provided that the linearization of the objective function near its minimum is valid over the domain of ß for which there is significant uncertainty. However, such an approximation is often not valid. Besides strong parameter interdependence, the surface of p(ß|D) can deviate remarkably from the multinormal distribution. It may have multiple optima, discontinuous first derivatives and curving multidimensional ridges in the response surface, that is often the case when dealing with soil hydrological models (e.g., Rasiah, 1992; Duan et al., 1992; Musters and Bouten, 2000). Moreover, the ellipsoid region, defined by Eq. [2] may represent a very poor approximation of parameter uncertainty, especially in the case of a strong hyperbolic-banana-shaped curvature in the p(ß|D) surface. This is indicative for strong parameter interdependency.

Monte Carlo Sampling of Posterior Distribution
A MCMC method for assessing parameter confidence intervals in nonlinear models is based on the idea that rather than compute a probability density, p(ß|D), it is sufficient to have a large random sample drawn from p(ß|D) that approximates the form of the density. Diagnostic measures of central tendency and dispersion of the posterior distribution can be computed, using the mean and standard deviation of the large sample. Using the profile likelihood ratio concept, p(ß|D) can be approximated by finding those parameter sets for which Eq. [4] holds. This directly leads to the question of how to efficiently sample from p(ß|D). Two generic approaches to sampling from the density are considered, namely uniform grid and a Metropolis sampling scheme.

Uniform Grid Sampling
One sampling strategy that is especially popular in the field of vadose zone hydrology is the uniform sampling strategy. To not favor any initial distribution of the parameters, a uniform prior distribution is used. The algorithm proceeds according to the following steps:

• Step 1: Discretize the domain for each parameter ß_i, i = 1,2,..., p into N equidistant discrete points. Thus, N^p parameter sets are sampled.
  
• Step 2: Set n = 1, and evaluate the posterior density for parameter set ßⁿ according to (Box and Tiao, 1973),
  

  [5]

  
  

[6]

where y_j is the measurement and _j the corresponding model prediction. The transformation g(·) allows to handle the nonconstant error variance situations which are common in soil hydrology, where the variance of the measurement error associated with y often varies with magnitude, that is heteroscedastic (Sorooshian and Dracup, 1980). Notice, that when squared, M is simply the familiar ‘sum of squared residuals’ function that is commonly minimized in model calibration. Specific information about the applied transformations can be found in the Case Study section.

• Step 3: Increment n. If n is less than N^p, return to step 2.
  
• Step 4: Use Eq. [4] to compute which of the N^p sampled parameter sets are acceptable within a prespecified confidence level (1 - ). These samples are used to estimate statistical measures of the posterior distribution for each of the parameters, such as mean, variance, etc.
  

The method of uniform grid sampling is a rather primitive method, as it requires massive computing times for large p. This may be overcome by using a small N, but this then results in inaccurate measures of central tendency and dispersion for the computed posterior probability distribution.

Markov Chain Monte Carlo Samplers
Markov Chain schemes represent a general method for sampling from the posterior distribution p(ß|D). A Markov Chain is generated by sampling ß⁽ⁿ⁺¹⁾ z(ß|ß⁽ⁿ⁾). This z is called the transition kernel of the Markov Chain. Consequently, ß⁽ⁿ⁺¹⁾ depends only on ß⁽ⁿ⁾ and not on ß⁽⁰⁾, ß⁽¹⁾,..., ß^(n-1). The MCMC sampler will converge to its stationary distribution as n .

The Metropolis Algorithm. The most general and earliest MCMC algorithm, known as the Metropolis algorithm (Metropolis et al., 1953), is given as follows,

1. Randomly start at a location in the feasible parameter space, ß = ß and compute the posterior density of this point, p(ßⁿ|D) using Eq. [5] and [6].
   
2. Sample a new parameter set ß⁽ⁿ⁺¹⁾ using the multi-normal distribution as proposal distribution,
   

   [7]

   where ß⁽ⁿ⁾ is the mean, _ß is the covariance matrix of ß and c is an scaling factor, typically in the range of 1 to 3 to ensure that one can sample from regions of p(ß|D) which are not adequately approximated by the multinormal distribution in Eq. [7].
   

3. Calculate p(ß⁽ⁿ⁺¹⁾|D) and compute the ratio = p/p,
   
4. Randomly sample a label Z from a uniform distribution over the interval 0 to 1.
   
5. If Z then accept the new configuration. However, if Z > then remain at the current position, that is, ß = ß.
   
6. Increment n. If n is less than a prespecified number of draws N, then return to Step 2.
   

The Metropolis algorithm will always accept candidate points (jumps) into a region of higher posterior probability but will also explore regions with lower posterior probability with probability . This algorithm is a MCMC sampler generating a sequence of parameter sets, {ß⁽⁰⁾, ß⁽¹⁾,.., ß⁽ⁿ⁾} that converges to the stationary distribution, p(ß|D) for large n (Gelman et al., 1997). Thus, if the algorithm is run sufficient long, the samples generated can be used to estimate statistical measures of the posterior distribution, such as mean, variance, etc. To speed up the convergence rate of the Metropolis sampler to the posterior target distribution, the covariance matrix of the proposal distribution, _ß was periodically updated using a sample of the ß's generated in the Markov Chain (Kuczera and Parent, 1998).

Creating a Proposal Distribution
To increase computational efficiency, it is recommended to initialize the Metropolis sampler with a proposal distribution, Eq. [7], which is as close as possible to the posterior target distribution, p(ß|D). We approximated such a distribution in two steps. First, we located the most likely parameter set of the multivariate posterior target distribution, using the Shuffled Complex Evolution (SCE) global optimization algorithm (Duan et al., 1992, 1993) to ensure that our initial ß⁽⁰⁾ values for the Metropolis sampler come from the posterior distribution. Second, once one or multiple optima are found using the optimization algorithm, at each mode the first-order approximation, Eq. [2], is used to generate the proposal distribution.

In brief, the SCE algorithm takes an initial population of parameter sets, randomly spread out over the feasible parameter space and proceeds by successively applying the evolution process, until the termination criteria is met. The sampled points constitute a population. The population is portioned into several communities (complexes), each of which is permitted to evolve independently using the downhill simplex geometric shape as a search direction. After a certain number of generations, the communities are mixed and new communities are formed through a process entitled shuffling. This enhances the survivability by a sharing of information about the search space, gained independently by each community. Duan et al. (1992)(1993) has shown that this series of operators results in a effective, efficient, and robust global optimization algorithm.

The two steps presented here for creating a proposal distribution help avoid pitfalls in a wide variety of applications. Finding the modes is helpful in enabling the MCMC samplers to begin roughly centered around the most likely parameter set of the posterior target distribution.

Monitoring Convergence of the MCMC Sampler
An important issue in MCMC sampling is convergence of the sampler to the stationary posterior distribution. Theoretically, the sampler converges in the limit as n , but in any applied problem one must determine how many draws to make with the sampler. The study by Kuczera and Parent (1998) has not explored the issue of convergence of the Metropolis algorithm.

Many authors have addressed the problem of drawing inferences from MCMC samplers, including Ripley (1987), Gelfland and Smith (1990), Geweke (1992), and Raftery and Lewis (1992) in the recent statistical literature. Gelman and Rubin (1992a) demonstrated that it is generally impossible to monitor convergence of an MCMC sampler using a single sequence (one random walk). A strategy recommended by Gelman and Rubin (1992b) is therefore to generate several independent parallel sequences, with starting points, ß⁽⁰⁾'s sampled from the proposal distribution. Convergence of the MCMC sampler can then be monitored using between sequences as well as within sequence information. Hence, early on, after performing only a limited amount of draws, distributional estimates (mean and variance) of the various parameters between the parallel sequences will differ substantially, as the individual sequences have explored different parts of the parameter space, presumably close to their initialized starting value. After increasing the number of draws within each sequence, a larger part of the posterior target distribution is explored. Consequently, distributional estimates of the various parameters within and between individual sequences will tend to convergence to stationary values.

Mathematically we proceed in five steps (Gelman and Rubin, 1992b),

1. Independently simulate q 2 sequences with the Metropolis sampler outlined in The Metropolis Algarithm section, each of length 2n, with starting points drawn from the proposal distribution. To diminish the effect of the initial draws, discard the first n draws of each sequence, and focus the attention on the last n.
   
2. For each parameter of interest, calculate the mean value, using the qn drawn values for parameter ß_i.
   
3. For each parameter ß_i of interest, calculate the variance between the q sequence means, ß_i, each based on n samples ß_i
   

   [8]

   and calculate the average of the q within-sequence variances, s²_ßi, each based on n - 1 degrees of freedom,
   

   [9, 4.]

   Monitor convergence of the MCMC sampler by estimating the factor by which the scale of the current distribution of ß_i might be reduced if n . This potential scale reduction is estimated by
   

   [10]

   and declines to 1 as n . is the ratio of the current variance estimate, , to the within-sequence variance, W. Because of its minor contribution, the factor to account for the extra variance of the Student's t distribution is omitted from Eq. [10]. If the potential scale reduction is high, then we have reason to believe that proceeding with further draws may improve our inference about the posterior target distribution.
   

4. Once R is near 1 for each of the parameters ß_i, we can conclude that each of the q sequences of n draws is close to the target distribution, p(ß|D). Since a score of 1 is typically difficult to achieve, Gelman and Rubin (1992b) recommended using a value of 1.2 and less to declare convergence.
   

Inferences Concerning the Error Variance Term ()
A problem that is often overlooked in the use of Bayesian estimators is the assignment of the value of the error variance of the measurements, , in Eq. [2] and [5]. Assigning each residual a weight relative to the reliability of the corresponding measurement preserves all statistically relevant information. However, there is no unambiguously correct way in which to select an appropriate value of the error variance. We briefly discuss three different approaches for assessing a value for ²,

1. The value of the measurements error variance is fixed to the instrumental error of the measuring device. This approach was used by Hollenbeck and Jensen (1998) and applied to retention and transient unsaturated outflow observations. However, for the outflow observations, none of the model simulations was within the band of the measurement errors of the outflow observations, suggesting that even the optimal set of parameters does not provide an adequate fit to the data. In such as case, parameter estimates and any statements on their confidence are meaningless.
   
2. The value of the error variance of the measurements is inferred posteriori from the quality of the fit. Approach 1 and 2 are classical and an implicit assumption of zero model error.
   
3. The value of the measurements error variance is taken as a fixed percentage of the mean value of the observations used for parameter identification. This approach is especially popular in the field of rainfall-runoff modeling (see for example, Thiemann et al., 2001).
   

On the basis of these propositions we are faced with the possibility that there may not exist an objective statistically correct choice for the value of and therefore no statistically correct size and shape of the posterior distributions of the model parameters. For each of the case studies, the error variance of the measurements was inferred posteriori from the quality of the most optimal fit (Option 2) derived using the SCE global optimization algorithm, outlined in Creating a proposal distribution section.

Case Studies
The first-order approximation, uniform grid, and MCMC sampling strategy for estimating parameter uncertainty were tested for two different case studies of increasing complexity. The first case concerns a simple water retention model, whereas the second case study explores the utility of the Metropolis algorithm for identifying unsaturated soil hydraulic parameters, using a laboratory multistep outflow experiment.

Case Study I: Water Retention Model
One of the most commonly used models of capillary pressure saturation in soils is the water retention function of VG (van Genuchten, 1980),

[11]

in which (L³ L^-3) denotes water content, _s and _r (L³ L^-3) are the water content at full and residual saturation respectively, (L) is the soil water pressure head, (L^-1) is the inverse of the air-entry value and n is a unitless pore-size distribution index (van Genuchten, 1980).

The parameters _s, _r, , and n were calibrated on observed water retention measurements for a sandy and clayey soil by maximizing the posterior density defined in Eq. [5] using the SCE-University of Arizona (UA) algorithm (Duan et al., 1992; Vrugt et al., 2001d). Experimental data of a sandy (4520) and a clayey soil (2362) were taken from the UNSODA soil hydraulic properties database (Leij et al., 1996). The number between parentheses refers to the UNSODA code. For the uniform grid and Metropolis algorithm, prior ranges for each of the retention parameters were defined in Table 1.

[12]

where K_s (LT^-1) is the saturated hydraulic conductivity and l (-) is the exponent parameter in Mualem's equation (Mualem, 1976). As direct methods for determination of the () and K() curve require rather restrictive initial and boundary conditions and are therefore, relatively tedious. Numerical inversion is an attractive alternative as a means of determining both curves from a single experiment. In such an approach, the optimal soil hydraulic parameters are found by fitting a numerical solution of Richard's equation to observations of measured variables during the experiment. When a joint parametric description of retention and unsaturated hydraulic conductivity is assumed, inversion of the Richard's equation will yield parameter estimates that apply to both characteristics simultaneously.

The experimental work by Eching and Hopmans (1993) and Eching et al. (1994) confirmed the suitability of the multistep outflow method in combination with soil water pressure head measurements during drainage of the soil core for assessing the soil water retention and unsaturated soil hydraulic conductivity curve by numerical inversion. Observed outflow measurements combined with soil water pressure head measurement within the soil core were taken from the multistep outflow project (Example 7) stored in the HYDRUS-1D project manager (imnek et al., 1998b). The parameters of the unsaturated soil hydraulic functions in Eq. [11] and [12] were inversely estimated using the HYDRUS-1D model (imnek et al., 1998b), which numerically solves Richard's equation in one-dimension using a Galerkin-type linear finite element scheme. As this case study explores the usefulness of additional experimental pressure head data for parameter identification purposes, we examined the posterior target distribution for each of the parameters in Eq. [11] and [12] using two different data sets. First, observed outflow measurements were translated to average soil water contents within the soil core using the initial soil water content, to increase the identifiability of the saturated and residual water content (Vrugt et al., 2001a). Subsequently, this data set was used separately and in a joint combination with the soil water pressure head measurements to infer the posterior target distribution of the parameters. We assumed that the soil water pressure head measurements errors have a heteroscedastic (nonconstant) variance that is linearly related to its magnitude, which can be stabilized by the Box-Cox transformation (Box and Cox, 1974),

[13]

where is a transformation parameter, assumed to be 0.8 (see Sorooshian and Dracup, 1980). For each measurement set, the measurement errors were set identical to the root mean square error (RMSE) of the most optimal fit derived using the SCE global optimization algorithm (Duan et al., 1992; Vrugt et al., 2001d). This approach seemed most appropriate as adapting the relatively small measurement errors reported in literature for either soil water pressure head (Romano and Santini, 1999) or water content measurements (Vrugt et al., 2001a), would lead to rejection of the applied HYDRUS-1D model, making even the optimal parameter values virtually meaningless (see Inferences concerning the error variance term []). To avoid that the Metropolis algorithm, progressively samples outside realistic physical bounds of the model parameters, prior parameter ranges for each of the soil hydraulic parameters were defined in Table 2.

	RESULTS AND DISCUSSION

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

View larger version (17K): [in this window] [in a new window]   Fig. 1. Optimized soil water retention curves for the sandy and clayey soil, using the van Genuchten model. Measured retention observations are indicated with symbols.

View larger version (25K): [in this window] [in a new window]   Fig. 2. (A) Relationship between the jumprate used for sampling from the proposal distribution for the sandy soil and the amount of model calls needed to achieve convergence of the Metropolis sampler, (B) Evolution of the Gelman-Rubin convergence diagnostic for each of the van Genuchten parameters of the sandy soil.

This is also evident in Fig. 3a and 3b , which presents the sampled _r–n parameter space for the sandy soil for each of the four sequences at two different stages during the evolution. Fig. 3a corresponds to the sampled parameter space after 4000 samples being generated (Location A in Fig. 2B), whereas Fig. 3B corresponds to the sampled parameter space after convergence has been achieved (Location B in Fig. 2B). We discarded the first n = 500 draws from each sequence (burn-in), as it is unlikely that these initial draws come from the stationary distribution needed to construct correct posterior estimates. The different sequences are coded with different symbols. The two-dimension scatter plots of the sampled parameter space demonstrate that at early stages during the evolution, the individual paths tend to occupy different regions of the posterior surface. This low mixing of the sequences is associated with a relatively high value of the scale reduction factor, indicating no convergence. Nevertheless, after a sufficient amount of draws (Stage II), each of the individual sequences have been able to fully explore the posterior target distribution, resulting in a scale reduction factor smaller than 1.2, indicating convergence to a stationary distribution.

View larger version (39K): [in this window] [in a new window]   Fig. 3. Scatter plot of r–n samples generated by the Metropolis algorithm at two different stages during the evolution. Each sequences generated with the Metropolis sampler is coded with a different symbol. For more explanation see text.

View larger version (17K): [in this window] [in a new window]   Fig. 4. Scatter plot of 2000 combinations of -n (A), and r–n (B) parameters sampled for the clayey soil using the Metropolis algorithm.

As compared with traditional statistical inference, the Metropolis algorithm not only generates a substantial better description of parameter uncertainty in soil hydrological models, but also provides powerful information about parameter interdependency in the full parameter space. Although the results presented here are restricted to two different sets of retention observations (sandy and clayey soil) similar conclusions were drawn using other sets of water retention observations, taken from the UNSODA database.

Case Study II: Laboratory Multistep Outflow Experiment
Table 5 presents the posterior mean, standard deviation and CV values for the unsaturated soil hydraulic parameters, derived from utilizing a classical first-order approximation and from Metropolis sampling using the observed soil water content within the soil core. Based on the previous results each individual sequence was generated with jumprate equal to 0.5. With this jumprate, the Metropolis-sampler arrived at a stationary posterior target distribution of the six soil hydraulic parameters after approximately 15000 model evaluations. Again, this confirms the high sampling efficiency of the Metropolis sampler as compared with hypercube sampling strategies. A comparison between the size of the linear and nonlinear confidence intervals suggests that the Metropolis-derived standard deviations are generally smaller than their corresponding first-order counterparts. This contradicts the idea that the first-order covariance matrix, used to construct the linear confidence intervals, always provides a lower bound estimate of the parameter confidence intervals. To understand why, consider Fig. 5 that presents the corresponding univariate marginal posterior probability distributions for each of the unsaturated soil hydraulic parameters in the Mualem–VG model. Also indicated in Fig. 5 are the final optimized values of the unsaturated soil hydraulic parameters using the SCE global optimization algorithm of Duan et al. (1992)(1993). Most of the histograms of the parameters not only deviate remarkably from the normal distribution, which is an explicit assumption when utilizing the first-order approximation, but also exhibit multimodality. The presence of multimodality in the posterior distribution of the parameters points at nonuniqueness problems of the parameters, when using observed soil water content data only. Nevertheless, the SCE algorithm has been able to successfully identify the global optimum within this high-density region. Hence, the modes of the univariate posterior target distributions of the various parameters, obtained using the Metropolis sampler, coincide with the most optimal parameter values, identified using the SCE algorithm.

View larger version (33K): [in this window] [in a new window]   Fig. 5. Univariate posterior probability distributions of the soil hydraulic parameters s, r, , n, Ks and l using observed soil water content dynamics during the outflow experiment. The arrows indicate the most likely parameter values, derived using the Shuffled Complex Evolution global optimization algorithm.

As the Bayes theorem provides a mathematical formulation of how previous knowledge maybe combined with new knowledge it allows to continually update information about a set of parameters ß as different sets of observations are taken. Therefore, the Bayes theorem is suited to study the effect of pooling observed soil water contents and soil water pressure head observations in one likelihood function on the identifiability of the soil hydraulic parameters. In this case, the conditional likelihood is defined as the multiplicative constraint of the probabilities for describing the observed soil water contents and soil water pressure head observations simultaneously, according to Eq. [5]. Figure 6 shows the observed and simulated soil water content dynamics (I) as well as the measured and simulated soil water pressure heads (II) during the MSO experiment, using either soil water contents (A), or a joint identification using soil water pressure head data and soil water contents simultaneously (B). When using Dataset A only for identification purposes, the HYDRUS-1D model generates an excellent fit to the observed soil water contents, at the expense of overestimating the measured soil water pressure heads in the lower water content range. Adding measured soil water pressure heads as an extra objective in the likelihood function, increases the quality of the fit to observed soil water pressure heads, but simultaneously worsens the quality of the fit to observed soil water contents. Because of errors in the model structure (and other possible sources), it is usually not possible to find a single unique solution that simultaneously minimizes all objectives. Instead, it is common to have a set of solutions with the property that moving from one solution to another results in the improvement of one objective while causing a deterioration in one or more others. This multi-objective equivalence of parameter sets is more commonly referred to as pareto optimal in literature.

View larger version (36K): [in this window] [in a new window]   Fig. 6. Observed and simulated soil water content dynamics (I) as well as the measured and simulated soil water pressure heads (II) during the multistep outflow (MSO) experiment, using either observed soil water contents (A), or a joint identification using soil water pressure head data and soil water contents simultaneously (B).

Finally, we would like to emphasize that one should be particularly careful to draw conclusions about parameter inequivalance without recourse to examination of the uncertainty associated with the parameters (Durner et al., 1997; Wildenschild et al., 2001). For instance, although the most optimal values for Scenario A and B outlined in Table 6 differ substantially the 95% confidence intervals of suggests that this parameter might be equivalent in describing both soil water content and soil water pressure head dynamics in the soil core simultaneously.

	SUMMARY AND CONCLUSIONS

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

 
In this paper, we have considered the Metropolis algorithm as an effective and efficient MCMC sampler for assessment of parameter uncertainty in nonlinear models. Parameter uncertainties obtained with traditional first-order approximations and uniform grid sampling strategies were compared with those obtained using the Metropolis algorithm, an effective and efficient Markov Chain Monte Carlo sampler. A diagnostic measure, based on multiple sequences generated in parallel, was usedto test whether convergence of the Metropolis sampler to the posterior distribution had been achieved.

Two case studies of increasing complexity served to illustrate the potentials of the Metropolis algorithm. The first case study considered a simple water retention and demonstrated that the Metropolis algorithm not only generates a substantial better description of parameter uncertainty in soil hydrological models, as compared with traditional statistical inference, but also provides powerful information about parameter interdependency in the full parameter space. We illustrated how the structure that is induced in the joint distributions of the Metropolis samples, can be used to assess the suitability of experimental designs for identification of the soil hydraulic parameters.

The second case study considered a more complex transient outflow experiment and served to illustrate the usefulness of the Metropolis algorithm for identification of the soil hydraulic parameters and value of additional pressure head data for parameter identification purposes. Samples generated with the Metropolis algorithm using observed soil water contents during the outflow experiment only, demonstrated multimodality in the univariate distribution of the soil hydraulic parameters, implying nonuniqueness. Augmenting the observed soil water contents with additional soil water pressure head data, did not reduce parameter uncertainty, but increased the identifiability of the global minimum in the parameter space. The practical advantage of using the Metropolis algorithm for identifying soil hydraulic parameters from outflow experiments is that it successfully finds the global minimum in the parameter space, based on one measurement type only (cumulative outflow or soil water pressure head observations) as it can cope with rough response surfaces generated by the objective function being used. Moreover, the algorithm infers the posterior probability distribution p(ß|D), which describes what is known about the model parameters given the Data D. Most significantly, all case studies illustrated the limitations of traditional inference regarding parameter uncertainty based on first-order approximations.

	ACKNOWLEDGMENTS

 
The ideas for this research arose out of discussions within the International Study Group on Inverse modeling (ISGIM) at the Meeting in Mobile, Alabama, Fall 2000. For this, we especially would like to acknowledge Jacob H. Dane. The Earth Life Sciences and Research Council (ALW) partly supported the investigations of the first author with financial aid from the Netherlands Organization for Scientific Research (NWO).

Received for publication November 19, 2001.

	REFERENCES

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

 

• Abbaspour, K.C., M.Th. van Genuchten, R. Schulin, and E. Schläppi. 1997. A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters. Water Resour. Res. 33:1879–1892.
• Abbaspour, K.C., M. Sonnleitner, and R. Schulin. 1999 Uncertainty in estimation of soil hydraulic parameters by inverse modelling: example lysimeter experiments. Soil Sci. Soc. Am. J. 63:501–509.[Abstract/Free Full Text]
• Box, G.E.P., and G.C. Tiao. 1973. Bayesian inference in statistical analyses. Addison-Weslay Publishers Co., Reading, MA.
• Box, G.E.P., and O.R. Cox. 1974. The analysis of transformations. J. Roy. Statist. Soc., Ser. B, 26:211–252.
• Carrera, J., and S.P. Neuman. 1986. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information. Water Resour. Res. 22:199–210.
• Christensen, S., and R.L. Cooley. 1999. Evaluation of confidence intervals for a steady-state leaky aquifer model. Adv. Water Resour. 22:807–817.
• Cooley, R.L. 1993. Exact Scheffé-type confidence intervals for output from groundwater flow models, 2. Combined use of hydrogeologic information and calibration data. Water Resour. Res. 29:35–50.
• Duan, Q., S. Sorooshian, and V.K. Gupta. 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28:1015–1031.
• Duan, Q., S. Sorooshian, and V.K. Gupta. 1993. A shuffled complex evolution approach for effective and efficient global minimization. J. Optimization Theory Application 76:501–521.
• Durner, W., B. Schultze, and T. Zurmühl. 1999. State of the art of inverse modeling of inflow/outflow experiments. p. 661–681. In M.Th. Van Genuchten et al. (ed.) Proc. of the Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, 22–24 Oct. 1997. University of California, Riverside, CA.
• Eching, S.O., and J.W. Hopmans. 1993. Optimization of hydraulic functions from transient outflow and soil water pressure data. Soil Sci. Soc. Am. J. 57:1167–1175.[Abstract/Free Full Text]
• Eching, S.O., J.W. Hopmans, and O. Wendroth. 1994. Unsaturated hydraulic conductivity from transient multi-step outflow and soil water pressure data. Soil Sci. Soc. Am. J. 58:687–695.[Abstract/Free Full Text]
• Gan, T.Y., and S.J. Burges. 1990a. An assessment of a conceptual rainfall-runoff model's ability to represent the dynamics of small hypothetical catchments, 1, Models, model properties, and experimental designs. Water Resour. Res. 26:1595–1604.
• Gan, T.Y., and S.J. Burges. 1990b. An assessment of a conceptual rainfall-runoff model's ability to represent the dynamics of small hypothetical catchments, 2, Hydrologic responses for normal and extreme rainfall. Water Resour. Res. 26:1605–1619.
• Gelfland, A.E., and A.F.M. Smith. 1990. Sampling-based approaches to calculating marginal densities. J. Am. Stat. Assoc. 85:398–409.[ISI]
• Gelman, A., J.B. Carlin, H.S. Stren, and D.B. Rubin. 1997. Bayesian data analysis. Chapmann and Hall, London.
• Gelman, A., and D.B. Rubin. 1992a. A single series from the Gibbs sampler provides a false sense of security. p. 625–632. In J.M. Bernardo et al. (ed.) Bayesian statistics 4. Oxford Univ. Press, New York.
• Gelman, A., and D.B. Rubin. 1992b. Inference from iterative simulation using multiple sequences. Statistical Sci. 7:457–511.
• Geweke, J. 1992. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. p. 169–193. In J.M. Bernardo et al. (ed.) Bayesian Statistics 4. Oxford Univ. Press.
• Gribb, M.M. 1996. Parameter estimation for determining hydraulic properties of a fine sand from transient flow measurements. Water Resour. Res. 32:1965–1974.
• Hollenbeck, K.J., and K.H. Jensen. 1998. Maximum-likelihood estimation of unsaturated hydraulic parameters. J. Hydrol. 210:192–205.
• Ibbitt, R.P. 1970. Systematic parameter fitting for conceptual models of catchment hydrology, Ph.D. dissertation, Imp. Coll. of Sci. and Techn., University of London, London.
• Johnston, P.R., and D. Pilgrim. 1976. Parameter optimization for watershed models. Water Resour. Res. 12:477–486.
• Kalbfleisch, J.D., and D.A. Sprott. 1970. Application of likelihood methods to models involving large number of parameters. J. Royal Stat. Soc., Ser. B, 32:175–208.
• Kirkpatrick, S., C.D. Gellat, and M.P. Vecchi. 1983. Optimization by simulated annealing. Science 220:671–680.[Abstract/Free Full Text]
• Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem for transient unsaturated soil water flow. Water Resour. Res. 24:817–830.
• Kuczera, G., and E. Parent. 1998. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. J. Hydrol. 211:69–85.
• Kuczera, G., and M. Mroczkowski. 1998. Assessment of hydrological parameter uncertainty and the worth of multiresponse data. Water Resour. Res. 34:1481–1489.
• Leij, F.J, W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. Unsaturated soil hydraulic database, UNSODA 1.0 user's manual. Report EPA/600/R-96/095. U.S. EPA, Ada, OK.
• Linssen, H.N. 1975. Nonlinearity measures: A case study. Statistica Neerlandica. 29:93–99.
• Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. 1953. Equations of state calculations by fast computing machines. J. Chem. Phys. 21:1087–1091.
• Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
• Musters, P.A.D., and W. Bouten. 2000. A method for identifying optimum strategies of measuring soil water contents for calibrating a root water uptake model. J. Hydrol. 227:273–286.
• Raftery, A.E., and S.M. Lewis. 1992. How many iterations in the Gibbs sampler? p. 763–773. In J.M. Bernardo et al. (ed.) Bayesian statistics, Vol. 4. Oxford University Press, Oxford, UK.
• Rasiah, V., G.C. Carlson, and R.A. Kohl. 1992. Assessment of function and parameter estimation methods in root water uptake simulations. Soil Sci. Soc. Am. J. 56:1267–1271.[Abstract/Free Full Text]
• Ripley, B.D. 1987. Stochastic simulation. Wiley, New York.
• Romano, N., and A. Santini. 1999. Determining soil hydraulic functions from evaporation experiments by a parameter estimation approach: Experimental verifications and numerical studies. Water Resour. Res. 35:3343–3359.
• imnek, J., O. Wendroth, and M.Th. van Genuchten. 1998a. A parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.[Abstract/Free Full Text]
• imnek, J., M. ejna, and M.Th. van Genuchten. 1998b. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, IGWMC-TPS-70, International Ground Water Modeling Center, Colorado School of Mines, Golden, CO.
• Sorooshian, S., and J.A. Dracup. 1980. Stochastic parameter estimation procedures for hydrological rainfall-runoff models: Correlated and heteroscedastic error cases. Water Resour. Res. 16:430–442.
• Sorooshian, S., and V.K. Gupta. 1983. Automatic calibration of conceptual rainfall-runoff models: The question of parameter observability and uniqueness. Water Resour. Res. 16:430–422.
• Thiemann, M., M. Trosset, H. Gupta, and S. Sorooshian. 2001. Bayesian recursive parameter estimation for hydrological models. Water Resour. Res. 37:2521–2536.
• Toorman, A.F., P.J. Wierenga, and R.G. Hills. 1992. Parameter estimation of hydraulic properties from One-step outflow data. Water Resour. Res. 28:3021–3028.
• van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1994. Inverse method to determine soil hydraulic functions from multi-step outflow experiments. Soil Sci. Soc. Am. J. 58:647–652.[Abstract/Free Full Text]
• van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
• Vrugt, J.A., W. Bouten, and A.H. Weerts. 2001a. Information content of data for identifying soil hydraulic parameters from outflow experiments. Soil Sci. Soc. Am. J. 65:19–27.[Abstract/Free Full Text]
• Vrugt, J.A., M.T. van Wijk, J.W. Hopmans, and J.