Vadose Zone Journal 1:239-251 (2002)
© 2002 Soil Science Society of America
Continuous Flow Method for Rapid Measurement of Soil Hydraulic Properties
I. Experimental Considerations
G. L. Butters* and
P. Duchateau
Department of Soil and Crop Sciences and Department of Mathematics, Colorado State University, Fort Collins, CO
* Corresponding author (gbutters{at}agsci.colostate.edu)
Received 23 April 2002.
 |
ABSTRACT
|
|---|
Soil hydraulic properties are important in many vadose zone processes, but the measurement of these properties is usually tedious and often difficult. We describe a continuous flow method that allows very rapid and accurate measurement of hydraulic conductivity, K(h), and moisture retention, 
(h), functions including hysteresis. The experimental design described here employs simultaneous tensiometry and water flow measurements that are easily automated. The design provides a highly versatile approach to controlling draining and/or wetting rates. The analysis uses a combination of direct Darcian analysis and numerical inversion of Richards' equation for estimation of the hydraulic properties. This combination allows for estimation of wetting and/or draining K(h) and (h) over the entire tensiometer range of measurement, while retaining the physical significance of the hydraulic parameter estimates. Estimation of K(h) and (h) to lower soil water pressures can be accomplished by supplying the inverse analysis with an independently measured water content such as (1.5 MPa). The key to applying the inverse analysis to variably saturated conditions unambiguously is identification of the air-entry pressure of the draining soil. We show how the air-entry pressure can be identified from an inflection in the soil water pressure and changes in drainage flux. The air entry identification is verified by measurement of the soil air pressure. The continuous flow approach is applied successfully to a variety of soil textures and structures. The effects of drainage rate and range of measurement are also demonstrated and discussed.
|
INTRODUCTION
|
|---|
WATER FLOW AND RETENTION in the vadose zone is of fundamental importance to soil scientists and hydrologists and is a critical element in assessing the environmental implications of soil management. The two key hydraulic properties, the water content–pressure head relationship, (h), and the conductivity–pressure head relationship, K(h), have been the subject of a great deal of research. Direct (Klute, 1986; Klute and Dirksen, 1986) and indirect (van Genuchten and Leij, 1992) laboratory and field methods abound in the soils, hydrology, and engineering literature. It is generally conceded that measurement of K(h) and (h) is typically time consuming and usually difficult. Recently, advances in numerical techniques have increased the attractiveness of inverse methods (Hopmans and Simunek, 1999) wherein K(h) and (h) are not explicitly measured but instead are implicitly deduced from the soil water flow in response to a well-defined perturbation of the soil water energy. The most common approach is to employ Richards' equation as the governing flow equation in combination with water flux and/or soil water pressure measurements to identify parameters of (h) and K(h) by a least-squares numerical iteration. An attractive advantage of the inverse approach is the time savings of measuring both (h) and K(h) in a single experiment (assuming a linkage of these functions). Kool et al. (1985) and Parker et al. (1985) showed the promise of the technique using records of cumulative drainage following a single large step of air pressure in laboratory soil cores. Later, Eching and Hopmans (1993), Eching et al. (1994), and van Dam et al. (1994) showed that including both water flux data and soil water pressure data as well as conducting the experiment in a series of small pressure steps over a few to several days reduced the uniqueness problems in identifying (h) and K(h). Using a gradual change in the soil air pressure, a continuous flow approach for (h) was demonstrated by Salehzadeh and Demond (1994). Their automated apparatus demonstrated good versatility in measuring moisture retention, but was deemed unreliable for K(h) due to uncertainty in estimation of hydraulic gradients. Wildenschild et al. (1997) suggested a two-stage approach to estimating soil hydraulic properties: a slow, continuous flow measurement of (h) for about 3 d, followed by resaturation and a one-step outflow measurement for inverse estimation of K(h). As rationale for this two-stage approach, Wildenschild and coworkers argued that the flow rates used for the (h) measurement were too low for either direct or inverse estimation of the hydraulic conductivity. However, estimation of the hydraulic conductivity at higher continuous flow rates allowing a single flow experiment capturing both (h) and K(h) was not tested due to uncertainty in the lower boundary condition at the higher flow rates.
An unresolved limitation of the inverse approach is the lack of reliability near saturation so that parameters of the estimated K(h), in particular the saturated hydraulic conductivity (KS), are fitting values of uncertain physical relevance. Durner et al. (1999) performed multistep and gradual pressure change (drainage or wetting in the pressure range of 0 to -60 cm over about 100 h) and suggested the continuous pressure change as an attractive alternative to multistepping provided that either KS or S are measured independently. Consistent with earlier research, the authors expressed difficulties with the inverse approach near saturation and concluded that "... KS can never be accurately estimated from an inflow–outflow experiment because the conductivity function becomes very insensitive near saturation." Despite the limitations, inverse analysis for estimation of soil hydraulic properties is a powerful tool and is gaining popularity (Kodesova et al., 1998; Simunek et al., 1998a; Inoue et al., 1998) in favor of tedious direct methods of assessing unsaturated flow properties.
The purpose of this research was to develop and apply a measurement technique for (h) and K(h) or K() that is both rapid and accurate. Specifically, the method detailed here combines direct and inverse analysis of tensiometry and continuous inflow–outflow data to measure both draining and wetting (h) and K(h) over the tensiometer range of soil water pressures, including saturation. Since this continuous flow technique does not employ equilibrium steps, it is uncommonly rapid but without the sudden, large pressure changes used in single and multistep approaches. Complete characterization of (h) and K(h) over the tensiometer range of soil water pressures, including hysteresis, can be accomplished in less than a day or two. The focus of this paper is the description of the experimental design and discussion of the direct and inverse analysis with particular attention to the near saturation dynamics. Demonstrating the direct and inverse recovery of hysteretic (h) and K(h), including scanning curves, from continuous flow data is the subject of ongoing research.
|
MATERIALS AND METHODS
|
|---|
Experimental Design
The laboratory design for a continuous flow measurement (Fig. 1) consists of three basic elements: (i) a gas phase barrier (ceramic or membrane plate) in contact with a pressure controlled reservoir, (ii) tensiometry to monitor soil water pressures at the boundaries (or within) the soil, and (iii) an electronic balance to monitor water flow (indicated by changes in weight). For the base plate and tensiometers, we routinely use 0.05- or 0.1-MPa (0.5- or 1-bar) high-flow ceramics. All plate surface exterior to the soil are covered with a layer of silicone to prevent evaporation. Pressure measurements are most conveniently performed using pressure transducers (we use Validyne Mdl. DP15-42; Validyne Engineering, Northridge, CA). The automated data acquisition that includes the PC-linked electronic balance is not a necessity, but it is strongly recommended. The most interesting aspect of the experimental design is the air pressure manipulation of the reservoir head space. This is the key to initiating flow in the soil sample and provides a versatile yet simple approach to manipulating the soil water pressure. Changing the air pressure in the reservoir is analogous to changing the elevation of the free water surface and the rate of air pressure change controls the rate of draining or wetting of the soil. Note that the soil air pressure is not explicitly manipulated but may change naturally with drainage or inflow. By adjusting the speed of the air pump to the reservoir, the design allows a great deal of versatility in controlling the rate of draining or wetting and the selection of turning points for reversing the process. Since the soil water pressure is measured at both soil boundaries, it is not necessary to record the reservoir level or the air pressure in the reservoir head space. This experimental design is similar to that of Salehzadeh and Demond (1994), but without the restricted sample length or small (1 cm) tensiometer spacing, which causes problems with resolving hydraulic gradients with even a small degree of electronic noise in the tensiometry. Further, the use of the air pump here greatly simplifies air pressure control, while the electronic balance provides exceptional, noise-free resolution of the outflow or inflow.
In a typical experiment, we begin by selecting the soil water pressure (in the tensiometer range of operation) at which to initiate wetting or draining. For example, suppose the soil sample of length L is saturated with the desired flow solution and allowed to drain to equilibrium with the free solution surface at the elevation of the plate. In this case, the initial soil water pressure (in head units) is -L/2 at the mid point of the sample, and further drainage is initiated by lowering the air pressure in the reservoir. While decreasing the air pressure at a desired rate (e.g., 1 cm min-1), the lower and upper boundary soil water pressures are recorded (typically at a 1-min interval) as is the weight of the soil sample. We have verified in separate experiments that the response time of the tensiometers with the null-type pressure transducers is very rapid and, relative to the soil water pressure changes we employ, provides an essentially lag-free record of soil water pressure with time. The drainage is allowed to proceed to the desired pressure in the tensiometer range. Wetting the sample may then proceed by reversing the air pump to the reservoir and slowly increasing the air pressure. The process (draining or wetting) may be reversed at any point, allowing recovery of drainage and wetting hydraulic functions, as well as the scanning curves connecting these functions. As a standard flow solution, we use 0.005 M CaCl2, with thymol (0.05 g L-1) as a microbial inhibitor (Klute,1986) for prolonged measurements, though the properties of the flow solution could easily be adjusted according to the purposes of the investigation.
As a caution, notice that the flowcell rests on the balance but is not free of connecting tubing during the experiment. The tensiometer and outflow tubing must have a degree of flexibility but should not deform over the pressure range of the experiment. We have found that thick-walled Masterflex tygon tubing is satisfactory. After the tubing is connected, it is usually necessary to adjust the tubing into a relaxed position so as to avoid drift in the weight measurements.
Analysis
The experimental methodology as described above affords the opportunity to combine direct and inverse analysis. The direct analysis is particularly useful in providing near-saturation data as input to the inverse estimation of the hydraulic properties. The benefit of this approach is that the hydraulic properties near saturation, particularly S and KS, retain physical significance. Further, as we will show, the direct analysis aids in identifying the air-entry soil water pressure, which will allow us to apply the inverse analysis near saturation without ambiguity.
Direct Analysis
The approach used here is essentially that suggested by Ahuja and El-Swaify (1976). For (h), the direct analysis is a straightforward mass balance where 
is determined from the change in weight of the soil sample of known volume. Using a small soil core (e.g., 50 mm in diameter by 30 mm long), the 0.01-g weight sensitivity of the digital balance allows virtually noise-free resolution of to about 0.0002 mm3 mm-3. Provided a single is known (usually the initial ), (h) is constructed from the weight changes and the corresponding arithmetic average h of the sample. Near saturation, provided the sample is short and the change in boundary pressure is slow (see discussion below), only a small difference in h between the upper and lower boundaries of the sample is observed and the average h is fairly precise. At lower hydraulic conductivities, h over the sample length becomes larger and consequently precision in the h estimate is lost. The decision to invoke inverse analysis for (h) instead of direct analysis thus depends on the desired range of h and tolerance to estimation errors in the average h of the sample.
The direct analysis for hydraulic conductivity is accomplished using Darcy's Law with the observed h and sample weight data. For a time interval t (e.g., 1 min.), the average lower boundary water flux, q(0,t), is determined from the change in sample weight ( wt) as
 | [1] |
