Soil Science Society of America Journal 66:1409-1423 (2002)
© 2002 Soil Science Society of America
DIVISION S-1—SOIL PHYSICS
Estimating Hydraulic Properties of a Fine-textured Soil Using a Disc Infiltrometer
R. C. Schwartz* and
S. R. Evett
USDA-ARS, Conservation and Production Research Lab., P.O. Drawer 10, Bushland, TX 79012-0010
* Corresponding author (rschwart{at}cprl.ars.usda.gov)
 |
ABSTRACT
|
|---|
Inverse optimization of parameters offers an economical means to infer soil hydraulic properties from in situ measurements of infiltration. We evaluated optimization strategies to inversely estimate soil hydraulic parameters using field measured tension disc infiltrometer data. We estimated the parameters n, 
, and Ks of the van Genuchten-Mualem (VGM) model, and a piecewise representation of conductivity near saturation using a numerical inversion of Richards' equation. In addition to cumulative infiltration, optimizations included in the objective function water retention data, water contents from cores extracted after termination of infiltration, or transient measurements of water contents using time domain reflectometry (TDR) probes. Three-parameter fits to field data were nonunique because of a positive correlation between and Ks. In contrast, fits of n and Ks with estimated from separate fits to retention data improved parameter identifiability while not compromising the fit to measured infiltration. Inverse optimizations that included in the objective function both water retention and cumulative infiltration, led to excellent fits of this data when initial volumetric water contents were >0.23 cm3 cm-3. Close fits to cumulative infiltration were also obtained at lower water contents, however, water retention data was underestimated likely because of hysteresis. Optimizations of cumulative infiltration with final soil core water content or TDR data led to estimates of final water contents that closely approximated measured water contents. However, measured TDR water contents were poorly matched by simulations at early times. A piecewise loglinear interpolation of hydraulic conductivity near saturation improved fits to measured cumulative infiltration and water retention data as compared with using the VGM model at all pressure heads.
Abbreviations: TDR, time domain reflectometry • VGM, van Genuchten-Mualem
|
INTRODUCTION
|
|---|
PARAMETERS DERIVED from in situ measurements of soil hydraulic properties are crucial to understanding and describing the dynamic processes of water flow in the field. The tension disc infiltrometer (Perroux and White, 1988) has become a valuable tool to investigate the hydraulic properties of soils at or near the surface. This infiltration-based method is particularly suitable for quantifying changes in near surface hydrology resulting from soil management activities such as tillage (Sauer et al., 1990; Logsdon et al., 1993). Although unconfined flow below the infiltrometer disc complicates the analyses of infiltration measurements, various methods have been devised to infer hydraulic properties from disc infiltrometer measurements. These techniques are based on quasi-analytical solutions of transient flow at early times (e.g., Smettem et al., 1994), Wooding's (1968) analysis of steady state infiltration from a disc source (Ankeny et al., 1991; Logsdon and Jaynes, 1993; Hussen and Warrick, 1993), or by inverse parameter optimization of the axisymmetric form of Richards' equation (
im
nek and van Genuchten, 1996). Angulo-Jaramillo et al. (2000) discuss many of the difficulties associated with both the transient and steady state analysis of unconfined infiltration. The focus of this paper is the estimation of soil hydraulic properties through inverse parameter optimization of the governing equations that describe water flow from a disc source. Inverse procedures tend to be less restrictive than direct analysis using quasi-analytical solutions and have the potential to yield information about conductivity and water retention over a wide range in pressure heads from a single infiltration experiment.
imnek and van Genuchten (1996)(1997) proposed an inverse method to estimate hydraulic properties using cumulative infiltration data from a disc infiltrometer. Based on the results of inverse fits to numerically generated data, they concluded that identifiability of parameters is improved when other information is included in the objective function. The most promising scenario was an objective function that included initial and final water contents as well as cumulative infiltration data. Final water contents were assumed to be in equilibrium with the supply pressure head and taken at the soil surface upon the termination of infiltration experiments. imnek et al. (1998b) later used this method in conjunction with multiple tension infiltrometer data to estimate hydraulic properties of two field soils.
Despite the advantages of using inverse optimization in conjunction with disc infiltrometer measurements, these methods are hampered by a number of practical problems that must be overcome so that they can be successfully used in the field. Presently, only a few researchers have described inverse optimizations of disc infiltrometer measurements in the field (imnek et al., 1998a, 1998b), likely because of the difficulty of obtaining and incorporating meaningful auxiliary data along with cumulative outflow data. For instance, errors can arise in the determination of volumetric water contents when sampling the soil surface after the removal of the disc infiltrometer because of the small sampling depth required, and because bulk density must be estimated for this thin layer (Angulo-Jaramillo et al., 2000). The most pertinent soil volume of interest directly beneath the disc is typically inaccessible to sensors. Steep gradients in water content near the soil surface require an accurate estimate of the initial water content profile for inverse estimation methods. Lastly, water retention curves obtained through numerical inversion of field infiltration experiments have typically compared poorly with laboratory retention data (imnek et al., 1998b; imnek et al., 1999b). A consistently workable method for combining inverse parameter estimation with field measured data is still elusive.
|
THEORY
|
|---|
Governing Equations
Isothermal water flow for a radially symmetric two-dimensional region in nonswelling, homogeneous, isotropic soils can be described with the following form of Richards' equation (Warrick, 1992):
 | [1] |
where 
is the volumetric water content (cm3 cm-3), t is time (s), z is the vertical coordinate taken positive downwards (cm), K is hydraulic conductivity (cm s-1), h is the pressure head (cm), and r is the radial coordinate (cm). Equation [1] can be solved subject to an initial water content depth profile i(z)
 | [2] |
and boundary conditions Eq. [3a] to [3e]
 | [3a] |
 | [3b] |
 | [3c] |
 | [3d] |
 | [3e] |
