| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Institute for Biodiversity and Ecosystem dynamics, University of Amsterdam, Nieuwe Prinsengracht 130, Amsterdam, 1018 VZ, the Netherlands
Corresponding author (j.vrugt{at}frw.uva.nl)
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
Abbreviations: MSO, multi-step outflow • MVG, Mualem–van Genuchten • OSO, one-step outflow • PIMLI, Parameter Identification Method based on the Localization of Information • SWIF, soil water in forested ecosystems model
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
The simplest form of parameter estimation is curve fitting in which measured data are represented by a static function with parameters that provide the best possible fit to the data (van Genuchten et al., 1991). A more complex form of parameter estimation is inverse modeling, where parameters are optimized while minimizing a suitable objective function that expresses the discrepancy between the output of a dynamic model and the measurements (Janssen and Heuberger, 1995). However, using this conventional approach, problems are often encountered with the non-uniqueness of the optimized parameters. Non-uniqueness leads to more than one set of parameters, each yielding minimum values for the objective function determined by local minima or by the same global minimum at more than one point in the parameter space (Gupta and Sorooshian, 1985; Duan et al., 1992).
Seemingly there was a widespread conviction that the best way to solve the non-uniqueness problem was to include additional and more accurate measurements (Eching and Hopmans, 1993; Van Dam et al., 1994). However, research into data requirements has led to the understanding that the information content of the data is far more important than the amount of data used for model calibration (Kuczera, 1982; Sorooshian et al., 1983; Gupta and Sorooshian, 1985; Yapo et al., 1996; Gupta et al., 1998). Therefore, we developed a parameter identification technique entitled PIMLI (Parameter Identification Method based on the Localization of Information) that treats the data in such a way as to preserve the specific information with respect to the various parameters and that uses the information content of data to identify the model parameters. For illustrating PIMLI we used the MSO experiment that only measures cumulative outflow, because of its well-known problems regarding the non-uniqueness of the identified parameters for unsaturated flow (Hopmans and 
im
nek, 1999).
The use of inverse methods for determining the unsaturated flow parameters from transient experiments was first reported by Zachmann et al. (1981). Kool et al. (1985a)(b) used this inverse modeling approach to determine soil hydraulic properties from one-step outflow (OSO) experiments, but experienced problems with the non-uniqueness of the parameter estimation. Further investigations of the inverse method demonstrated the need for additional 
(h) data (Hudson et al., 1991; Van Dam et al., 1992; Bohne et al., 1993) or tensiometer measurements inside the soil sample (Kool and Parker, 1988; Toorman et al., 1992; Eching and Hopmans, 1993) to overcome the problem of non-uniqueness of parameters. The benefit of including tensiometer measurements is apparent, as the optimized soil water retention is forced to match the observed (h) data. Another way to realize more reliable parameter estimates is to incrementally increase the pneumatic pressure in several steps, the MSO experiment, instead of a single pressure increment (Van Dam et al., 1994).
The objective of this study was to analyze the temporal variability of model sensitivity for the various soil hydraulic parameters for a MSO experiment that only measures cumulative outflow. These results were used to assess transient flow experiment potentials and limitations for parameter identification. For both purposes we used the PIMLI algorithm.
First, we evaluated the uniqueness of the inverse problem for a sandy soil using two-dimensional response surfaces of the objective function. These response surfaces were obtained by perturbing two selected soil hydraulic parameters around their true values, while maintaining the additional parameters constant at their true values (Toorman et al., 1992; imnek et al., 1998). Then, we illustrated the PIMLI for the same sandy soil, a silty soil, and a clayey soil. Our study was carried out with numerically generated error-free data. These artificial measurements are preferred because the soil hydraulic parameters are then known beforehand (Toorman et al., 1992; imnek et al., 1998).
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
 | (1) |
is the volumetric water content (L3 L-3), t is time (T), x is the spatial coordinate (L) (positive upward), and K(h) is the unsaturated hydraulic conductivity (L T-1). Because the user can control the upper and lower boundary in SWIF, it is possible to apply a zero flux upper boundary and suction at the lower boundary (Dirichlet condition), as required in the outflow experiment.
The soil hydraulic functions in SWIF are described by the Mualem–van Genuchten (MVG) model (Mualem, 1976; van Genuchten, 1980):
 | (2) |
 | (3) |
 | (4) |