where 
fs is the density of the flow solution and A is the cross-sectional area of the sample. Given the 0.01-g resolution on wt, the minimum observable flow per unit area for a typical field core (e.g., diameter = 50 mm) is 0.0005 cm (assuming fs = 1 g cm-3). Assuming no evaporative water losses at the upper boundary, q(L,t) = 0, the water flux at the midpoint of the sample (L/2) is estimated as the average of the boundary fluxes, q(0,t)/2. The hydraulic gradient at any time is determined from the boundary tensiometry, H/z = [h(L,t) + L - h(0,t)]/L. Using the average H/z over the time interval, K = -q(L/2,t)/(H/z) and is associated with the average h or average of the time interval. We reemphasize that h is measured at the lower soil boundary (z = 0) so that the conductivity of the plate is irrelevant in the direct calculations and the inverse estimation. Naturally, a low conductivity plate reduces response time of the soil water pressure to changes in the reservoir pressure, but it is only necessary that h(0,t) is known. We have found that inferring the lower soil boundary pressure using the plate conductance as is sometimes suggested (Wildenschild et al., 1997) unnecessarily limits flow rates or is unreliable due to changing plate conductance. Since both the direct and inverse calculations are sensitive to uncertainty in the lower boundary pressure, we favor the trade-off of accuracy in h(0,t) over the disturbance penalty of installing a lower boundary tensiometer.
The accuracy of the direct K estimate is sensitive to the length (L) of the soil sample and to the rate of pressure change at the lower boundary (P) relative to the hydraulic conductivity. If P and L are both large in a soil of low K, the hydraulic gradient in the sample quickly becomes nonlinear, and H/z and q(L/2,t) are poorly approximated. Under these conditions in a draining soil, much of the h occurs over a fairly short distance near the lower boundary and the overall H/z interpreted from the boundaries is smaller than the actual gradient near the lower boundary. Further, the nonlinear hydraulic gradient in the draining soil near saturation is accompanied by a nonlinearly increasing flux (moving downward from z = L to z = 0) so that the magnitude of the actual midpoint flux (|q(L/2,t)|) is less than the estimate from the boundaries (|q(0,t)/2|). The compounding effect of the nonlinear h and q profiles causes overestimation of K by the direct calculation.
Guidance as to combinations of L and P allowing accurate recovery of K by the direct calculation is given by the contours in Fig. 2. The contours of the ratio (Kestimated/Kactual) were constructed using Hydrus-1D (Simunek et al., 1998b) to solve the forward problem for multiple combinations of L and P. In each case, the simulated soil was drained from equilibrium (h = 0 at z = 0) by imposing negative P at the lower boundary. The resulting q(0,t) and h(L,t) record was used to calculate K directly by the Darcian analysis discussed above. The ratio of the estimated K to the actual K at a given pressure head (in the range L/2 < |h| < 15 cm) is depicted in the figure for hypothetical soils of contrasting K and pore-size distributions. It is apparent in Fig. 2 that short samples (< 
5 cm) with slow pressure change (< 1 cm min-1) favor the accurate recovery of near saturation K. The fairly limited range of permissible (L, P) pairs for the coarse-textured soil (Fig. 2c) is a pore-size distribution effect. The very sharp decrease in K with decreasing h in the coarse soil leads to nonlinear hydraulic gradients even close to saturation so that small P is advised. In practice, however, if P is too small, the hydraulic gradient may remain very small and difficult to calculate within the noise of the tensiometry.
It follows from the discussion above that the range of soil water pressures where the direct analysis may be useful is also texture (conductivity) dependent. Consider, for example, three initially saturated soil samples of contrasting properties (see Table 1) subject to a lower boundary pressure decrease P of 1 cm min-1. As illustrated in Fig. 3, the accuracy of the direct K(h) and (h) approximations deteriorate as K decreases, particularly in the coarsest soil. The very sharp decrease in K(h) with drainage in the sand (Fig. 3c) results in a large pressure difference across the soil length, signified as h = [h(0,t) - h(L,t)] in the figure, compared with the loam (Fig. 3b) and especially the clay (Fig. 3a). At an average soil water pressure of -40 cm, for example, h is about -45 cm in the sand compared with about -7 cm in the clay. Again, because the pressure gradient is nonlinear with much of the pressure change occurring near the lower boundary, the average pressure as [h(L,t) + h(0,t)]/2 is lower than exists in the majority of the sample. Consequently, the direct estimate of becomes less accurate as |h| increases and overestimation of at the arithmetic h is eventually observed. As a general rule, the deviation between the actual and estimated (h) and K(h) will occur closer to saturation for longer samples, and faster drainage rates and estimation errors are greatest for soils with a narrow pore-size distribution (large n in Eq. [2] or [5]). These limitations of the direct analysis motivate an inverse approach.
Inverse Analysis
Hydrus-1D provides numerical solvers for both direct and inverse variably saturated flow problems. Assuming Richards' equation as governing for water flow, the inverse analysis uses the flux and/or pressure response data (from well-defined boundary perturbations) in a weighted least-squares parameter estimation algorithm for K(h) and (h). Three functional representations for K(h) and (h) may be selected for parameter estimation:
van Genuchten (1980):
 | [2a] |