where h0 is the inlet pressure head at the soil surface. Equation [3a] is a prescribed head surface boundary below the disc source with radius r0 and Eq. [3b] describes a zero flux boundary at the surface for r > r0. The lower boundary condition, Eq. [3c], permits free drainage at an effectively infinite distance from the source and Eq. [3d] and [3e] specify zero flux boundaries. The radial flux term in Eq. [1] is indeterminate at r = 0 and must be transformed to apply the boundary condition [3e]. Application of l'Hospital's rule to the radial flux term and Eq. [3e] gives
 | [4] |
The right-hand side expression in Eq. [4] was applied at r = 0 to implement the zero flux boundary for a numerical solution to Eq. [1].
The VGM model (Mualem, 1976; van Genuchten, 1980)
 | [5] |
 | [6] |
can be used to describe the constitutive soil hydraulic properties of Eq. [1] at pressure heads less than hp0. Here r and s are the residual and saturated water contents (cm3 cm-3), respectively, Ks is the saturated hydraulic conductivity (cm s-1), S is the fluid saturation ratio 
/
, m = 1 - 
, and n and (cm-1) are empirically fitted parameters. Analogous to imnek and van Genuchten (1997), Ks is considered as a fitted parameter that may differ substantially from the true saturated conductivity. Moreover, at pressure heads very near saturation, K(h) for fine-textured soils is overestimated by Eq. [6] when fitted to unsaturated conductivity data. Likewise, unsaturated conductivity for fine-textured soils is underestimated by Eq. [6] when Ks is forced to match measured values during parameter estimation (Assouline and Tartakovsky, 2001). Consequently, K(h) must be modified near saturation to correctly describe infiltration into dry fine-textured soils. At pressure heads greater than hp0, K(h) can be described using piecewise continuous loglinear interpolation
 | [7] |
where L0, L1, L2, and L3 are the Lagrangian coefficients for linear interpolation and hp0, hp1, hp2, and hp3 are monotonically increasing pressure heads for which K(h) is known or can be estimated.
Steady State Flow from a Disc Source
Wooding (1968) demonstrated that by linearization of the governing partial differential equation steady state outflow Q(h) (cm3 s-1) from a circular source at a supply pressure h0 (cm) could be approximated as
 | [8] |
where r0 is the radius of the infiltrometer (cm), K(h) is the unsaturated hydraulic conductivity, and hi is the pressure head corresponding to the initial water content. Typically Gardners' conductivity relationship (Gardner, 1958) is substituted into Eq. [8] to obtain a closed-form solution to the integral with the assumption that K(hi) is negligible. This permits the saturated conductivity and the exponent of Gardners' conductivity function, g, to be estimated by piecewise interpolation (Ankeny et al., 1991) or by least squares nonlinear regression (Logsdon and Jaynes, 1993; Hussen and Warrick, 1993). The K(h) relationship is not usually loglinear at pressure heads near saturation (e.g., h < 20 cm) and this can cause difficulties in the determination of K(h0) by the regression method (Logsdon and Jaynes, 1993). Although the piecewise estimation method partially removes the dependency of an assumed loglinear relationship, conductivities at the lowest and highest pressure heads are poorly estimated because g is extrapolated.
Parameter Optimization
The parameters of the constitutive relationships in Eq. [5], [6], and [7] can be estimated by minimization of the objective function (imnek and van Genuchten, 1996)
 | [9] |