(L3 L-3) denotes water content, s is the saturated water content (L3 L-3), r is the residual water content (L3 L-3), Ks is the saturated conductivity (L T-1), and 
, and 
are curve shape parameters.
Artificial Measurements
The artificial outflow measurements were calculated for a soil core with a diameter of 5.04 cm and a height of 5 cm, corresponding with a soil volume of 100 cm3. At the bottom of the soil core, we simulated a 0.7-cm-thick porous ceramic plate with a saturated conductivity of 0.0034 m d-1. Because the ceramic plate was considered part of the porous system in the numerical simulations, parameter values were chosen such that the ceramic remained fully saturated for the range of pressure head steps. As the initial condition, hydraulic equilibrium was assumed with a pressure head 
at the bottom of the soil core, and 
at the top. The following pressure head steps and time periods (in parentheses) were applied in the experiment: 
, and 
.
Response Surface Analysis
The commonly used objective function OF(b), which is normally minimized with the help of a classical parameter optimization algorithm, is defined as
 | (5) |
s, r, 
, n, Ks, and ), j represents the different sets of measurements (cumulative outflow, water content, and flux density), nj is the number of measurements within a particular set, q*j(ti) are measurements of type j at time ti, qi(ti,b) are the corresponding model predictions using the parameters in b, and wj and wi,j are weighting factors associated with measurement type and individual measurements, respectively. Because all the measurements were generated with the numerical model, we assumed that the errors associated with individual measurements of all types were identical. Therefore, wi,j was set equal to 1 for all the measurements. Differences in weighting between measurement types, as caused by differences in magnitude and their number nj, were normalized by dividing each data point by both the variance of the measurements of type j and the nj (Clausnitzer and Hopmans, 1995)
 | (6) |
j and nj denote the standard deviation and the number of j-type measurements, respectively.
We investigated the posedness of the conventional inverse solution using two-dimensional response surfaces of the objective function. We approached this question in a way similar to that done previously by Toorman et al. (1992) and imnek et al. (1998). The response surfaces are obtained by solving the objective function in Eq. [5] for many possible combinations of the selected pair of parameter values, while keeping the additional four parameters constant at their true values. These response surfaces reveal the occurrence of local minima, the presence of a well-defined global minimum, the parameter sensitivity, and parameter correlations. If response surfaces do not display a well-defined global minimum in the two-dimensional parameter planes, the conventional inverse parameter estimation technique may certainly be expected to be unsuccessful in a multidimensional plane. The response surfaces were calculated on a rectangular grid with parameter values corresponding with the sandy soil in Table 1. Each parameter domain was discretized into 40 equidistant discrete points with domain ranges as presented in Table 1, resulting in 1600 grid points for each response surface.
|
imnek and van Genuchten, 1996). Nevertheless, response surfaces are a suitable tool to obtain insight into the behavior of the objective function in the full parameter space.
PIMLI
PIMLI is a parameter identification method that for each parameter uses a separate subset of measurements with the highest information content for that particular parameter. If we want to identify those different subsets, we need a stochastic approach (Hollenbeck and Jensen, 1998) since analytical solutions for transient water flow are not available and cannot be derived. Therefore, the numerical model SWIF was first used to generate 1000 Monte-Carlo simulations of a MSO experiment with given initial and boundary conditions by randomly selecting combinations of the MVG parameters s, r, , n, Ks, and . The parameter sets were selected with a Latin-Hypercube sampling method (McKay et al., 1979) in which parameter ranges in Table 1 were used. We assumed no correlation between the different soil hydraulic parameters. The outcome of the Monte Carlo analysis resulted in a population of generated model outputs at every discrete time value. Thereafter, a reference run of measurements was simulated with SWIF using hydraulic properties corresponding with the sandy soil (Table 1). The cumulative outflow, its first derivative (flux density), and the average water content of this reference run as deduced from the cumulative outflow and the water content at the end of the experiment are presented in Fig. 1A
.
|
 | (7) |
 |
 |
