Vadose Zone Journal 2:201-211 (2003)
© 2003 Soil Science Society of America
Estimating Soil Hydraulic Parameters of a Field Drainage Experiment Using Inverse Techniques
Z. Fred Zhang*,
Andy L. Ward and
Glendon W. Gee
Pacific Northwest National Laboratory, Hydrology Group, MSIN K9-33, P.O. Box 999, Richland, WA 99352
* Corresponding author (fred.zhang{at}pnl.gov)
Received 7 March 2002.
 |
ABSTRACT
|
|---|
Accurate assessment of water flow and contaminant transport in unsaturated porous media at the field scale is often hindered by difficulties associated with obtaining reliable estimates of soil hydraulic properties. The unsteady drainage-flux method is one of the commonly used methods to measure in situ unsaturated hydraulic properties of soils. However, the properties obtained by this method using instantaneous profile data analysis may not be the best estimation of actual values of hydraulic properties. We present an improved analysis of the data from drainage experiments using inverse modeling, which uses nonlinear regression methods to estimate hydraulic parameters. Parameter identifiability is evaluated through sensitivity and uniqueness analyses. We used the combination of the inverse modeling program, UCODE, with the flow simulator, STOMP, for inverse modeling. Applying the inverse method to a field drainage experiment in sandy soil showed that all the van Genuchten (1980) hydraulic parameters could be estimated uniquely when both water content (
) and pressure head (h) data were used. The parameter estimates by inverse technique using both and h data simulated the flow better than the parameter values obtained by the conventional instantaneous-profile analysis method. After the spatial and temporal sensitivities were analyzed, a more rational experimental design was recommended.
|
INTRODUCTION
|
|---|
AN ADEQUATE PHYSICAL DESCRIPTION of water flow and contaminant transport in the vadose zone requires a priori knowledge of soil hydraulic property data. One of the field methods to determine unsaturated hydraulic properties in situ is the unsteady drainage-flux method (Green et al., 1986), which is based on the Darcian analysis of transient soil water content and hydraulic head profiles during vertical drainage following a heavy irrigation. The soil is wetted beyond the maximum depth for which the determinations are desired. The groundwater table should be sufficiently below this depth to obtain maximum possible unsaturated drainage. Richards et al. (1956) first used the drainage-flux method in the field. Nielsen et al. (1964), Rose et al. (1965), and van Bavel et al. (1968) developed it further. Watson (1966) improved the analysis of data by replacing the computation of differences in time and depth by the presumably more accurate instantaneous profile method. Many investigators have used the method to determine the hydraulic conductivity of well-drained soils. Fluhler et al. (1976) evaluated the error associated with hydraulic conductivity (K) values determined by this method. In their analysis, the relative errors were 20 to 30% of the measured K in the wet range but much greater in the drier range. Errors in K values obtained during the early stages of drainage are primarily due to errors in determining the hydraulic gradient, while at later times, errors in water content measurement are more serious. Falleiros et al. (1998) tested the repeatability of the instantaneous profile method and in their analysis found it unsatisfactory for describing soil water dynamics because of its spatial and temporal variability. The unsteady drainage-flux method requires that there be no evaporation occurring at the soil surface during the whole experiment. Due to the nature of the instantaneous profile method, the failure of even one tensiometer, which is not uncommon in the field, will significantly affect the calculation of hydraulic properties. Therefore, it may be useful to have a method that applies under flexible boundary conditions and can also give the optimized hydraulic property.
The solution of an inverse problem entails determining unknown causes on the basis of their effects (Hopmans and Simunek, 1999). The inverse method can be used to estimate hydraulic parameters by minimizing an objective function containing the sums of squared deviations between observations and predicted flow variables. Inverse methods have the benefits of (i) calculating parameter values that produce the best fit between observed and simulated values, (ii) quantifying the confidence limits in parameter estimates and the predictions, (iii) providing diagnostic statistics that quantify the quality of calibration and data shortcomings and needs, and (iv) not restricting the initial and boundary-flow conditions, the constitutive relationships, or the treatment of heterogeneity. The use of inverse methods for determining the hydraulic parameters of unsaturated soil from transient experiments was first reported by Zachmann et al. (1981). Dane and Hruska (1983) used an inverse method to estimate soil hydraulic parameter during drainage. Due to the lack of adequate analytical tools, only two parameters were estimated. Zijlstra and Dane (1996) inversely estimated the hydraulic parameters of homogeneous and layered soils using a quasi-Newton algorithm with Richardson extrapolation and generated water content data of various sizes and degrees of perturbation from drainage experiments in synthetic soils. The method resulted in successful simultaneous identification of up to nine parameters. Kool and Parker (1988) analyzed the inverse problem of determining unsaturated sol hydraulic properties from hypothetical one-dimensional transient infiltration and redistribution events. Their results show that simultaneous use of water contents and heads as input data for the inverse problem is partially redundant, although both types of data are required to obtain a reasonably well-posed problem. A number of laboratory and field applications (e.g., van Dam et al., 1992; Parkin et al., 1995; Simunek and van Genuchten, 1996; Lehmann and Ackerer, 1997; Abbaspour et al., 1997; Inoue et al., 2000; and Zhang et al., 2000) have shown the potential of the inverse method for improved designs and analyses of vadose zone flow and transport experiments. Hence, improved results may be obtained using the inverse method to analyze the data from a drainage experiment.
A requirement of using the inverse method is that the inverse problem not be ill posed. The ill-posedness is generally characterized by the instability and nonuniqueness of the identified parameters (Carrera and Neuman, 1986). In an instable solution small errors in the observations will cause serious errors in the parameter estimates. When nonuniqueness occurs, the criterion to be minimized has local minima or a global minimum at more than one point in the parameter space. The parameter estimates will differ according to the initial guess values of the parameters. If a problem is ill posed, some of the parameter estimates will be much different than the true values. Kool et al. (1985) and Parker et al. (1985) applied the inverse method to the laboratory-based one-step outflow method for determining hydraulic parameters. They concluded that nonuniqueness could be minimized if the experiment were designed to cover a wide range in water contents. Zhang et al. (2000) inversely estimated field hydraulic parameters using a surface line source. They found that measurements of only pressure head or only tracer travel time (T) did not give unique estimates of the saturated hydraulic conductivity, Ks, and the inverse macroscopic capillary length scale, 
. Any combination of the measurements of , h, or T gives a unique estimation of Ks and .
In this research, the data from a drainage experiment (Rockhold et al., 1988) were reanalyzed to estimate hydraulic parameters using an inverse method. The results were then compared with those obtained with conventional instantaneous profile analysis and with those measured in the lab. The inverse modeling was performed using the combination of the inverse model UCODE (Poeter and Hill, 1998, 1999) and the STOMP flow simulator (White et al., 1995; White and Oostrom, 2000). Detailed description of the two programs is given below. Parameter stability was examined through sensitivity analysis. Low sensitivity indicates high parameter instability. The parameter uniqueness was evaluated by the correlation coefficients between parameters. High correlation indicates nonuniqueness of parameter estimates.
|
MATHEMATICAL RELATIONS
|
|---|
The hydraulic parameter values are estimated by minimizing an objective function, O(ß), where ß denotes the parameter set to be estimated. The objective function is a measure of the fit between simulated values and observations and is typically defined as the sum of the weighted squared residuals:
 | [1] |