where ß is the vector of optimized parameters, 
represents m vectors of independent variables, Nj is the length of the jth vector, 
j is the standard deviation associated with measurement errors, wi,j is the weight, and y is the measured response at each observation point i,j, and f (ß, i,j) is the predicted response as evaluated using Eq. [1] with the appropriate initial and boundary conditions.
Using numerically generated data simulating infiltration from a disc source, imnek and van Genuchten (1996) established that optimization of the cumulative infiltration alone results in a relatively intractable problem of nonunique parameter estimation. They concluded that other auxiliary information such as pressure head or water content measurements are required to improve parameter identifiability and convergence properties. Moreover, when using field data, parameter identifiability problems could be exacerbated further because of inherent errors in measuring cumulative infiltration and water contents in addition to errors caused by deviations of the flow from the invoked theoretical model (e.g., nonisotropic flow, nonstationarity of hydraulic properties with depth and time, temperature induced variations, air entrapment, etc.).
In this paper, we develop and evaluate several inverse optimization strategies and associated field methods for use with tension infiltrometers to estimate the hydraulic parameters of a fine-textured soil. Specifically, we compare parameter identifiability and the resultant fit to measured data among optimizations that, in addition to cumulative infiltration I(t), include as vectors in of the objective function (i) laboratory water retention data, LAB(h), from undisturbed soil cores; (ii) volumetric water contents, SC(z, T), from cores extracted at two depths, z, below the disc after termination of infiltration at time T; and (iii) transient measurements of volumetric water contents, TDR(t), using TDR probes inserted at the soil surface. In addition, we develop inverse methodology to estimate near-saturated hydraulic conductivity using multiple tension infiltration experiments.
|
MATERIALS AND METHODS
|
|---|
Field and Laboratory Experiments
Infiltration experiments were carried out on a fallowed no-tillage field and a native pasture (fine, mixed, superactive, thermic Torrertic Paleustoll) at the USDA-ARS Conservation and Production Research Laboratory, Bushland, TX. In the no-tillage field, infiltration experiments were completed at the surface and at a 20-cm depth. Single tension infiltration measurements were completed at a nominal potential of -15 cm H2O using a 0.2-m diam. disc infiltrometer. This infiltrometer permits infiltration to continue undisturbed while water is being resupplied (Evett et al., 1999). For multiple tension infiltration experiments, cumulative outflow was measured over a range of pressure heads, nominally -15, -10, -5, and -0.5 cm H2O. All measurements on cropland were made in nonwheel-tracked interrows. Infiltration plots were prepared by removing all vegetation and residues that would interfere with achieving a level surface. A layer of fine sand approximately 7- to 10-mm thick was placed over the surface to fill small depressions and facilitate contact between the soil and the nylon membrane of the infiltrometer. For some of the infiltration experiments, six three-rod, 20-cm TDR probes (Dynamax, Inc., Houston, TX, model TR-100)1 were inserted into the soil surface at a distance of 5 to 7 cm from the perimeter of the tension disc. Three of the probes were inserted vertically, and the remaining were inserted into the soil at a 45° angle downward from horizontal towards the disc center. Deionized water was permitted to infiltrate at each tension for at least 0.5 h. Water level in the infiltrometer tube was monitored with a pressure transducer at no more than 7.5 s intervals. Water contents were measured every 300 s using a TDR cable tester (Tektronix, Inc., Beaverton, OR, model 1502C)1 connected to the TDR probes through a coaxial multiplexer (Dynamax, Inc., Houston, TX, model TR-2001; Evett, 1998), both of which were controlled by a laptop computer running the TACQ program (Evett, 2000a,b).
Three sets of undisturbed soil samples (3 cm length by 5.4 cm diam.) were extracted 0.5 to 0.75 m from the disk center at depths of 1 to 4, 5 to 8, 11 to 14, and 15 to 18 cm to estimate initial water content. Two additional cores were extracted below the center of the disc at depths of 1 to 4 and 5 to 8 cm upon the termination of single-tension experiments to estimate the final water content. In the laboratory, after permitting these undisturbed soil cores to come to saturation, water retention curves were obtained using a hanging water column (0.2–15 kPa) and pressure plate apparatus (30–100 kPa). Equation [5] was fitted to retention data to estimate n, , and s using an adaptive, model-trust region method of nonlinear, least-squares parameter optimization (Dennis et al., 1981; Dennis and Schnabel, 1983). For these fits, the value of r was set to 0.005 cm3 cm-3 because it otherwise tended to take on values larger than water contents measured in the field.
Wooding's Solution
Many of the earlier described difficulties with the analysis of steady state infiltration using Wooding's solution can be overcome by substituting a K(h) function into Eq. [8] that is more flexible than Gardners' relationship. The VGM conductivity relationship is one such function that does not a priori assume log-linearity near saturation. We substituted Eq. [6] into Eq. [8] giving an expression with three unknowns (n, , and Ks). Steady state volumetric fluxes at each of the four supply pressures, Q(h0), were calculated using the final 300 s of outflow data. Parameters n, , and Ks were estimated by fitting Eq. [8] to the four steady state volumetric fluxes. The integral in Eq. [8] was numerically integrated using guassian quadrature. The fitted water retention characteristic curve in conjunction with measured initial water contents were used to estimate the initial pressure head hi. Hydraulic conductivities at the four supply pressures were obtained by substituting the optimized values of n, , and Ks into Eq. [6]. It should be noted that the fitted parameter values embody little physical meaning and only serve to calculate the hydraulic conductivity for a particular inlet pressure within the applicable experimental range.
Numerical Solution of Richards' Equation
A second-order, finite difference numerical method of lines procedure similar to that of Tocci et al. (1997) was used to solve the pressure head based form of Richards' equation in two-dimensions. The set of ordinary differential equations resulting from the spatial discretization of Eq. [1] was integrated over time using DASPK, a variable step-size, variable order integrator for differential algebraic systems of equations (Brown et al., 1994). A generalized minimum residual (GMRES) method (Saad and Schultz, 1986) was used to solve the nonlinear system at every time step. Backward differentiation formulas of orders one through five are used by DASPK to advance the solution in time and a local error control strategy is used to select the step-size and order of the integration. Although the mass-conserving mixed form of Richards' equation (e.g., Celia et al., 1990) is used in most codes to ensure mass balance, these algorithms integrate in time using low-order methods. Recent work using higher order time-stepping methods, however, has demonstrated that the pressure-head form can be accurate, economical, and numerically stable in the presence of sharp wetting fronts (Tocci et al., 1997).
A relative and absolute error tolerance of 1 x 10-3 cm for local error control and an initial time step of 3.6 x 10-9 s were used to obtain all numerical solutions. The finite difference grids for the method of lines solution were selected to ensure that mass balance errors within the solution domain remained <0.5% at all observed times throughout each infiltration experiment. The lower and right boundaries were normally set at 30 or 40 cm and the number of nodes along each axis ranged from 60 to 80. The initial soil water content profile was approximated in the model by a third-order b-spline interpolant so that average water contents integrated over depth corresponded closely (±0.001 m3 m-3) with water contents obtained from extracted soil cores. At depths >20 cm, water content was assumed to be constant to satisfy the lower boundary condition of free drainage. Although true initial water contents below 20 cm may not have been reflected by this assumption, this did not influence simulated and measured infiltration since wetting fronts did not extend beyond 15 cm from the surface. For experiments that employed TDR probes, interpolated initial water contents integrated over depth agreed closely (±0.02 m3 m-3) with average initial water contents measured using the TDR probes.
Predicted cumulative infiltration depth I(t) (cm) was calculated as
 | [10] |