Q (cm3), (cm3 cm-3), and F (m d-1) denote the error of cumulative outflow, water content, and flux density, respectively; ,end is the error in water content of the soil sample as caused by weighing at the end of the experiment, r denotes the radius of the soil core, and t is the time interval between two subsequent measurements, which was chosen as 0.05 d. The error of the outflow was determined from the resolution and accuracy of pressure transducers that are used for automated monitoring of the outflow dynamics. We assumed ,end to be 0.01 (cm3 cm-3).
A parameter set was accepted if its corresponding simulation fits a particular measurement of the reference run within an accuracy range of either 0.1 cm3 (2Q) for cumulative outflow 0.021 cm3 cm-3 for water content (2) or 2.10-3 m d-1 in terms of flux density (2F). For example, in Fig. 1B at 
, the cumulative outflow of the reference run equals 21.82 cm3, so all Monte Carlo simulations are accepted that have a cumulative outflow at 
that lie between 21.72 and 21.92 cm3. Hence, at , the Monte Carlo Runs 1, 2, and 3 are accepted and Monte Carlo Run 4 is rejected. At every point in time the measurement of the sandy soil is represented by a number of Monte Carlo simulations, with known soil hydraulic parameters. These vectors of accepted parameters are then used to calculate the standard deviation of each parameter.
PIMLI uses the following measure to quantify the sensitivity of a specific measurement for one of the model parameters
 | (8) |