and
 | [2b] |
where S is the saturated water content, r is the residual water content, KS is the saturated hydraulic conductivity, and 
, m(= 1 - 1/n), and 
are fitting parameters.
Modified van Genuchten (Vogel and Cislerova, 1988):
 | [3a] |
and
 | [3b] |
with
 | [3c] |
 | [3d] |
where m and a are fictitious water contents (m 
S, a 
r), Kk is a measured hydraulic conductivity at (k,hk), and he is the air-entry pressure. While comparably cumbersome, the modified van Genuchten forms provide more accurate representations than the standard equations in soils with a measurable air-entry pressure. In our application of the modified forms, we let a = r, Kk = KS, and k = S, so the conductivity function simplifies to
 | [4] |
Brooks and Corey (1966):
 | [5a] |
and
 | [5b] |
where 
is a fitting parameter termed the pore-size distribution index. For the inverse parameter estimation, if the modified van Genuchten or Brooks and Corey model is selected, the user must deselect the internal interpolation table (used to speed convergence of the numerical solver) by setting the interpolation limits to zero (set the "Iteration Criteria" input file in Hydrus-1D).
The continuous flow approach we use establishes a variable water pressure lower boundary condition, h(0,t), while monitoring (i) the pressure response at the upper boundary, h(L,t), and (ii) the water flux response at the lower boundary, q(0,t). It is straightforward to input these data files into Hydrus-1D (see Simunek et al., 1999, for details). Additional data such as S, KS, or (1.5 MPa) may be included in the objective function. The user chooses the weighting of the data and the method of weighting the data in the objective function. We typically select equal weighting (e.g., 1) of the pressure and flux data, greater weighting (e.g., 10) of directly measured data such as S, and choose weighting by standard deviation for the parameter optimization.
A Simulated Experiment for Inverse Recovery
As a test of the robustness of Hydrus-1D to recover hydraulic properties from the combination of experimental observations suggested, simulated experimental data sets were generated for three test soils (sand, loam, and silty clay loam). For the simulations, a test soil (L = 2.5 cm), initially with h = 0 at z = 0, was subjected to a 1 cm min-1 pressure decrease at the lower boundary for 500 min. Hydrus-1D was used to generate the h(L,t) and q(0,t) records. Since an actual experimental data set will contain a degree of experimental error, the generated q(0,t) and h(L,t) were both contaminated with varying degrees of random error (2, 5, and 15%). The random error was introduced as described by Kool et al. (1985)(see their Eq. [9]). The results of the inverse analysis (Table 2) show that the parameter recovery is remarkably stable with errors in the parameter estimates typically <2%. Neither KS nor S were used as data for the inverse solution. To further challenge the recovery, we intentionally selected somewhat poor initial parameter estimates. For example, the initial estimates for the loam test soil are the Hydrus-1D suggested parameters for a clay. The tabulated results are for the case where we set S at its actual value and floated all other parameters. The case of floating S was also evaluated (not shown) with a finding of a small increase in recovered parameter error. Importantly, the water content parameters, S and r, were the most susceptible to error, and in some cases the recovery of these parameters was uncertain even with zero error in the simulated data. We observed that the inverse routine in Hydrus-1D returns optimized (S–r) and concluded that it is inadvisable to float both S and r, unless a measured (,h) pair such as S and/or (1.5 MPa) is included in the inverse data set. To illustrate the pitfall when a measured is not included in the inverse analysis, consider recovery of the simulated experiment in the silty clay loam test soil (Fig. 4). When simultaneously floating S and r, the recovered (h) has the correct n and values but is shifted on the -axis because of the incorrect S and r found in the inverse analysis. By contrast, recovery of (h) is very good despite 15% random error in the simulated data when S is set at the true value in the inverse analysis. In the parameter estimation, the simulated and optimized q(0,t) and h(L,t) were in excellent agreement in all cases of Fig. 4; r2 = 1.000 when S was floated with zero data error, and r2 = 0.9925 when S was set with the 15% random error data. It is evident (see also Hollenbeck and Jensen, 1998) that an excellent r2 in the inverse recovery does not guarantee accurate representation of (h), and we recommend conditioning the inverse approach with direct measurement of parameters when possible.
View this table:
[in this window]
[in a new window]
|
Table 2. Actual and inverse recovered parameters of van Genuchten (h) and K(h) for three hypothetical soils. The inverse recovery is compared for a simulated continuous flow experiment contaminated with varying degrees of random error (R.E.) in h(L,t) and q(0,t).
|
|
|
RESULTS
|
|---|
Near Saturation Dynamics
The measurement of soil hydraulic properties near saturation is a challenging problem for inverse methods. Boundary effects may lead to large errors in estimation of the soil water flux, and uncertainty as to the air-phase pressure may confuse interpretation of soil water pressure measurements. Further, it is assumed in applying Richards' equation that the air phase is at atmospheric pressure, such that the soil water tension is equivalent to the capillary pressure. A major motivating factor of this work was to develop a method that could be applied at and near saturation (as well as lower pressures) with a minimum of ambiguity, and thereby retain a physical interpretation of the parameter estimates. To elucidate the dynamics of water flow near saturation, we employed a combination of water pressure, air pressure, and water flow measurements as illustrated in Fig. 5 for three soils. The soil air pressure was measured in a small (2-mm diam.) artificial cavity about 15 mm within the soil at z (cm) above the lower boundary of the core. The cavity was created using a small drill and accessed through the sidewall by installing a small (10 mm) length of narrow diameter tubing. In a typical result, as the soil water pressure is decreased at z = 0, the water pressure at the upper boundary parallels the decrease and the air pressure in the artificial cavity decreases. During this period, there is a measurable weight change in the flowcell, which, we argue below, is water loss from the soil–flowcell boundaries. The most interesting feature in Fig. 5 is the nearly simultaneous change in the water pressure and water flux followed shortly by relaxation of the subatmospheric air pressure in the artificial cavity (e.g., note changes at about 6 min in Fig. 5a). Taken together, these observations suggest that air entry at the soil boundary and initiation of soil drainage, presumably from the top down, is signaled by the inflection in the soil water pressure. As observed by Corey and Brooks (1975)(Fig. 4), we see that the soil water tension reaches a critical value necessary to initiate drainage of the largest pores. As air enters these pores, the tension in the surrounding pore water relaxes slightly and an inflection in the tensiometry is evident. Air continuity over the length of the sample may not be instantaneous, as illustrated in the repacked loam with air sensors at three positions (Fig. 5c). The lag in air continuity is less pronounced in structured soil and is related to the air permeability of the soil. Multiple air pressure measurements in the structured sandy clay loam of Fig. 5a, where biopores extended the length of the sample, found simultaneous establishment of atmospheric pressure air pressure near the upper and lower boundaries (not shown). For repacked soil, vents near the lower boundary may hasten air continuity (though an inflection in the water pressure data occurs at air entry with or without venting).
The observation of the sudden air pressure change in the artificial soil cavity is consistent with the air entry or bubbling pressure concept long used in describing soil hydraulic properties (e.g., Brooks and Corey,1966). A challenge here is reconciling the bubbling pressure concept with the notable loss in water from the flowcell (soil-plate-column assemblage) prior to establishment of air continuity. Strictly speaking, if water is treated as incompressible, then the only discharge (albeit small) from the soil prior to the bubbling pressure would be due to deformation of air–water interfaces at the soil boundary as the soil water pressure is reduced (White et al., 1972). For the magnitude of the weight loss we observe, we believe that the apparent contradiction is primarily a boundary artifact of the plate. Through very careful measurement of the flowcell tare weight and the known S of the test soils, we found that the flow-cell assemblage routinely contained more water (as much as 5% by volume) than could be accounted for based on the S of the soil. Moreover, at the air-entry pressure, the of the soil in the flowcell closely matched the S measured independently, suggesting that the apparent soil drainage prior to air entry is actually from boundary voids. To reduce this problem, one may begin (after saturating the soil) with the water table a small distance below the lower plate (the distance used should be less than the air-entry suction head to maintain saturated soil).
The inverse recovered (h) and K(h) for the sandy clay loam of Fig. 5a, a repacked loam, and a repacked fine sand are shown in Fig. 6 together with the directly estimated values. To accomplish the inverse recovery, we first identified the air-entry pressure as discussed above and then, for the Hydrus 1-D input, set q(0,t) = 0 for all times prior to air-entry, while retaining the entire h(0,t) and h(L,t) record. In addition to supplying the known S and (1.5 MPa), we supplemented the inverse data set with = S for a few h he values. We then fixed S and KS at the direct estimates while floating all other parameters (, n, , and r) in Hydrus-1D. The weighting of the = S values in the air-entry region affects the parameter estimates when using the standard van Genuchten models for K(h) and (h); greater weight will force the model to better represent the constant KS and S prior to air entry, but at the expense of the representation at lower pressures. As with the initial parameter estimates, it is generally a case of trial and error adjustment to minimize the objective function in the inverse estimation. We typically assign weights of 1 to 5 to the (h > he) data, weight of 1 on all q(0,t) and h(L,t) data, and weight of 10 on the supplemented s and (1.5 MPa) data.
Figure 6 illustrates the successful application of the inverse parameter estimation to soils at and near saturation using a continuous flow approach. While |h| is small, the agreement between the direct estimates and the inverse recovery is generally very good, particularly for (h). As expected from the numerical simulations (Fig. 3), the direct analysis is a poor estimator as |h| gets large, so it is not alarming that the inverse and direct estimates diverge in the sand where |h| grows to several multiples of the sample length. In general, the modified van Genuchten models are closer to the direct estimates than the standard van Genuchten forms though the differences are slight. An interesting feature to note in Fig. 6 is early in the drainage where open markers are used to symbolize the apparent and K values prior to air continuity in the soil. In this region, we hypothesize that the apparent change in is primarily a boundary artifact as discussed above. As such, the calculated K values are not physically relevant prior to air entry, and the open markers are shown for illustration purposes only. The fluctuations in the apparent K prior to air entry is a common feature and it is not unusual in this region to calculate negative K values. The chaos in the apparent K probably arises from noise in the tensiometry, which causes sign and magnitude oscillations in the small hydraulic gradient. Though the tensiometry noise is small (typically less than ± 0.5 mm pressure head), recall that the drainage begins with an equilibrium profile; that is, h(L) = h(0) - L. Prior to air entry, the zero hydraulic gradient persists since the soil water is stationary but elastic. During this no-flow period, work done at the lower boundary is undiminished by viscous or adhesive forces so that the lower boundary pressure change is mimicked in time at the upper boundary. Even slight deviations in the tensiometry introduce magnitude (positive or negative) to the apparent hydraulic gradient and erroneous fluctuations in the apparent K are observed. Fortunately, it is not necessary to calculate K prior to air entry, but it is sometimes helpful in identifying or confirming the air-entry pressure to look in the K record for a sudden stabilization of values. It is clear that recognition of air entry is important in applying both the direct and inverse analysis, and failure to do so may lead to large errors and great uncertainty in estimation of the hydraulic conductivity near saturation.
The Rate and Range of Pressure Change
In conducting the flow experiments described here, the experimentalist selects the rate of pressure change at the lower boundary and the range of soil water pressures to assess. Aside from the rate issues associated with error in the direct estimates (discussed above), the question arises, To what extent does the drainage rate and range of observation affect the estimation of the hydraulic properties? Because of the variety and complexity of porous material, there is not a simple answer to this question. In the discussion that follows, we suggest why the rate and range of pressure change might affect the hydraulic property estimation and illustrate our points with experimental observations.
Rate Effects
As a practical consideration, the selection of the pressure change rate (P) is bounded only by the tensiometer response time. Unlike the use of manometers, with the null type pressure transducer, there is very little flow across the ceramic and tensiometer response times are very rapid. We tested |P| as high as 8 cm min-1 and observed tensiometer response lag of <1 min. Since we typically use |P| < 2 cm min-1 to retain accuracy in the direct analysis, no correction for tensiometer lag is suggested. Of more relevance, it could be argued that the hydraulic properties are rate dependent, and not just in highly structured soil as intuition might suggest. Topp et al. (1967) compared retention curves obtained by transient, steady-state, and equilibrium methods for a draining sand and found lower water contents for a given potential by the equilibrium approach. Similar rate dependence in (h) of silica beads was reported by Salehzadeh and Demond (1994). Rate dependence in the hydraulic properties of a sand using outflow methods with sudden pressure steps was demonstrated by Wildenschild et al. (2001). One explanation for the rate dependence of hydraulic properties is that a portion of the wetted pore space is isolated or cut-off from the draining pore network. Analogous to the mobile–immobile water concept used in solute transport, the isolated water moves only slowly to the draining network and is not realized as drainage at a particular soil water potential when |P| is large but is observed when |P| is smaller. Surprisingly, there are also reports of exactly the opposite effect, that is, lower water contents with higher rates of drainage. Constantz (1993) measured (h) in a variety of soils and found a 1 to 10% difference depending on drainage rate, with the larger associated with the slower drainage rate. Constantz speculated that the higher drainage rates resulted in lower pore-water salt concentrations and consequently lower pore-water surface tension and lower water contents. Wildenschild et al. (1997) also reported lower water contents with more rapid drainage (in comparing transient and equilibrium methods in a sand and a clay) but attributed the finding to differences in packing between the samples used in the comparison.
Another factor that might lend apparent rate dependence is the frequency of observation relative to the rate of pressure change. Unlike the pore-connectivity effect, this is not based on the physics of the flow system per se but instead is related to the information content of the data. This is probably most critical near saturation, where the functions are most sensitive to changes in pressure. An obvious remedy is to increase the observation frequency with increasing P.
We are currently evaluating rate effects on a variety of soils using the continuous flow approach. Retention curves for two soils with contrasting rates of pressure change at the lower boundary are illustrated in Fig. 7. The retention curves for the silty clay loam are nearly identical over the range of rates evaluated. Some rate dependence appears for the structured field soil and the slight decrease in (at a specific h) with increasing drainage rate is consistent with the findings of Constantz (1993). The results shown for the field soil are from a series of seven rate experiments (showing only the slowest, fastest, and intermediate rates in Fig. 7b) where the results of the inverse parameter estimation (Table 3) were very similar except for a linear trend (r2 = 0.58) in the n parameter of the van Genuchten (h) representation. For example, n = 1.57 and 1.85 at |P| = 0.51 and 2.72 cm min-1, respectively. The decreasing n with decreasing |P| suggests that the slower drainage rate reveals a slightly broader pore-size distribution than the faster rate. For each rate experiment, the pressure and outflow data was recorded at a 1-min interval. We repeated the inverse parameter estimation on the data of the slowest rate case culled to a 3-min sampling interval without affecting the results. Thus, the trend in n does not appear to be an artifact of the observation frequency relative to the change in the soil water potential.
View larger version (25K):
[in this window]
[in a new window]
|
Fig. 7. Effect of drainage rate on the estimated (h) for a nonstructured (repacked) silty clay loam and (b) a structured sandy clay loam. Parameter estimates for each soil and drainage rate are listed in Table 3.
|
|
View this table:
[in this window]
[in a new window]
|
Table 3. Inverse recovered parameters of van Genuchten (h) and K(h) for two soils at contrasting rates of pressure change (P) at the lower boundary.
|
|
In another evaluation of rate effects, we selected two structured field soils (with numerous biopores evident) for a comparison of the continuous flow estimation of (h) with equilibrium measurements. To accomplish this comparison, we followed the continuous flow measurement with resaturation of the soils and drainage to equilibrium at selected potentials. The results (Fig. 8) show lower water contents for the equilibrium measurements compared with the continuous flow result, although the differences are typically less than about 1% volumetric water content. Once again, we see that the inverse and direct estimation of (h) are in good agreement while |h| is small, but the overestimation error of by the direct estimation increases as |h| becomes large.
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 8. The direct, inverse, and equilibrium draining (h) in a structured sandy clay loam field sample (L = 3.5 cm) and a structured clay loam field sample (L = 3.5 cm).
|
|
Range Effects
As an experimental variable, the range of pressure change is bounded by the tensiometer operation, typically reported as no lower than about -700 cm soil water pressure. It is commonly required, however, to estimate (h) and K(h) over the entire range of plant-available water. To accomplish this, the inverse data, h(L,t) and q(0,t), may be supplemented with direct measurements outside the tensiometer range, such as (1.5 Mpa). In Fig. 9, the inverse recovered (h) for a sandy clay loam field soil is compared using data from a single experiment truncated at different soil water pressures. In both Fig. 9a and 9b, we floated r in the inverse analysis (along with and n), but in case of Fig. 9a, we supplemented the inverse data with the independently measured (1.5 Mpa). For this soil, we find that the inverse recovered (h) is nearly identical for the three pressure ranges provided that the (1.5 Mpa) is supplied in the inverse data (Fig 9a). If the (1.5 Mpa) is not supplied (Fig. 9b), the inverse recovered r becomes sensitive to the range of the experimental observation and the curves diverge at the lower soil water pressures. We also found that as the pressure range became smaller, the inverse parameter estimation became more sensitive to the initial parameter estimates.
While Fig. 9 illustrates the benefit of including a (h) measurement outside the tensiometer range, a more dramatic range effect may occur depending on the selection of the hydraulic property model, particularly in a soil with a sizable air-entry pressure. In Fig. 10, the optional model types for (h) and K(h) in Hydrus 1-D are compared for three ranges of pressure change in a repacked silty clay loam. With an air-entry pressure of -32 cm, the (h) and K(h) of this silty clay loam are better represented by the modified van Genuchten and Brooks and Corey models than the standard van Genuchten model. Of the three models, the standard van Genuchten model is the most sensitive to the data range used in the parameter estimation. This occurs because fitting the standard van Genuchten model to the air-entry region compromises the representation at lower soil water pressures. The parameter estimation results are listed in Table 4. The greatest range effect occurs in the estimation of the n parameter in the standard van Genuchten model (n = 6.84 and 2.16 for the 0–50 and 0–300 cm suction ranges, respectively). In all cases, the inverse estimation readily converged and the regression (r2) between the observed and predicted q(0,t) and h(L,t) was very close to 1, indicating excellent agreement. Contributions to the total sum of squares error were not consistently largest for a particular input data type. The large differences in the standard van Genuchten model with experimental range emphasizes the importance of identifying the air-entry pressure and selecting the hydraulic property model accordingly.
View this table:
[in this window]
[in a new window]
|
Table 4. Inverse parameter estimates for the hydraulic functions of a draining silty clay loam as a function of the observation range of the soil water tension. The number N of points in the inverse data increases with range as N = 184 for the 0- to 50-cm range, N = 352 for the 0- to 100-cm range, and N = 1104 for the 0- to 300-cm range.
|
|
|
SUMMARY
|
|---|
The experimental design and analysis presented here provides a simple and rapid method for simultaneous measurement of (h) and K(h). A principal advantage of the continuous flow method as described here is that the combination of measurements allows identification of the air-entry pressure and hence permits application of inverse analysis at and near saturation with minimal ambiguity. Additionally, the direct analysis compliments the inverse approach by providing physically based KS and S that can be used to anchor these parameters in the inverse analysis. As a limitation, the direct K estimate applies to fairly short samples and is susceptible to tensiometer errors, especially systematic errors. Even a small offset in the tensiometers may lead to errors in the estimated K, especially if the error in h at the boundaries is contrasting. For example, in calculating the direct K for the silty clay loam of Fig. 9, a +0.2-cm offset in the lower tensiometer simultaneous with a -0.2-cm offset in the upper tensiometer would increase the estimated KS by a factor of two. Meticulous calibration and positioning of the pressure transducers is necessary for repeatable estimates of K by the direct method, particularly in coarse-textured soils. However, with care, the repeatability of the K estimation is very good. For example, in the series of seven P rates on the same sandy clay loam field sample discussed in conjunction with Fig. 8, the direct estimate of KS was 0.0039 ± 0.0015 cm min-1 (mean ± 95% confidence interval), which is much less than the order of magnitude variation cited as typical by Wildenschild et al. (2001).
Another attractive feature of the continuous flow approach is the exacting control over the rate of pressure change. Unlike single and multistep approaches, which typically employ sudden changes in the air pressure to induce flow, the continuous flow method uses gradual changes in the soil water pressure that more closely mimic most natural flow situations. For the repacked and structured field soils used in this study, the influence of the drainage rate on the estimated (h) was noted, but the effect was typically small. More significant differences in the inversely estimated (h) or K(h) can occur depending on the measured range of soil water pressure and the hydraulic property function selected for the soil. In future work we will further demonstrate the flexibility of the continuous flow approach by looping and nesting pressure changes to measure the main branches of hysteretic hydraulic functions as well as primary and higher-order scanning curves.
|
ACKNOWLEDGMENTS
|
|---|
The authors would like to thank Dr. Art Corey for his insightful discussions regarding the soil water dynamic near air entry. We also wish to express our appreciation to the anonymous reviewers for their discerning and helpful comments.
|
REFERENCES
|
|---|
- Ahuja, L.R., and S.A. El-Swaify. 1976. Determining both water characteristics and hydraulic conductivity of a soil core at high water contents from a transient flow experiment. Soil Sci. 121:198–204.
- Brooks, R.H., and A.T. Corey. 1966. Properties of porous media affecting fluid flow. J. Irrig. Drainage Div., ASCE 72:61–88.
- Constantz, J. 1993. Confirmation of rate-dependent behavior in water retention during drainage in nonswelling porous materials. Water Resour. Res. 29:1331–1334.
- Corey, A.T., and R.H. Brooks. 1975. Drainage characteristics of soils. Soil Sci. Soc. Am. Proc. 39:251–255.
- Durner, W., B. Schultze, and T. Zurmühl. 1999. State-of-the-art in inverse modeling of inflow/outflow experiments. p. 661–681. In M.Th. van Genuchten et al. (ed.) Proc. of the Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. 22–24 Oct. 1997. Univ. of California, Riverside, CA.
- Eching, S.O., and J.W. Hopmans. 1993. Optimization of hydraulic functions from transient outflow and soil water pressure data. Soil Sci. Soc. Am. J. 57:1167–1175.[Abstract/Free Full Text]
- Eching, S.O., J.W. Hopmans, and O. Wendroth. 1994. Unsaturated hydraulic conductivity from transient multistep outflow and soil water pressure data. Soil Sci. Soc. Am. J. 58:687–695.[Abstract/Free Full Text]
- Hollenbeck, K.J., and K.H. Jensen. 1998. Experimental evidence of randomness and non-uniqueness in unsaturated outflow experiments designed for hydraulic parameter estimation. Water Resour. Res. 34:595–602.
- Hopmans, J.W., and J. Simunek. 1999. Review of inverse estimation of soil hydraulic properties. p. 643–659. In M.Th. van Genuchten et al. (ed.) Proc. of the Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. 22–24 Oct. 1997. Univ. of California, Riverside, CA.
- Inoue, M., J. Simunek, J.W. Hopmans, and V. Clausnitzer. 1998. In-situ estimation of soil hydraulic functions using a multi-step soil water extraction technique. Water Resour. Res. 34:1035–1050.
- Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: I. Theory and numerical studies. Soil Sci. Soc. Am. J. 49:1348–1354.[Abstract/Free Full Text]
- Klute, A. 1986. Water retention: Laboratory methods. p. 635–662. In A. Klute (ed.) Methods of soil analysis. Part 1. 2nd. ed. Agron. Monogr. 9. ASA and SSSA, Madison, WI.
- Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity: Laboratory methods. p. 687–734. In A. Klute (ed.) Methods of soil analysis. Part 1. 2nd. ed., Agron. Monogr. 9. ASA and SSSA, Madison, WI.
- Kodesova, R., M.M. Gribb, and J. Simunek. 1998. Estimating soil hydraulic properties from transient cone permeameter data. Soil Sci. 163:436–453.
- Parker, J.C., J.B. Kool, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation. II. Experimental studies. Soil Sci. Soc. Am. J. 49:1354–1359.[Abstract/Free Full Text]
- Salehzadeh, A., and A.H. Demond. 1994. Apparatus for the rapid automated measurement of unsaturated soil transport properties. Water Resour. Res. 30:2679–2690.
- Simunek, J., R. Angulo-Jaramillo, M. Schaap, J.-P. Vandervaere, and M.Th. van Genuchten. 1998a. Using an inverse method to estimate the hydraulic properties of crusted soils from tension infiltrometer data. Geoderma 86:61–81.[ISI]
- Simunek, J., M. Sejna, and M.Th. van Genuchten. 1998b. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat and multiple solutes in variably-saturated media. Version 2.0, IGWMC-TPS70. Int. Ground Water Modeling Center, Colorado School of Mines, Golden, CO.
- Simunek, J., and M.Th. van Genuchten. 1999. Using the Hydrus-1D and Hydrus-2D Codes for estimating unsaturated soil hydraulic and solute transport parameters. p.1523–1536. In M.Th. van Genuchten et al. (ed.) Proc. of the Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. 22–24 Oct. 1997. Univ. of California, Riverside, CA.
- Topp, G.C., A. Klute, and D.B. Peters. 1967. Comparison of water content-pressure head data obtained by equilibrium, steady-state, and unsteady-state methods. Soil Sci. Soc. Am. Proc. 31:312–314.
- van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1994. Inverse method to determine soil hydraulic functions from multistep outflow experiments. Soil Sci. Soc. Am. J. 58:647–652.[Abstract/Free Full Text]
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
- van Genuchten, M.Th., and F.J. Leij. 1992. On estimating the hydraulic properties of unsaturated soils. p. 1–14. In M.Th. van Genuchten et al. (ed.) Indirect methods for estimating the hydraulic properties of unsaturated soils. Univ. of California, Riverside, CA.
- Vogel, T., and M. Cislerova. 1988. On the reliability of unsaturated hydraulic conductivity calculated from the moisture retention curve. Transp. Porous Media 3:1–15.
- Wildenschild, D., K.H. Jensen, K.J. Hollenbeck, T.H. Illangasekare, D. Znidarcic, T. Sonnenborg, and M.B. Butts. 1997. A two-stage procedure for determining unsaturated hydraulic characteristics using a syringe pump and outflow observations. Soil Sci. Soc. Am. J. 61:347–359.[Abstract/Free Full Text]
- Wildenschild, D., J.W. Hopmans, and J. Simunek. 2001. Flow rate dependence of soil hydraulic characteristics. Soil Sci. Soc. Am. J. 63:35–48.
- White, N.F., H.R. Duke, D.K. Sunada, and A.T. Corey. 1972. Boundary effects in the desaturation of porous media. Soil Sci. 113:7–12.
This article has been cited by other articles:
|
|
|
|
 
M. J. Thomasson, P. J. Wierenga, and T. P. A. Ferre
A Field Application of the Scaled-Predictive Method for Unsaturated Soil
Vadose Zone J.,
October 3, 2006;
5(4):
1093 - 1109.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
M. G. Schaap and M. Th. van Genuchten
A Modified Mualem-van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation
Vadose Zone J.,
De |
TOP
ABSTRACT
INTRODUCTION