where yi is the ith observation, 
i is the ith simulated value, wi is the weight associated with the ith observation and is defined as the reverse of the variance of the measurement error, and N is the total number of observations. UCODE calculates the weights based on the statistics (i.e., variance, or standard deviation, or coefficient of variation of the error) of the observations (Hill, 1998).
Parameter identifiability can be evaluated through sensitivity and uniqueness analyses and two-dimensional response surfaces of the objective function. Parameter sensitivity measures how sensitive an observation is to the changes of a parameter. The dimensionless scaled sensitivities (SS) are defined as:
 | [2] |
where ßj denotes the jth parameter. The overall sensitivity of a parameter is described by the composite scaled sensitivities (CSS), which is the average of the SS values associated with the same parameter. Hill (1992), Anderman et al. (1996), and Hill (1998) calculated the composite scaled sensitivity for the jth parameter, CSSj, as
 | [3] |
Given the same observation error, a lower CSS value indicates a larger uncertainty of parameter estimate. According to Hill (1998), if one or more parameters have CSS values less than about 0.01 times the maximum CSS value, it is likely that the regression will not converge. For easy comparison of the CSS values associated with different parameters, a CSS ratio (
) for parameter j is defined as
 | [4] |
where max(CSS) is the maximum of the CSS values for all the parameters. Hence, the maximum value of is unity, which is associated with the parameter with the maximum CSS value. A parameter with very small value (e.g., <0.01; Hill, 1998) is not likely to be identifiable using corresponding observations. Different types of observations usually have different CSS values for the same parameter.
The above dimensionless sensitivities are needed to compare the importance of different types of observations or the overall effects of the observations with the estimation of parameter values. Calculating dimensionless scaled sensitivities requires that the weights associated with each observation be known. For other purposes, dimensional sensitivities are required, which do not need the predetermination of the weights associated with each observation. For example, dimensional sensitivities are useful to compare the importance of different parameters in predicting a single type of observations. Dimensional sensitivities also can be used to examine the spatial and temporal effects on parameter sensitivity. The dimensional 1%-scaled sensitivities (PSS) is defined as the change of the simulated value when a parameter value increases by 1%:
 | [5] |