(p)m is the standard deviation of parameter, p, in the accepted sets at measurement, m, and (p) is the standard deviation of the parameter, p, in the initial ranges. If ICm(p) is close to zero, this implies that this measurement is hardly sensitive for parameter, p, whereas an information content close to one indicates that the information content or sensitivity of measurement, m, for parameter, p, is very high. The time series of information content is then used to split the total time series of measurements into disjunctive subsets, with only measurements of high information content for a specific parameter. Hence, robust information for the least sensitive parameters will only appear if the most dominant parameters are close to their true value. Therefore, we followed an iterative procedure in which each iteration is used to identify a subset of measurements with highest information content for a not yet constrained parameter. Once these subsets are identified, then the histogram of the accepted parameters in this set is used for sampling in the next iteration. For this, the histogram is divided into 10 different equidistant classes in which the frequency of every class is used in a Latin Hypercube sampling.
To reduce the number of Monte Carlo simulations and still guarantee that enough simulations are accepted to ensure a representative value of the information content of a particular measurement, we temporarily adjusted Q to 5Q for the localization of the different subsets after the first 1000 Monte Carlo simulations. Every next iteration, Q was reduced with 0.05 cm3 until the final value of 0.05 cm3 for Q was achieved in the fourth iteration. Once all subsets are localized, then the parameters are more constrained in the sampling procedure. This leads to more accepted runs in the next iteration.
Normally, sensitivity analyses are used to calculate sensitivity coefficients that characterize the behavior of the objective function at a particular location in parameter space, presumably in the vicinity of the true parameter values (Kool and Parker, 1988; imnek and van Genuchten, 1996). In contrast to sensitivity analysis, PIMLI uses the variability in time of this sensitivity, expressed as the information content of a measurement with respect to the various model parameters, to select subsets of data that each contain explicit information for a particular parameter. Once identified, these subsets are used to constrain its corresponding parameter.
|
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
r–n and Ks– response surfaces indicate a correlation between both parameters. If r or Ks is perturbed positively, then n or will also be perturbed positively (Toorman et al., 1992). The relatively low sensitivity of the objective function for r is due to the fact that the effective saturation (Se) decreased only to about 0.08 within the range of measurements for the sandy soil used here. This makes the estimation of r dependent on extrapolation beyond this range. A similar problem was also encountered by Parker et al. (1985) and imnek et al. (1998). The relatively large area by the objective function = 0.10 contour in Fig. 2e and 2f indicates that there are many parameter combinations in the Ks– and r–Ks parameter plane that all provide reasonably good predictions for the cumulative outflow, flux density, and the average water content curve in Fig. 1A. If a conventional parameter optimization technique is used to simultaneously estimate s, r, , n, Ks, and based on the defined objective function, it cannot be successful because some of the response surfaces, like the Ks– and r–Ks parameter planes, do not display a well-defined global minimum. Hence, it would depend on the sensitivity of the optimizer and the parameter start values, and it is unlikely to find a unique parameter set. The experimental work by Eching and Hopmans (1993) and Eching et al. (1994) showed that these uniqueness problems could be minimized if both cumulative outflow and soil water pressure head data were included in the objective function.
|
s) of the soil sample was the most sensitive parameter and that there is a distinct difference in information content of the various sets of measurements for the same parameter. Moreover, at the beginning of the outflow experiment, the average water content of the soil sample contains the most information for the s parameter (Fig. 3B1). This is understandable because the water content is close to saturation at the beginning of the experiment and the sample is almost dry at the end. The distinct difference in information content of the average water content and cumulative outflow for the saturated water content is due to the fact that the cumulative outflow only provides information about a difference in water content and not about the absolute water content of the soil sample (i.e., cumulative outflow between 
and 
is comparable with cumulative outflow between 
). Simultaneous estimation of s and r on the basis of cumulative outflow only is therefore impossible (Van Dam et al., 1992).
|
s is identified (0 
t 0.5) and the histogram of the accepted parameters in this set is used for sampling in the next iteration, the information for the parameter at hydraulic equilibrium in the cumulative outflow measurements at the beginning of the experiment becomes more apparent (Fig. 3A2). Thus, a pressure step that passes the air-entry value of the soil sample contains most information for the parameter . The cumulative outflow measurements during hydraulic equilibrium after the first pressure increment (1.0 t 1.5) are therefore used to constrain in the next iteration. When the uncertainty in the hydraulic parameters s and diminishes, then the information for the parameter n appears at hydraulic equilibrium in the cumulative outflow at intermediate pressure steps in the next iteration (5.0 t 5.5). Although the error in water content is much larger than the error in outflow, the water content of the soil sample still provides reasonable information for the parameter n.
In Iteration 3, after constraining s, , and n in the sampling procedure, the information content for is highest immediately after pressure increments in the low water content range of the soil sample. This makes sense because the factor Se in the hydraulic conductivity function, Eq. [3], is most dominant at lower water contents. Hence, the flux density measurements between 
are used to constrain the parameter in the next iteration. Once the parameter is adjusted in the next iteration then the information of the residual water content (r) appears in the final cumulative outflow at the end of the experiment (18 t 20). Hence, if the shoulder and transition part of the water retention (Hollenbeck and Jensen, 1998) are constrained, then the final cumulative outflow of the soil sample provides a reasonable estimate for the parameter r. If the error in water content would be reasonably lower than 0.021 cm3 cm-3, the average water content of the soil sample can be used to further constrain r. Once the uncertainty in the water retention of the soil sample and are reduced, the information for the saturated conductivity (Ks) appears in the flux density out of the soil sample (Fig. 3C6) after the third suction step (2.0 t 2.25). In general, the information for the parameters and Ks is not found in the entire time series of cumulative outflow or water content, but only in the flux density measurements immediately after a suction increment. This is also demonstrated if we compare the response surfaces of Fig. 2e with the response surfaces that are calculated with the selected sets of measurements for and Ks in the objective function (Fig. 4)
. The mutual dependency of the parameters Ks and , as found in case of the entire time series, is reduced if we focus on disjunctive subsets. The sensitivity of the OF(b) to strongly increases and the identifiability of therefore increases (Fig. 4B). As a result of this the identifiability of Ks also increases. Although not presented, the r–n and r–Ks parameter planes yield identical results.
|
s, , and n, 10 iterations of 1000 Monte Carlo simulations are sufficient to accurately identify these flow parameters, whereas for Ks, , and r, more iterations still improve the results.
|
s, , n, Ks, and for a sandy and silty soil. The parameter values of the sandy and silty soil as used in the reference run lie within the 95% confidence interval as found with PIMLI. Moreover, standard deviations of parameter ranges are reduced with at least 75% for r (silt) to 99% for n (silt). Identification problems only occur for (silt) and various parameters of the clayey soil where the effective saturation only decreased to 0.70 and therefore the range of measurements is not sufficient to select truly disjunctive subsets of measurements. Parameters then remain too much correlated, as found with the response surfaces in Fig. 2d. The overestimation of r is then compensated by a slight overestimation of n. The relatively large standard deviation of for the clayey soil and partly the silty soil is due to the fact that most information for this parameter occurs at the dry end (see Fig. 3C4), also beyond the range of measurements of both soils.
|
log(-h) 2.7]. Divergence occurs only beyond the range of measurements. This inverse problem is therefore ill posed or nonidentifiable because more than a single parameter set leads to the same model response within the range of measurements of this MSO experiment. This points at the major limitation for parameter identification. Only unique parameter estimates can be obtained if the range of measurements is sufficient to divide the entire data set into disjunctive subsets that each contain the most information for the different unknown parameters. So, with this MSO experiment, we showed that there is enough information in the cumulative outflow, flux density, and water content to enable a unique parameter combination with PIMLI, if the range of measurements is large enough.
|
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
s, r, , n, Ks, and on the basis of cumulative outflow, water content, and flux density measured during a multi-step outflow experiment. Therefore, we propose PIMLI, a Parameter Identification Method based on the Localization of Information, that uses the variation in time of the model sensitivity for the various parameters to split the total set of measurements into disjunctive subsets that each contain the most information on one of the model parameters.
Results of PIMLI analysis showed that the average water content of the soil sample contains the most information for the saturated water content at the beginning of the outflow experiment. The distinct difference in information content of the average water content and cumulative outflow for the saturated water content is due to the fact that the cumulative outflow only provides information about a difference in water content and not about absolute water content in the soil sample. Simultaneous estimation of the saturated and residual water content on the basis of cumulative outflow only is therefore impossible (Van Dam et al., 1992). The cumulative outflow at hydraulic equilibrium after the pressure step that passes the air-entry value of the soil sample contains the most information for the parameter . Later in the experiment, also during hydraulic equilibrium, the information for n is highest. Once the water retention of the soil sample is more constrained, then the information for the parameters and Ks appears in the flux density of the outflow of the soil sample. The cumulative outflow measurements at the end of the experiment contain the most information for the residual water content. The information content for Ks and is highest immediately after a pressure step in the high and low water content range, respectively. PIMLI analysis also shows the limitations of experimental data of which the range of measurements is not large enough to split the data set into truly disjunctive subsets.
Using the different localized subsets in an iterative Monte Carlo simulation procedure, we showed, using PIMLI, that there is enough information in the cumulative outflow, water content and flux density to enable a unique set of soil hydraulic parameters (s, r, , n, Ks, and ), if the range of measurements is sufficient.
|
ACKNOWLEDGMENTS |
|---|
Received for publication September 17, 1999.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
imnek. 1999. Review of inverse estimation of soil hydraulic properties. p. 643–659. In M.Th. van Genuchten, et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. Univ. of California, Riverside.imnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour. Res. 32:2683–2696.imnek, J., O. Wendroth, and M.Th. van Genuchten. 1998. Parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.This article has been cited by other articles:
|
|
 |