where qz(r, 
) is the Darcy vertical flux density at t = . Nodal fluxes across the inlet surface boundary were integrated over time and radial distance using the trapezoidal rule. Mass balance error was calculated by summing all integrated boundary fluxes, dividing this value by the change in water volume, and subtracting this quotient from unity.
Optimization Strategies
Table 1 summarizes the field sites, experiment type, and the corresponding objective function that was minimized to fit hydraulic parameters using observed data. The residual water content r was set to a constant value of 0.005 as noted earlier. The saturated water content obtained from the fit of Eq. [5] to laboratory retention data or estimated from bulk density measurements was fixed in all subsequent inverse fits of cumulative infiltration to improve the identifiability of the remaining parameters.
View this table:
[in this window]
[in a new window]
|
Table 1. Summary of disc-infiltrometer experiments and components of the objective functions for each of the experimental plots.
|
|
Optimizations Using Multiple Tension Infiltration Data and Water Retention Measurements
Water retention data is often collected in conjunction with tension infiltrometer data. Yet in the majority of analyses they are treated independently of one another. The retention curve furnishes static information about the soil matrix whereas infiltration measurements contain dynamic information related to the capillary drive (Morel-Seytoux, 2001). Inclusion of LAB(h) along with multiple tension cumulative infiltration measurements in the objective function offers a means of incorporating both sources of information in the optimized parameters.
Cumulative infiltration is typically curvilinear at early times for the first of a series of ascending pressure heads, especially under dry soil conditions. This results from the absorption of water by the soil matrix and eventual filling of the available pore space. If near steady-state has been attained at the initial supply pressure head (e.g., -15 cm), then cumulative infiltration at subsequent higher pressure heads is nearly linear with time since the increase in available pore space within the wetted perimeter is typically negligible. These observations suggest that the shape parameters n and would be more sensitive to the cumulative infiltration curve at early times for the initial and lowest supply pressure head. In addition, K(h) calculated using the optimized value of Ks obtained from this infiltration stage would be most representative of unsaturated conductivities at potentials less than the supply pressure head. In contrast, only a small portion of the K(h) function near saturation would be sensitive to the cumulative infiltration at late times.
Accordingly, we used a stepwise strategy for inverse fitting of water retention and conductivity parameters using Richards' equation to estimate cumulative infiltration over multiple tensions. Initially n, , and Ks are optimized using both the laboratory water retention data and cumulative infiltration at the lowest imposed pressure head hp0 with the VGM functions to describe the constitutive relationships. At each of the succeeding and incrementally higher pressure heads (i.e., hp1, hp2, and hp3) Richards' equation was used to inversely fit these respective segments of the cumulative infiltration curve and successively estimate the piecewise conductivities K(hp1), K(hp2), and K(hp3). When the wetting front is contained within a homogeneous soil layer, the maximum principle guarantees that the maximum pressure head within the solution domain will be achieved at the surface boundary for this particular infiltration problem (Celia et al., 1990). This signifies that fitting K(h) in a piecewise manner can be achieved without extrapolating K(h) beyond the imposed inlet pressure head h0. For each consecutive fit, the (h) and K(h) relationships in Eq. [5] and [6], respectively, had already been optimized and are used to describe unsaturated flow at pressure heads less than the initial supply pressure head hp0. Again, we emphasize that the optimized value of Ks permits the description of unsaturated conductivity at pressure heads less than hp0 and may not be reflective of the true saturated hydraulic conductivity of the soil. This methodology assumes that soils are homogeneous to the depth of the wetting front penetration (typically <15 cm). Some hysteretic behavior is accommodated using this method since the fit of the (h) and K(h) relationships at h 
hp0 is carried out using cumulative infiltration that may reflect a certain degree of hysteresis that consequently would be exhibited in the optimized parameter values of and n. Completion of infiltration measurements with the disc infiltrometer at multiple tensions as described above yielded an optimized K(h) function that was compared with independently calculated conductivities using Wooding's Eq. [8].
Optimizations Using Final Volumetric Water Contents of Soil Cores
In situ real-time measurement of water contents during infiltration in the field is difficult without some soil disturbance adjacent to or beneath the disc infiltrometer. Typically, a thin layer of soil immediately beneath the infiltrometer disc is collected after the termination of the infiltration experiment (t = T) to determine the water content and facilitate the analyses aimed at estimating sorptivity or conductivity (Clothier and White, 1981; Smettem et al., 1994; imnek and van Genuchten, 1997). The volumetric water content estimated from these surface samples is considered to be in equilibrium with the inlet pressure head. Errors can result because of the small sampling depth required and the fact that the bulk density must be estimated from other measurements (Angulo-Jaramillo et al., 2000). We sampled water contents by extracting a 10-cm length soil core under the disc after the termination of the infiltration experiment. Cores were taken at the radial origin, where changes in water content with depth and horizontal distance are smallest, so that water content errors associated with positioning of the coring device would be minimized. In addition, the time at which cores were sampled, T, was recorded to permit the calculation of water content changes because of drainage after termination of infiltration. Once extracted, the cores were dissected to procure the 1- to 4- and 6- to 9-cm increments for water content determination and water retention measurements. Soil core water contents at the two depth increments SC(z, T) were included in the objective function by imposing a zero-flux surface boundary upon termination of infiltration and integrating (r, z) over space to numerically calculate average cylindrical water contents at t = T. This methodology has the advantage of measuring water content using a known soil volume and provides for a better description of water contents within the wetted soil volume. Because of the difficulties of extracting core samples from saturated or nearly saturated soils, we used this method only for infiltrometer experiments carried out at supply pressure heads less than about -15 cm.
Optimizations Using TDR Water Contents
In controlled laboratory settings fast response tensiometers and TDR have often been used to supplement cumulative outflow data with pressure head and water content measurements. Under field conditions, such auxiliary data is difficult to obtain simultaneously with outflow data without some soil disturbance caused by the installation of sensors. Nonetheless this information may greatly improve the identifiability of fitted hydraulic parameters. We used TDR to measure water contents over time TDR(t) below the infiltrometer disc. To minimize soil disturbance while maximizing the contact of the probe with the wetted soil volume, we inserted three TDR probes diagonally into the soil a few centimeters from the disc edge and oriented towards the origin. Topp et al. (1980) demonstrated that water contents measured by the TDR technique in the presence of wetting fronts are essentially equivalent to average water contents within the measurement volume. Water contents obtained using the TDR probes oriented at a 45° angle downward from horizontal towards the disc center TDR(t) were included in the objective function by integrating (r, z) over space to numerically calculate average water contents that would be detected by a TDR probe.
Minimization of the Objective Function
Minimization of the objective function was implemented using an adaptive, model-trust region method of nonlinear, least-squares parameter optimization (Dennis et al., 1981; Dennis and Schnabel, 1983). Derivatives with respect to each fitting parameter were calculated using forward differencing. Iterations of the nonlinear least-squares estimation procedure were continued until both the maximum scaled relative change in the parameters and the ratio of forecasted change in the residual sum of squares were <1 x 10-3. Combinations of three or fewer parameters were fitted to cumulative outflow data, water retention measurements, TDR water contents, and soil core water contents to identify the parameters sets that yielded convergence and the lowest sum of squared residuals (SSR).
Initial parameter estimates were selected based on fitted values from the water retention data and final steady state infiltration rates. Additionally, optimizations were always restarted using different initial guesses to evaluate the possibility that previous fits converged on local minima and to ascertain if optimized parameters converged to similar values. The initial pressure head distribution with depth required for the solution of Richards' equation was always calculated from the measured initial water content distribution at each iteration of the nonlinear, least-square solver. Using numerically generated data, imnek and van Genuchten (1997) found that expressing the initial condition in terms of water content led to a better identifiability of parameters as compared with using initial pressure head data.
Predicted and measured cumulative infiltration I(t) in the objective function were expressed as increments of volume per unit area of the infiltrometer base (cm3 cm-2). Cumulative infiltration recorded in the field typically consisted of over 1000 data points. To reduce the storage requirements in the nonlinear fitting routine, data included in the objective function were limited to those collected at 180-s time increments and at times when TDR water contents were measured. Cumulative infiltration at these specified times were calculated by using a 9-point centered linear fit of measured outflow volumes. The three-rod TDR probes were assumed to measure the average water content within a 3 by 8 by 20 cm3 right rectangular prism. Predicted water contents for soil cores and TDR probes were calculated by integrating over their respective soil volumes using the trapezoidal rule. These values were fitted directly to water contents obtained from soil cores or average water contents measured by the three TDR probes oriented at 45° for the recorded times in the field.
Squared residuals in Eq. [9] were weighted equally (wi,j = 1) for all data sets. Residuals for the data sets I(t), LAB(h), SC(z, T), and TDR(t) were normalized using the number of observations in each data set and a measurement variance of unity in Eq. [9]. A variance of unity was used since measured and estimated standard deviations associated with each of the data sets were similar in magnitude. For instance, average standard deviations for a surface no-till plot were 0.02 cm3 cm-3 (N = 3) for TDR(t) at 45°, 0.03 cm (N = 24) for I(t) based on pressure transducer signal variations, and 0.03 cm3 cm-3 (N = 12) for LAB(h) at each pressure head. Based on initial water content measurements, estimated standard deviations for volumetric water contents from extracted soil cores ranged from 0.01 to 0.03 cm3 cm-3 (N = 3). For this narrow range in standard deviations among all data types, limited justification exists for assuming other than equal error variances.
The time and depth coordinates of measured infiltration data were transformed to account for the presence of a layer of contact sand. As per Vandervaere et al. (2000), the infiltration depth associated with the sand layer, Is, was calculated as 
x sand
where Vs is the measured volume of sand and sand(h0) is the associated available pore space. The time period for measured cumulative infiltration to achieve Is was subsequently defined as ts. Accordingly, cumulative infiltration measured in the field was transformed to (I - Is) and time was transformed to (t - ts) for inclusion in the objective function. The contact sand used in this study (Ottawa F-110)1 had a water content of 0.38 ± 0.02 for the range in supply pressure heads of -16 to 0 cm H2O.
|
RESULTS AND DISCUSSION
|
|---|
Water Retention Fits
Results of the fit of the parameters in Eq. [5] to retention data are shown in Table 2 for each of the sites. The fitted values of , n, and s for each site and depth represent the aggregate fit of all soil cores obtained from the particular site since fitted parameters did not significantly differ among plots. In fits to retention data obtained from both the surface and subsurface of the no-tillage plots, s was set equivalent to the porosity estimated from mean bulk density because fitted values were overestimated in these plots. The estimated retention parameters, which represent drying curves, were used as initial or fixed values in subsequent optimizations using cumulative infiltration data. The use of drying curves permits the establishment of an upper limit on retention curves fitted using infiltration data, which contains information about wetting. Therefore, retention curves fitted using infiltration data (represented by a wetting branch) should be equivalent to the laboratory curve or offset toward lower water contents when hysteresis is manifested.
View this table:
[in this window]
[in a new window]
|
Table 2. Results of inverse optimizations obtained by including both cumulative infiltration I(t) and laboratory retention measurements LAB(h) in the objective function. (Values in italic signify fitted parameters and values in parenthesis below signify the 95% confidence interval as calculated from asymptotic standards errors.)
|
|
Optimization of K(h) and (h) Using Infiltration and Water Retention Measurements
Results of inverse optimizations for parameters when fitted to both laboratory retention measurements and cumulative infiltration at a single supply pressure head are summarized in Table 2. For all the optimizations in which n, , and Ks were simultaneously fitted, the minimization algorithm had difficulty converging. Moreover, Ks converged to substantially dissimilar values (>50%) when the fitting procedure was restarted with different initial estimates. Inspection of the response surface of the objective function within the –Ks parameter plane (Fig. 1)
demonstrates why convergence problems were problematic for these fits. The long narrow valley exhibited by the response surface suggests that Ks is not uniquely defined (computationally at least) for these optimizations when fitted simultaneously with . Toorman et al. (1992) also demonstrated that there was a positive correlation between and Ks for optimizations using numerically generated one-step outflow data. They attributed this identifiability problem to the small importance of gravity for the short cores used in the study. This is made evident by and Ks appearing only as a ratio of each other when the VGM relationships are substituted into Eq. [1] and the gravity term is dropped. The geometry-induced enhancement of capillarity over gravity for unconfined infiltration in fine-textured soils suggests that the simultaneous identification of and Ks would be exacerbated for our optimization problems. To avoid this identifiability problem, we fixed to the value obtained by the fit of Eq. [5] to laboratory water retention data. Inverse optimizations with fixed resulted in small 95% confidence intervals for the estimates of n and Ks (Table 2) and a well-defined minimum as demonstrated by the response surface within the n–Ks parameter plane (Fig. 2)
. Moreover, the minimum converged value of the objective function for these two-parameter fits was not more than 6% greater than the minimum achieved with the three-parameter fits. The fitted values of n show remarkable consistency for plots within the same field and fall within the expected range for clay and silty clay soils (Yates et al., 1992).
Measured and optimized cumulative infiltration depths for each of the sites and the corresponding water retention functions are plotted in Fig. 3 and 4
. Fitted cumulative infiltration depths corresponded closely to measured depths except at very early times, especially for the no-tillage plots (Fig. 3 and 4). For the no-tillage at 20-cm depth and the native pasture plots, cumulative infiltration was satisfactorily fitted with only slight modifications to the fitted water retention parameters n and . In contrast, n converged to values larger than that obtained from fits to retention data from the no-tillage surface plots. This resulted in a water retention curve that was displaced from the drying retention data towards lower water contents (Fig. 4). We believe that hysteresis of the soil hydraulic functions contributed to the disparity between the retention curve predicted from infiltration and those measured in the laboratory. Hysteresis was not apparent in the other plots (Fig. 3) possibly as a result of relatively larger initial water contents (i > 0.23, see Table 1) that narrowed the range in water contents inside the wetted region during the experiment. For the no-tillage surface plots, however, water contents ranged from 0.05 to 0.40 during the course of the infiltration experiments.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 3. Measured cumulative infiltration and water retention data and corresponding two parameter optimized curves for the subsoil no-tillage and native pasture plots (Table 2). Error bars represent 95% confidence limits for the mean water content from soil cores sampled at the 1- to 4- and 6- to 9-cm depth increments in the native pasture plots and 11- to 14- and 15- to 19-cm in the no-tillage plots. The dotted line in no-tillage graphs represents the fitted curves obtained by fitting only the VGM constitutive relationships to the entire infiltration curve at all supply pressure heads (see section, Optimization of K(h) Near Saturation Using Multiple Tension Infiltration Data).
|
|
View larger version (29K):
[in this window]
[in a new window]
|
Fig. 4. Measured cumulative infiltraiton and water retention data and corresponding two-parameter optimized curves for the no-tillage surface plots (Table 2). Error bars represent 95% confidence limits for the mean water content from soil cores sampled at the 1- to 4- and 6- to 9-cm depth increments in both plots.
|
|
Optimization of K(h) Near Saturation Using Multiple Tension Infiltration Data
The optimized parameters obtained from the fit to both the cumulative infiltration and laboratory retention measurements (Table 2) were next used in the constitutive Eq. [5] and [6] as constants to solve Richards' equation and fit the cumulative infiltration measurements obtained at supply pressure heads greater than hp0. This permitted the sequential, one-parameter fits of K(hp1), K(hp2), and K(hp3) (Table 3). For these fits, only the cumulative infiltration data falling within the supply pressure head (hp1, hp2, or hp3) were weighted to unity in the objective function. The weights of all other cumulative infiltration data (at earlier times) were set to zero. Optimizations carried out in this manner yielded estimates of conductivity at each supply pressure head with relatively narrow range (less than ± 8% normalized) in the 95% confidence limits (Table 3). Excellent agreement between measured and optimized cumulative infiltration depths were obtained using the piecewise method (Fig. 5)
that could not otherwise be achieved using the VGM model to define K(h) throughout the entire range in pressure heads. For instance, a single inverse fit of the cumulative infiltration and water retention data for Plot 1 of the no-tillage subsoil using only the VGM relationship yield SSRs for infiltration data and retention data six and 50 times greater, respectively, than those obtained using the piecewise method (see Fig. 3 and 5). imnek et al. (1998b) also used the VGM equations over the entire range in pressure heads to inverse fit hydraulic parameters to multiple tension infiltrometer data. They obtained a good fit to cumulative infiltration but the predicted water retention curve seriously underestimated retention data obtained from laboratory measurements. Our results suggest that if a single fit of the VGM model is used over the entire pressure range then large values of Ks and smaller values of n are required to adequately describe the conductivity and water retention relationships of these fine-textured soils at high potentials, which in turn poorly represents K(h) near saturation.
View larger version (16K):
[in this window]
[in a new window]
|
Fig. 5. Measured cumulative infiltration data and corresponding fitted curves for the subsoil no-tillage and native pasture. Inverse fits to cumulative infiltration after the first supply pressure head were obtained using a loglinear piecewise description of the K(h) function (Table 3). Inverse fits for the first supply pressure head are shown in Table 2. Step changes in the supply pressure head are indicated by the symbol  in the cumulative infiltration plots. The dotted line in the no-tillage graph represents estimated cumulative infiltration obtained by fitting only the VGM constitutive relationships to the entire infiltration curve at all supply pressure heads.
|
|
The K(h) functions derived from the four optimizations at the four supply pressure heads are presented for each plot in Fig. 6 and compared with values calculated for each supply head using Eq. [8]. Although we estimated the integral in Wooding's Eq. [8] using the initial pressure head hi calculated from the initial water content and the water retention curve, setting hi to a large negative value (-1000 cm H2O) produced essentially identical results. The conductivity at each supply pressure head calculated using Eq. [8] compared closely with the optimized K(h) function for all plots (Fig. 6). Such close agreement between estimated conductivities implies that Wooding's analysis of steady state infiltration rates is valid even for the silty clay soil used in this study which, based on numerical studies (Warrick, 1992), should approach steady state flow conditions at times far in excess of the approximately 1 to 1.5 h we used in this study. imnek et al. (1998b) also demonstrated good correspondence between K(h) obtained by inverse optimization and Wooding's analysis except at higher pressure heads. We obtained better agreement between these two analyses, especially at the highest pressure head, probably because a piecewise description of K(h) was used near saturation.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 6. Unsaturated hydraulic conductivities at each supply pressure head calculated using Wooding's analysis (symbol) and the corresponding optimized hydraulic conductivity function (line) obtained from four sequential inverse fits (Tables 2 and 3).
|
|
Optimizations Using Final Volumetric Water Contents of Soil Cores
Results of parameter optimizations that included in the objective function both infiltration data and the volumetric water contents from cores extracted after termination of each experiment are shown in Table 4. As with the previous optimization results, three-parameter fits of , n, and Ks, led to convergence problems because of nonuniqueness in the –Ks parameter plane (Fig. 7) . To investigate this problem further, we numerically generated infiltration and water content data and subsequently used these data for inverse optimizations. Doing so led to excellent convergence properties for these three-parameter fits, similar to results obtained by imnek and van Genuchten (1997). We also completed inverse optimizations using the generated final volumetric water content data with added or subtracted deterministic errors (±0.02) (e.g., imnek and van Genuchten, 1997). These optimizations also converged to parameter estimates close to the true values, however asymptotic standard errors of the estimates were significantly larger than error-free data.
View this table:
[in this window]
[in a new window]
|
Table 4. Results of inverse optimizations obtained by including both cumulative infiltration I(t) and volumetric water contents of extracted soil cores SC(z, T) in the objective function. (Values in italic signify fitted parameters and values in parenthesis below signify the 95% confidence interval as calculated from asymptotic standards errors.)
|
|
Experimentation with changing the standard deviation for the water contents in the objective function indicate that all three parameters become identifiable using field measured data if is decreased 10-fold. However, these results indicate that to attain identifiability, SC(z, T) would need to be measured with a standard error of 0.003 cm3 cm-3, a level of accuracy that is, in practice not attainable considering that only one sample for each depth increment can be extracted after termination of the infiltration experiment. We speculate that deviation of the infiltration process from the invoked theoretical model may be influencing the optimization results as do measurement errors in the data.
All two-parameter fits of n and Ks converged to estimates with relatively small 95% confidence intervals (Table 4), and with good agreement between measured and optimized cumulative infiltration depths. Fitted volumetric water contents underestimated measured water contents (Table 4), especially for the 6- to 9-cm depth increment. However, all but one of the estimated water contents had acceptably small residuals within the expected range of sampling error of about ±0.03 cm3 cm-3. Simulated drainage after termination of the infiltration experiment and before soil core extractions (about 1–2 min) indicated only a minor reduction (<0.017) in the volumetric water contents of extracted soil cores.
Optimizations using SC(z, T) data resulted in significantly larger parameter estimates of n and smaller errors in the fitted cumulative infiltration for Plots 1 and 2 of the no-tillage field (Table 4) as compared with the optimizations that used the water retention data (Table 2). Optimizations using SC(z, T) data led to a lowering of the fitted water retention curve below that of the laboratory retention data (Fig. 8)
, likely because of hysteresis. Hysteresis was manifested by an increase in the fitted value of n and was more strongly expressed in infiltration experiments with lower initial water contents. Since initial water contents varied with depth, the fitted water retention function represents a lumped scanning curve rather than any single scanning curve. The identification of parameter estimates that could describe hysteretic relationships will require optimizations using modified retention and conductivity functions that account for hysteresis such as imnek et al. (1999a).
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 8. Measured water retention data and corresponding optimized retention curves obtained from the inverse fit to infiltration data and final volumetric water contents from extracted soil cores (Table 4). Error bars represent 95% confidence limits for the mean water content from soil cores sampled at the 1- to 4- and 6- to 9-cm depth increments in all plots within a given field.
|
|
Optimizations Using TDR Water Contents
Results of inverse optimizations that included TDR water contents as well as cumulative infiltration in the objective function are summarized in Table 5. The lower portion of the Ap horizon at depths greater than about 10 cm in no-tillage surface plots possessed hydraulic properties that differed from the surface layer and were more representative of retention characteristics for the Bt horizon at 20 to 30 cm. The second layer did not significantly influence cumulative infiltration because the wetting front was contained principally in the upper 10 cm. But, the TDR probes inserted at 45° extend into this second layer; and TDR measurements of water content did reflect the mean water content of both layers. Simulated drainage using the best fit parameters and a single layer caused predicted water contents to decrease within soil volumes measured by TDR probes at early times. The simulated water content decrease did not agree with TDR data that indicated stable water contents prior to infiltration. To address this difficulty, we simulated infiltration in no-tillage surface plots using two soil layers at 0 to 10 and 10 to 40 cm. Parameters from the two-parameter fits (Table 2, no-tillage at the 20-cm depth, Plot 1) were used to simulate water flow in the lower layer; and the parameters for the 0- to 10-cm layer were obtained by inverse parameter estimation. We emphasize that the hydraulic properties of the second soil layer had an insignificant influence on cumulative infiltration and fitted parameters. For example, a ten-fold decrease in the saturated conductivity in second layer yielded inverse fitted parameters that varied only 3 to 7% from the estimates using the unmodified hydraulic properties of the second layer.
View this table:
[in this window]
[in a new window]
|
Table 5. Results of inverse optimizations obtained by including in the objective function both cumulative infiltration I(t) and water contents measured by time domain reflectometry (TDR) probes TDR(t) oriented 45° from perpendicular. (Values in italic signify fitted parameters and values in parenthesis below signify the 95% confidence interval as calculated from asymptotic standards errors.)
|
|
For the optimizations in which n, , and Ks were simultaneously fitted, the minimization algorithm had difficulty converging. As with previous optimizations, these three-parameter fits led to nonunique solutions as indicated by influence plots in the –Ks parameter plane. Two-parameter fits of n and Ks using TDR data (Table 5) converged to estimates with values similar to those obtained for optimizations using final water contents from soil cores (Table 4). Also, final simulated water contents deviated from TDR data by 0.01 to 0.03 cm3 cm-3, again similar to the results for optimizations that used final water contents from soil cores. However, optimized TDR water contents were significantly underestimated at early times (Fig. 9)
. Correspondence between simulated and measured TDR water contents was especially poor for Plot 2 of the no-tillage site despite the fact that TDR data comprised 20% of the error in the objective function. Optimizations using a ten-fold increase in the weight of TDR data in the objective function led to significantly larger parameter estimates of n and Ks but failed to yield any significant improvement in the simulated water contents (Fig 9). Simulations predicted sharp wetting fronts. In contrast, the measured data showed a much earlier arrival of the wetting front and a gradual increase in water contents thereafter. Differences between simulated and measured water contents at early times were probably due to physical nonequilibrium. At later times in the infiltrometer experiments, water contents probably began to attain near equilibrium conditions that led to a better agreement between measured and simulated water contents. The greater weighting of TDR water contents, and hence the greater emphasis of nonequilibrium conditions at early times, resulted in fitted parameters that were more representative of coarser-textured soils.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 9. Mean measured (symbol) and fitted (line) water contents for TDR probes for surface no-tillage plots. Fitted water contents were obtained by inverse optimization of infiltration data and 45° TDR water contents. Optimized parameters were then used to predict water contents measured by 90° TDR probes. Error bars represent the average 95% confidence interval about the mean water contents measured by TDR.
|
|
|
SUMMARY AND CONCLUSIONS
|
|---|
The insensitivity of the objective function over a wide range in Ks for the three-parameter fits make simultaneous identification of Ks, , and n very difficult, if not impractical, for all optimizations investigated. We attribute a portion of this identifiability problem to the enhancement of capillarity over gravity for unconfined infiltration in fine-textured soils. Inclusion of measured soil core water contents at the termination of infiltration experiments did not improve the identifiability of Ks and for these three-parameter fits. We speculate that deviations of the flow from the invoked theoretical model are also influencing the optimization results as much as, or more than, unavoidable measurement errors in water content and cumulative infiltration. Based on these results, we recommend that be estimated using water retention data and thereafter be fixed at this value for inverse fits to cumulative infiltration data. For these soils, the two-parameter fits with held constant improved the identifiability of Ks and n while not compromising the fit to measured infiltration. We emphasize, however, that the optimization strategies developed for the fine-textured soils in this study may not necessarily be appropriate for coarser-textured soils.
For two-parameter fits, minimizations of the objective function that included both cumulative infiltration and drying water retention data led to excellent fits for those experiments that had relatively high initial water contents (i > 0.23 m3 m-3). At lower initial water contents, good fits to cumulative infiltration were obtained only with an estimated retention curve that exhibited hysteresis as compared with measured water retention data. In those cases where initial water contents were low, even better fits to cumulative infiltration could be obtained by minimizing the objective function that included both cumulative infiltration data and volumetric soil water content measured upon the termination of the outflow from the disc infiltrometer. These optimizations also led to estimates of final water contents that closely approximated measured water contents (within 0.03 cm3 cm-3) irrespective of whether water contents were measured using soil cores or TDR.
Optimizations that included both cumulative infiltration and volumetric water contents measured by TDR with probes oriented 45° from horizontal returned parameter estimates that were similar to values obtained for fits that used final soil water contents. These optimizations resulted in final simulated water contents that were within 0.03 cm3 cm-3 of those measured by TDR. However, transient flow near the margins of the wetting front was not well described by Richards' equation probably because of physical nonequilibrium processes. As a consequence, water contents measured by TDR in these regions were poorly estimated early in the simulations; and the identifiability of hydraulic parameters for these minimization problems was not improved over that obtained using soil cores for final water contents. Final water contents measured by TDR and from soil cores extracted under the tension infiltrometer were fairly well predicted, likely because of equilibrium having been achieved.
Inverse optimizations over multiple supply pressure heads, with K(h) defined using piecewise loglinear interpolation at pressure heads near saturation and the VGM function at lower pressure heads, resulted in excellent fits to measured cumulative infiltration. In contrast, inverse optimizations using theVGM function to describe K(h) at all pressure heads resulted in much poorer fits to infiltration data and unsatisfactory predictions of the water retention characteristic curve. Wooding's analysis of multiple tension infiltration experiments yielded estimates of K(h) near saturation that compared closely with inverse optimization estimates.
Based on the results of this study, both single tension infiltration experiments with final soil water contents and multiple tension experiments would be suitable to characterize hydraulic properties of fine-textured soils. Water retention measurements from soil cores would be required to estimate . Optimizations of multiple-tension infiltration experiments would provide estimates of conductivity near saturation. In contrast, optimizations using cumulative infiltration data at approximately -15-cm potential, in conjunction with final volumetric water contents, would furnish more precise information pertaining to hydraulic conductivities at larger potentials. Ordinarily, only a single soil core can be removed at the termination of infiltration experiments. Obviously, this sample must be taken with care to ensure that errors in measured volumetric water contents are small. The use of TDR can reduce errors in water content measurements, however probes should be placed well within the final wetted perimeter to maximize the collection of relevant data for use in optimizations.
|
ACKNOWLEDGMENTS
|
|---|
We acknowledge the helpful discussions and exchanges with J. imnek.
|
NOTES
|
|---|
TOP
ABSTRACT
INTRODUCTION
[I(t), 