After the PSS values at different positions and times are calculated, the distribution of the sensitivities can be contoured. The sensitivity maps then can be used to identify where additional observations would be most important to the estimation of different parameters. The spatial and temporal distribution of PSS can be used for determining the best experimental design.
Correlation of parameter estimates can be analyzed by using the variance–covariance matrix for the final estimated parameter values. Parameter uniqueness was evaluated by the correlation coefficient (R) between the parameters. A high correlation between parameters indicates that the parameters may not be identified uniquely. The absolute values of the correlation coefficients |R| 
0.95 indicate that the parameter values cannot be uniquely estimated with the observations used in the regression (Hill, 1998). One possible cause for the existence of a high correlation between parameters is that two parameters are mathematically correlated under some condition. For example, the saturated water content (s) and residual water content (r) are perfectly correlated if only h data are used in the inverse modeling. This problem may be overcome by taking observations of different types, including water content, pressure head, and water flux.
Parameter identifiability can also be evaluated in terms of two-dimensional response surfaces of O(ß) of soil hydraulic parameter pairs. The occurrence of a well-defined global minimum indicates the uniqueness of the parameter estimates.
The parameter uncertainty shows the precision of the parameter estimates and was expressed as a confidence interval (CI). A narrow CI indicates greater precision. If the model correctly represents the system, the CI also can be thought of as a measure of the probable accuracy of the estimate. During the inverse modeling, some of the parameters may be log-transformed as Bj = log(ßj), where Bj is the regressional parameter of physical parameter ßj. UCODE calculates the linear confidence interval (LCI) for each regressional parameter. Hence, the LCIs for the physical parameters can be calculated as follows:
 | [6a] |
 | [6b] |
where a parameter with a hat symbol above it denotes the best-fit value, t(N, 1.0 - 
/2) is the Student t statistic for N degrees of freedom and a significance level of , and s is the standard deviation. Equation [6b] means that, when a parameter is log-transformed, its LCI is expressed as the optimized value multiplied or divided by (x/÷) a factor, which has the minimum value of unity.
|
MATERIALS AND METHODS
|
|---|
Hydraulic Functions
The relationship between and h was described by the commonly used van Genuchten (1980) function:
 | [7] |
where is a fitting parameter that is inversely proportional to the air-entry pressure value, n is a parameter related to the width of pore-size distribution of the medium, and m is a constant that is commonly approximated by m = 1 - 1/n (van Genuchten, 1980). The Mualem (1976) model described the hydraulic conductivity function. Combining Eq. [7] with the Mualem (1976) model yields the hydraulic conductivity function as
 | [8] |