|
  J. A. Vrugt, P. H. Stauffer, Th. Wohling, B. A. Robinson, and V. V. Vesselinov Inverse Modeling of Subsurface Flow and Transport Properties: A Review with New Developments Vadose Zone J., May 27, 2008; 7(2): 843 - 864. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 T. Wohling, J. A. Vrugt, and G. F. Barkle Comparison of Three Multiobjective Optimization Algorithms for Inverse Modeling of Vadose Zone Hydraulic Properties Soil Sci. Soc. Am. J., January 25, 2008; 72(2): 305 - 319. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. Mertens, R. Stenger, and G. F. Barkle Multiobjective Inverse Modeling for Soil Parameter Estimation and Model Verification Vadose Zone J., August 24, 2006; 5(3): 917 - 933. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
T. J. Heimovaara, J. A. Huisman, J. A. Vrugt, and W. Bouten Obtaining the Spatial Distribution of Water Content along a TDR Probe Using the SCEM-UA Bayesian Inverse Modeling Scheme Vadose Zone J., November 1, 2004; 3(4): 1128 - 1145. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
M. Javaux and M. Vanclooster In Situ Long-Term Chloride Transport through a Layered, Nonsaturated Subsoil. 1. Data Set, Interpolation Methodology, and Results Vadose Zone J., November 1, 2004; 3(4): 1322 - 1330. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
Z. F. Zhang, Z. F. Zhang, A. L. Ward, and G. W. Gee Estimating Soil Hydraulic Parameters of a Field Drainage Experiment Using Inverse Techniques Vadose Zone J., May 1, 2003; 2(2): 201 - 211. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. A. Vrugt and W. Bouten Validity of First-Order Approximations to Describe Parameter Uncertainty in Soil Hydrologic Models Soil Sci. Soc. Am. J., November 1, 2002; 66(6): 1740 - 1751. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. A. Vrugt, J. W. Hopmans, and J. Simunek Calibration of a Two-Dimensional Root Water Uptake Model Soil Sci. Soc. Am. J., July 1, 2001; 65(4): 1027 - 1037. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| The SCI Journals | Agronomy Journal | Crop Science | |||
| Journal of Natural Resources and Life Sciences Education |
Vadose Zone Journal | ||||
| Journal of Plant Registrations | Journal of Environmental Quality |
The Plant Genome | |||