where Ks is the saturated hydraulic conductivity.
Computer Programs and Flow Simulation
Two computer programs were used to carry out inverse modeling. The forward model used to simulate water flow was the STOMP numerical simulator (White et al., 1995; White and Oostrom, 2000). STOMP is designed to solve a variety of nonlinear, multiple-phase, flow and transport problems for unsaturated porous media. The UCODE model (Poeter and Hill, 1998, 1999) was used to minimize the objection function and calculate all kinds of diagnostic statistics, such as sensitivity and correlation coefficients and confidence intervals of the parameter estimates. To perform inverse modeling, STOMP was coupled with and controlled by UCODE. UCODE also provides a number of other options, such as executing the application model at initial parameter values and calculating parameter sensitivities and response surfaces.
The procedures to estimate the hydraulic parameters are
- UCODE runs STOMP at the initial or updated parameter values.
- UCODE extracts the model predictions from the STOMP outputs and calculates O(ß).
- UCODE updates the parameter values using the modified Gauss–Newton method.
- Steps 1 to 3 are repeated until the user-specified convergence criterion is met.
Experiment and Parameter Inversion
A drainage experiment was conducted in a steel caisson at the U.S. Department of Energy's Hanford Site, WA (Rockhold et al., 1988). The caisson was 7.6 m long and 2.7 m in diameter. It was filled with a relatively uniform material, consisting of approximately 97% sand, 2% silt, and 1% clay. Soil water contents were measured with a neutron probe. Neutron-probe readings were taken every 0.3 m to a depth 2.4 m, with an additional reading at 0.45 m. The measurement volume of neutron probe varies with soil water content. Van Bavel et al. (1956) indicated that the resolution of a neutron probe ranges from about a 16-cm radius at saturation to 70 cm at near zero water content. Lack of high resolution makes it impossible to detect accurately any discontinuity or sharp change in water content gradient in a soil profile (McHenry, 1963). Depending on the distribution of , the average over the measurement volume may be higher or lower than the water content at the center of it. Since the data were treated as point observations in inverse modeling, the differences of the average over measurement volume and the values at the center will have the same effect as observation errors on parameter estimation. At the drainage stage in a uniform soil, since the distribution of is nearly symmetric at the center of the measurement volume, we assume that the water content measured using a neutron probe represents the value at the center of the measurement volume. Soil water pressure heads were measured with tensiometers and a Tensimeter pressure transducer (Soil Measurement Systems, Tucson, AZ). Tensiometer readings were taken every 0.3 m to a depth of 1.8 m, with two additional readings at 0.15 and 0.45 m. Rockhold et al. (1988) calculated the hydraulic functions using the instantaneous profile method (Green et al., 1986). They also measured the hydraulic functions of the same soil in the laboratory by measuring the steady-state and h at multiple constant water fluxes on undisturbed soil samples.
Three cases were considered to invert the hydraulic parameters through the field drainage experiment: (i) using observations only, (ii) using h observations only, and (iii) using both and h observations. During the inverse modeling, the spatial step size was 0.01 m. When only the water content data were used, the boundary condition at soil depth z = 0 was zero flux and that at z = 1.80 m was unit gradient. When the simulation included h data, the values of h at z = 0.15 and 1.80 m were set as time-variable boundary conditions.
In the inverse procedure, the convergence criterion was set to 0.01, which means the nonlinear regression converges if the relative changes of the O(ß) are less than 0.01 for three sequential iterations. For the inverse procedure, Ks, , and n were log-transformed as this enhances the rate of convergence and constrains the solution to a positive parameter space (Carrera and Neuman, 1986). The error of measured with a neutron probe is approximately 0.01 m3 m-3. The error h using a Tensimeter pressure transducer is about 0.02 m (Marthaler et al., 1983). Hence, the standard deviation (SD) of 0.01 m3 m-3 is assumed for the observations and 0.02 m for the h observations. Although the assignment of the observation errors are subjective, in most circumstances the estimated parameter values are not very sensitive to moderate changes in the weights used (Hill, 1998). Our inverse simulations using different guessed standard deviations of and h show that for the infiltration experiment mentioned above, the changes of the optimized parameters are no more than 5% if the SD value of either or h changed by 50%.
Any physically reasonable guessed values of the parameters to be estimated may be used as initial values of the inverse model. For a well-posed inverse problem, different initial guess values will produce the same optimization. However, the inverse model generally converges faster if the initial guessed values are more close to the actual parameter values. We determined the initial values of the parameters on the basis of the results of Rockhold et al. (1988), who determined the hydraulic properties using the instantaneous profile method. The inverse model was run multiple times using different initial parameter values. The hysteresis and air entrapment were ignored.
For comparison with the inverse method, we also determined the hydraulic parameter values by simultaneously fitting K() and K(h) to corresponding lab measured values or the field hydraulic data determined using the instantaneous profile method.
|
RESULTS AND DISCUSSIONS
|
|---|
Parameter Identifiability
We used the composite scaled sensitivity (CSS) to evaluate the sensitivities of flow to changes of each parameter. According to Hill (1998), if one or more of the CSS ratio () values is less than about 0.01, it is likely that the nonlinear regression will not converge. However, we found that the nonlinear regression may or may not converge if one or more of the j values is between 0.01 and 0.1. Thus, we set the sensitivity ratio value of 0.1 as the criterion to determine the identifiable and unidentifiable parameters. Parameters with < 0.1 were treated as unidentifiable parameters.
The composite scaled sensitivity values () of flow to the hydraulic parameters are depicted in Fig. 1. Figure 1a shows that soil water content is most sensitive to the changes in s and least sensitive to changes in . The value for parameter was much less than 0.1 and indicated that is not identifiable using water content data only. The low sensitivity of to changes in was also noted by Dane and Hruska (1983) and Zijlstra and Dane (1996). Soil water pressure head was most sensitive to (Fig. 1b). The values for both s and r were less than 0.1 and indicated that the two parameters are not identifiable using pressure head data only. When and h data were combined in the objection function, the composite scaled sensitivity values of and h for the same parameter compensated each other (Fig. 1c). For example, the low sensitivity of was compensated by the high sensitivity of h to the changes in . The low sensitivities of h were compensated by the high sensitivities of to the changes of s and r. Thus, when both and h data were used, the difference between the maximum and minimum values of CSS was smaller than any of the cases when only one data type was used. When both and h data were used, none of the values was less than 0.1; hence, all the parameters were identifiable if no high correlations between parameters existed.
After sensitivity analysis, the parameters with values less than 0.1 were set to be constant during inverse modeling. The remaining parameters were estimated by nonlinear regression. The uniqueness of the parameter estimates was evaluated by the correlation coefficients (R) between the parameters (Table 1) at the optimized parameter values. A smaller value of |R| indicates the higher uniqueness. When only data were used in the inverse modeling, all the parameters except were estimated. The largest value of |R| was 0.771, which is much smaller than the critical value of 0.95 (Hill, 1998) and suggests that all the parameter estimates were unique. When only h data were used in the inverse modeling, parameters Ks, , and n were estimated. However, a high correlation existed between Ks and , with R = 0.968. This means that the estimates of Ks and were not unique. Thus, parameters Ks and were not identifiable using only h data, although h was very sensitive to both of the parameters. When both and h data were used in the inverse modeling, all the five parameters were estimated. The largest value of |R| was 0.803, which suggests that all the five parameter estimates were unique.
View this table:
[in this window]
[in a new window]
|
Table 1. Correlation coefficient between parameters when different types of data were used to inversely estimate the hydraulic parameter.
|
|
Considering both the sensitivity and uniqueness, for the drainage experiment conducted at the Hanford Site, parameter was not identifiable when only data were used, four of the five parameters (i.e., Ks, , s, and r) were not identifiable when only h data were used, and all the five parameters were identifiable when both the and h data were used. The results are consistent with findings by Kool and Parker (1988) in that both data types were required to obtain a reasonably well-posed problem.
Since the sensitivity coefficients are dependent on the parameter values, it is suggested that the sensitivity for all the parameters be reexamined at the final parameter values. The values of R depend on the total number of observations and observation error. A high correlation may occur if the number of observations is small and/or the observation errors are large.
The parameter identifiabilities were also examined using two-dimensional response surfaces of the objective function (Eq. [1]). There were three types of typical response surfaces:
- The contour lines are nearly parallel to one of the axes. This indicates that the observations are not sensitive to the changes in the parameter whose associated axis is nearly parallel to the contour lines. For example, soil water content was insensitive to the changes in according to the sensitivity analysis; thus, the value of O(ß) showed little change as changed (Fig. 2a). With this type of response surface, it is likely that the inverse modeling will not converge.
- The contour lines are nearly parallel but not parallel to any of the two axes. This indicates that the two parameters are highly correlated, while the flow is still sensitive to either of the parameters. Figure 2b is an example that shows the two-dimensional response surface at the log(Ks)–log() plane when only h data were used. Although pressure head is quite sensitive to both Ks and (Fig. 1b), there is not a well-defined minimum in the response surface due to the high correlation between the two parameters. With this type of response surface, the inverse modeling may converge, but the convergence is dependent on the starting values. Hence, the associated parameters are not identifiable.
- A minimum point is well defined, and the contours are nearly concentric ellipses. This indicates that the flow is sensitive to both of the parameters, and the correlation between the parameters is also low, and hence both of the parameters are identifiable. Figure 2c shows an example of the response surface at the log(Ks)–log() plane when both the and h data were used in the inverse modeling.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 2. The two-dimensional response surfaces in the log(Ks)– log() planes using (a) only , (b) only h, and (c) and h observations. Parameter Ks is in centimeters per hour, and is per centimeter.
|
|
Thus, the parameter sensitivity and uniqueness can be shown visually by response surfaces. However, since the influence of only two parameters on the O(ß) is investigated, the resulting response surface represents only a cross section of the full parameter space. By limiting the number of parameters within a single response surface analysis, the behavior of the O(ß) in the different parameter planes can only suggest how the O(ß) might behave in the full parameter space. For example, local minima of the O(ß) could exist, but not appear in the cross-sectional planes (Simunek and van Genuchten, 1996).
Vrugt et al. (2001) proposed a parameter identification method based on the localization of information (PIMLI), which uses the variability in time of the model sensitivity for the various parameters to split the total set of measurements into disjunctive subsets, each of which contains the most information on one of the model parameters. Each distinguished subset is used to constrain its corresponding parameter. Their results based on the synthetic data show 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. However, as pointed out by the authors, since PIMLI only uses a relatively small number of measurements with highest information content for a specific parameter, problems might occur if PIMLI is applied to real laboratory or field data sets that are corrupted with errors. Comparing with the approaches based on the Gauss–Newton method (e.g., UCODE), the PIMLI method requires much more time to obtain the optimized parameters since it needs 1000 Monte Carlo simulation for each iteration. For example, there is a simulation and one forward run of it takes 10 min. Five parameters need to be estimated and the parameter estimation needs 10 iterations to converge. Then, the total runs are 10(5 + 1) = 60, and it will take about 10 h if the Gauss–Newton method is used. When the PIMLI method is used, the total runs are 10 x 1000 = 10 000, and it takes 69 d, which is 166 times that of the Gauss–Newton method. Hence, considering the simulation time, it may be difficult to apply the PIMLI method to many field experiments.
Parameter Estimates and Comparisons between Simulations and Observations
The optimized parameter values with 95% linear confidence interval (LCI) using five methods are shown in Table 2. Three of the five methods (i.e., 1, 2, and 5) gave the estimates of all the five parameters. Methods 3 and 4 can only estimate some of the parameters because of the high correlation between parameters and/or low parameter sensitivities. Among all the methods, the inverse method using only h data (Method 4) produced the poorest results since four of the five parameters could not be identified. The inverse method using both and h data (Method 5) gave the best results shown by the smallest 95% LCI for the estimates of parameters Ks, , and s and relatively small 95% LCI for n and r.
View this table:
[in this window]
[in a new window]
|
Table 2. The optimized parameter values with 95% linear confidence interval (LCI). Parameters Ks, , and n were log-transformed when they were estimated.
|
|
The optimized values of s from different methods are similar, ranging from 0.286 to 0.335. For the remaining parameters, the difference between values estimated by the different methods was, on average, larger than the uncertainty (95% LCI) of the estimates. Note that the parameters that were kept constant and those with nonunique estimates were not considered when the comparisons were made. For example, the 95% LCI of parameter Ks given by different methods varies from 17 to 53%, while the differences between the optimized values of Ks vary from 16% (Methods 1 and 2) to 108% (Methods 2 and 5). The 95% LCI of parameter given by different methods varies from 5 to 13%, while the differences between the optimized values of vary from 4% (Methods 2 and 5) to 52% (Methods 1 and 2). The 95% LCI of parameter n varies from 5 to 52%, and the differences between the optimized values of n vary from 4% (Methods 3 and 4) to 126% (Methods 1 and 3). For parameter r, the uncertainty of the estimates is no larger than ± 0.004 m3 m-3, while the best-fit estimates vary from 0.065 to 0.102 m3 m-3. Therefore, selecting an appropriate method is more important than trying to improve the resolution of an inappropriate method. For example, to predict the flow in the field, the field instantaneous profile method may be more appropriate than the lab method that measures hydraulic parameters using small (e.g., 5-cm diam.) soil cores, although the lab method may be better controlled in flow boundary conditions than the field method.
Comparisons of the observed values and simulated values using these parameters are shown in Fig. 3. Using the parameters measured in the laboratory (Method 1, Table 2) significantly overestimated (aFig. 3-1a) and underestimated h (bFig. 3-1b). When the parameter values measured using the instantaneous profile method (Method 2), the simulations of both and h were improved significantly (bFig. 3-2a and 3-2b), although values were overestimated at the medium water content range (e.g., between about 15 and 25%). In Method 3, four of the five parameters were inversely estimated using data only. Consequently, the simulation of matched the observations well (Fig. 3-3a), while the simulations of h were significantly different from the observations except when the soil was very wet (e.g., h > -10 cm) (Fig. 3-3b). Similarly, Method 4 used the h data only to inversely estimate the parameters, and only parameter n may be uniquely determined. As a result, the pressure head was simulated well (bFig. 3-4b), while water content was simulated poorly (aFig. 3-4a). Method 5 used both and h data to inversely estimate the parameters, which were all uniquely determined (Table 2). Using these parameter values, both and h were simulated well (bFig. 3-5a and 3-5b).
View larger version (40K):
[in this window]
[in a new window]
|
Fig. 3. Comparisons between the observed and the simulated values using the parameter values determined by different methods. The plot numbers are corresponding to Methods 1 through 5, and the parameter values are listed in Table 2.
|
|
In summary, when the parameter values were obtained by fitting to laboratory observations, both the field and h were poorly simulated. The parameter values obtained by fitting to the results of the instantaneous profile method simulated the flow reasonably well. The parameter values obtained by the inverse method using both and h data produced the best simulation of water flow, which indicates an improved analysis of data from a drainage experiment over the traditional instantaneous profile method.
Spatial and Temporal Sensitivities and Experimental Design
To accurately estimate the hydraulic parameters, an experiment should be designed in such a way that observations have relatively large sensitivity to changes in parameters. Hence, it is suggested that sensitivity and uniqueness analyses be performed before the design of an experiment for inverse modeling if similar analyses cannot be found from references. The sensitivity of flow to changes in a hydraulic parameter is a function of both soil depth (z) and time (t) and is different for a different parameter. Figure 4 shows the contours of the 1%-scaled sensitivity (PSS) of to changes in the hydraulic parameters in the t–z plane. The absolute values of the 1%-scaled sensitivity of to changes in Ks, n, and s increase with z (Fig. 4a, 4c, and 4d), while the sensitivity to changes in and r decrease with z (Fig. 4b and 4e). These effects of z are more pronounced at smaller times than at larger times. The sensitivities to changes in all the parameters except r decrease with t (Fig. 4). These observations suggest that, for accurate estimation of all the parameters, should measured under both wet and dry soil conditions, since the soil at the upper profile is drier than that at the lower profile and at larger times than that at smaller times.
The spatial and temporal effects on |PSS| of h are similar for all the parameters (Fig. 5). The values of |PSS| of h decrease with z and increase with t. This suggests that h observations should be taken at relatively dry soil conditions to estimate parameters more accurately. Tensiometer response time may vary from minutes to hours or even days (Klute and Gardner, 1962) depending on the conductance of the porous cup and the conductivity and pressure head of the surrounding soil. Therefore, considering both the low sensitivity at very wet soil conditions and the equilibrium time needed for tensiometers, one may measure only at the early stage (e.g., within a few hours since the start of the drainage process) of a drainage experiment since the soil water content is high. The start of h measurements may be taken when there is enough time for tensiometers to reach equilibrium.
The above analyses of the sensitivities of and h to changes in soil hydraulic parameters suggest two approaches for experimental design. One is to measure over relatively large vertical scale (e.g., 1 m or larger) and at numerous soil depths such that the soil profile will include both dry and wet soil conditions during the drainage process. The measurements of h may be taken at multiple positions at the upper part of the soil profile that is relatively dry compared with the lower part. In this way, the overall experimental time may be kept relatively short. This approach is suitable for soils with a deep and uniform unsaturated profile. The second approach is to measure for relatively long duration but on a relatively small vertical scale. Pressure head can be measured at relatively larger times. This approach is suitable for soils with a relatively shallow profile or within a single layer for layered soils. The analysis of Zijlstra and Dane (1996) with error-free data of shows that, for simultaneous estimates of , Ks, and n using their proposed quasi-Newton method, the data with an extent of 10 depths and five times were needed. Our first approach is consistent with their results.
The application of the two approaches was tested using part of the observations (Rockhold et al., 1988) of the drainage experiment. Table 3 shows the parameter estimates using and h observations at early times (t 
7.6 h) or at shallow depths (z 0.45 m). The results are similar to those when all the data of and h are used (Method 5, Table 2).
|
CONCLUSIONS
|
|---|
The unsteady drainage method has been used to measure in situ unsaturated hydraulic properties. However, the properties obtained by the instantaneous-profile data-analysis method may not be the best estimation of the actual values. The inverse method has advantages of calculating parameter values that minimize the difference between observed and simulated values and quantifying the confidence limits in parameter estimates. We have presented an improved analysis of the data from a drainage experiment using inverse modeling, which uses nonlinear regression methods to estimate hydraulic parameters. Inverse modeling was performed using the combination of the inverse model UCODE and the STOMP flow simulator. Application of inverse modeling to a field drainage experiment showed that some of the hydraulic parameters could not be estimated inversely when only water content data or only pressure head data were used because of low sensitivity and/or high correlation between parameters. All the parameters could be inversely estimated uniquely when both water content and pressure head data were used. Comparisons of parameter estimates showed that, for a given parameter, the difference between parameter values obtained by different methods was larger than the uncertainty of the estimates. Hence, selecting an appropriate method for parameter estimation is more important than trying to improve the resolution of an inappropriate method. The spatial and temporal sensitivity analyses showed that, for accurate parameter estimation, water content should be measured under both wet and dry soil conditions, and pressure head measured under relatively dry soil condition. Hence, two approaches are recommended for a drainage experiment. One is to measure water content over a relatively large vertical scale and at numerous soil depths, with pressure head measured at multiple positions at the upper part of the soil profile that is drier than the deeper part. The second approach is to measure water content for relatively longer duration, but on a relatively small vertical scale; pressure head is measured at relatively larger times.
|
ACKNOWLEDGMENTS
|
|---|
Funding for this research is provided by the U.S. Department of Energy as part of the Hanford Science and Technology Vadose Zone Initiative. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DE-AC06-76RL01830.
|
REFERENCES
|
|---|
- Abbaspour, K.C., M.Th. van Genuchten, R. Schulin, and E. Schlappi. 1997. A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters. Water Resour. Res. 33:1879–1892.
- Anderman, E.R., M.C. Hill, and E.P. Poeter. 1996. Two-dimensional advective transport in ground-water flow parameter estimation. Ground Water 34:1001–1009.[ISI]
- Carrera, J., and S.P. Neuman. 1986. Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, stability, and solution algorithms. Water Resour. Res. 22:211–227.
- Dane, J.H., and S. Hruska. 1983. In-site determination of soil hydraulic properties during drainage. Soil Sci. Soc. Am. J. 47:619–624.[Abstract/Free Full Text]
- Falleiros, M.C., O. Portezan, and J.C.M. Oliveira. 1998. Spatial and temporal variability of soil hydraulic conductivity in relation to soil water distribution, using an exponential model. Soil Tillage Res. 45:279–285.
- Fluhler, H., M.S. Ardakani, and L.H. Stolzy. 1976. Error propagation in determining hydraulic conductivities from successive water content and pressure head profiles. Soil Sci. Soc. Am. J. 40:830–836.[Abstract/Free Full Text]
- Green, R.E., L.R. Ahuja, and S.K. Chong. 1986. Hydraulic conductivity, diffusivity, and sorptivity of unsaturated soils: Field methods. p. 771–798. In A. Klute (ed.) Method of soil analysis. Part 1. 2nd ed. SSSA Book Ser. 5. ASA and SSSA, Madison, WI.
- Hill, M.C. 1992. A computer program (MODFLOWP) for estimating parameters of a transient, three-dimensional, ground-water flow model using nonlinear regression. U.S. Geological Survey Open-File Rep. 91-484, 358. U.S. Geol. Survey, Reston, VA.
- Hill, M.C. 1998. Methods and guidelines for effective model calibration. U.S. Geological Survey Water-Resources Investigations Rep. 98-4005. Available at: http://water.usgs.gov/software/ucode.html [modified 4 Jan 2002; verified 14 Feb. 2003]. U.S. Geol. Survey, Reston, VA.
- Hopmans, J.W., and J. Simunek. 1999. Review of inverse estimation of soil hydraulic properties. p. 643–660. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media, Part 1. University of California, Riverside, CA.
- Inoue, M., J. Simunek, S. Shiozawa, and J.W. Hopmans. 2000. Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments. Adv. Water Res. 23:677–688.
- Klute, A., and W.R. Gardner. 1962. Tensiometer response time. Soil Sci. 93:204–207.
- Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem for transient unsaturated flow. Water Resour. Res. 24:817–830.
- Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties from one-step experiments by parameter estimation. I. Theory and numerical studies. Soil Sci. Soc. Am. J. 49:1348–1354.[Abstract/Free Full Text]
- Lehmann, F., and P. Ackerer. 1997. Determining soil hydraulic properties by inverse method in one-dimensional unsaturated flow. J. Environ. Qual. 26:76–81.[Abstract/Free Full Text]
- Marthaler, H.P., W. Vogelsanger, F. Richard, and P.J. Wierenga. 1983. A pressure transducer for field tensiometers. Soil Sci. Soc. Am. J. 47:624–627.[Abstract/Free Full Text]
- McHenry, J.R. 1963. Theory and application of neutron scattering in the measurement of soil moisture. Soil Sci. 95:294–307.
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
- Nielsen, D.R., J.M. Davisson, J.W. Biggar, and R.J. Miller. 1964. Water movement through Panoche clay loam soil. Hilgardia 35:491–506.
- Parker, J.C., J.B. Kool, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties from one-step experiments by parameter estimation. II. Experimental studies. Soil Sci. Soc. Am. J. 49:1354–1359.[Abstract/Free Full Text]
- Parkin, G.W., R.G. Kachanoski, D.E. Elrick, and R.G. Gibson. 1995. Unsaturated hydraulic conductivity measured by time domain reflectometry under a rainfall simulator. Water Resour. Res. 31:447–454.
- Poeter, E.P., and M.C. Hill. 1998. Documentation of UCODE, a computer code for universal inverse modeling. U.S. Geological Survey Water-Resources Investigations Rep. 98–4080. Available at: http://water.usgs.gov/software/ucode.html [modified 4 Jan 2002; verified 14 Feb. 2003]. U.S. Geol. Survey, Reston, VA.
- Poeter, E.P., and M.C. Hill. 1999. UCODE, a computer code for universal inverse modeling. Comput. Geosci. 25:457–462.
- Richards, L.A., W.R. Gardner, and G. Ogata. 1956. Physical processes determining water loss from soil. Soil Sci. Soc. Am. Proc. 20:310–314.
- Rockhold, M.L., M.J. Fayer, and G.W. Gee. 1988. Characterization of unsaturated hydraulic conductivity at the Hanford site. PNL-6488. Pacific Northwest Laboratory, Richland, WA.
- Rose, C.W., W.R. Stern, and J.E. Drummond. 1965. Determination of hydraulic conductivity as a function of depth and water content for soil in situ. Aust. J. Soil Res. 3:1–9.
- Simunek, 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.
- van Bavel, C.H.M., G.B. Stirk, and K.J. Brust. 1968. Hydraulic properties of a clay loam soil and the field measurement of water uptake by roots: 1. Interpretation of water content and pressure profiles. Soil Sci. Soc. Am. Proc. 32:310–317.
- van Bavel, C.H.M., N. Underwood, and R.W. Swanson. 1956. Soil moisture measurement by neutron moderation. Soil Sci. 82:29–41.
- van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1992. Inverse method for determining soil hydraulic functions from one-step outflow experiments. Soil Sci. Soc. Am. J. 56:1042–1050.
- 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. 2001. Information content of data for identifying soil hydraulic parameters from outflow experiments. Soil Sci. Soc. Am. J. 65:19–27.[Abstract/Free Full Text]
- Watson, K.K. 1966. An instantaneous profile method for determining the hydraulic conductivity of unsaturated porous materials. Water Resour. Res. 2:709–715.
- White, M.D., and M. Oostrom. 2000. User's guide of STOMP—Subsurface transport over multiple phases. PNNL-12034 UC-2010. Pacific Northwest National Laboratory, Richland, WA.
- White, M.D., M. Oostrom, and R.J. Lenhard. 1995. Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator. 1. Nonvolatile three-phase model description. Adv. Water Resour. 18:353–364.
- Zachmann, D.W., P.C. DuChateau, and A. Klute. 1981. The calibration of the Richards flow equation for a draining column by parameter identification. Soil Sci. Soc. Am. J. 45:1012–1015.[Abstract/Free Full Text]
- Zhang, Z.F., R.G. Kachanoski, G.W. Parkin, and B. Si. 2000. Measuring hydraulic properties using a line source: II. Field test. Soil Sci. Soc. Am. J. 64:1563–4569.[Abstract/Free Full Text]
- Zijlstra, J., and J.H. Dane. 1996. Identification of hydraulic parameters in layered soils based on a quasi-Newtonian method. J. Hydrol. 181:233–250.
This article has been cited by other articles:
|
|
|
|
 
A. L. Ward and Z. F. Zhang
Effective Hydraulic Properties Determined from Transient Unsaturated Flow in Anisotropic Soils
Vadose Zone J.,
November 20, 2007;
6(4):
913 - 924.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
H. Vereecken, R. Kasteel, J. Vanderborght, and T. Harter
Upscaling Hydraulic Properties and Soil Water Flow Processes in Heterogeneous Soils: A Review
Vadose Zone J.,
January 24, 2007;
6(1):
1 - 28.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
M. J. Singleton, E. L. Sonnenthal, M. E. Conrad, D. J. DePaolo, and G. W. Gee
Multiphase Reactive Transport Modeling of Seasonal Infiltration Events and Stable Isotope Fractionation in Unsaturated Zone Pore Water and Vapor at the Hanford Site
Vadose Zone J.,
August 1, 2004;
3(3):
775 - 785.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
