Published online 26 May 2006
Published in Vadose Zone J 5:764-776 (2006)
DOI: 10.2136/vzj2005.0129
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
An Air–Water Interfacial Area Based Variable Tortuosity Model for Unsaturated Sands
Raziuddin Khaleela,* and
K. Prasad Saripallib
a Fluor Government Group, P.O. Box 1050, Richland, WA 99352
b Pacific Northwest National Laboratory, Richland, WA 99354
* Corresponding author (raziuddin_khaleel{at}rl.gov)
Received 9 November 2005.
 |
ABSTRACT
|
|---|
A new variable tortuosity definition is introduced that is based on the immiscible fluid (air–water) interfacial area. Unsaturated media tortuosity (
a) is defined as the ratio of aaw to aaw,o where aaw is the estimated air–water interfacial area in a real unsaturated medium (i.e., a soil sample), and aaw,o is the same variable for the corresponding, idealized capillary bundle. The air–water interfacial area for both real and idealized media is directly proportional to the area under their respective retention curves. With being the saturated tortuosity, we relate the variable tortuosity ratio (/a) to the Se
term in Mualem's ( = 0.5) and Burdine's ( = 2) pore-size distribution models. Thus, instead of using tortuosity and pore connectivity formulations, which have empirical exponents of either 0.5 or 2, the new model depends on a variable interfacial area for varying saturation and soil texture, as reflected in the measured retention data. We tested the new definition of tortuosity to predict unsaturated hydraulic conductivity, K, as a function of volumetric moisture content, 
, for 22 repacked Hanford sediments that are comprised of mostly coarse and fine sands but some also contain a sizeable fraction (as high as 27%) of fines (silt and clay). Replacing the Se term in van Genuchten–Mualem (VGM) model by the tortuosity ratio /a, and still using saturated hydraulic conductivity and moisture retention parameters as used in the conventional approach, we obtained a-based K() predictions that are nearly identical to the conventional VGM model predictions. We also compared the a-based K() predictions with the standard Brooks–Corey–Burdine (BCB) model predictions. In comparison to the VGM model predictions, a-based BCB K() predictions appear to be less biased relative to the measured K for the coarse-textured samples.
Abbreviations: BC, Brooks–Corey capillary pressure desaturation function • BCB, Brooks–Corey–Burdine model • MRC, moisture retention curve • VG, van Genuchten model • VGM, van Genuchten–Mualem model
|
INTRODUCTION
|
|---|
THE PRESENCE of a distinct fluid–fluid interface characterizes the flow of immiscible (air–water) fluids in unsaturated porous media. The interfaces play a key role on the dynamics of multiphase fluid flow and contaminant transport in porous media. The interfacial area per bulk volume of porous medium (aaw, cm2/cm3) is recognized as a fundamental variable necessary to describe the complexities of the pore-scale distributions of immiscible fluids (Gray and Hassanizadeh, 1991). For liquid saturations above residual, several models have been developed to analyze pore-scale saturation distribution and to predict the fluid–fluid interfacial area. These include models based on the pore-scale analysis of spatial distribution of immiscible fluids (Gvirtzman and Roberts, 1991), a pore-scale network model (Reeves and Celia, 1996), capillary tube analogy (Cary, 1994), and thermodynamics (Bradford and Leij, 1997; Gray and Hassanizadeh, 1991; Miller et al., 1990; Morrow, 1970; Leverett, 1941). While the theoretical work has provided new insights, an experimental verification of changes in aaw with changes in relative saturation of immiscible fluids has evolved only recently. The fluid–fluid interfacial area has been measured in the laboratory by surfactant adsorption methods (Karkare and Fort, 1996), light transmission techniques (Niemet et al., 2002; Niemet and Selker, 2001), and interfacial tracer methods (Peng and Brusseau, 2005; Brusseau et al., 2003; Costanza-Robinson and Brusseau, 2002; Anwar et al., 2000; Schaefer et al., 2000; Kim et al., 1999, 1997; Saripalli et al., 1997). Results of various experimental investigations are in general agreement with theoretical studies and, for liquid saturations above residual, suggest an overall decrease of interfacial area with an increase in the wetting-phase saturation (Fig. 1
).
The interfacial area distribution is dictated by the soil texture (Peng and Brusseau, 2005) as well as by the morphology of the pores and solid surfaces that define the pores, and the resulting soil moisture retention characteristics (Leverett, 1941). For a soil sample, the laboratory-measured moisture retention curve (MRC) is an effective index of pore size distribution as well as the thermodynamic status of fluid retention (Leverett, 1941; Morrow, 1970). Fluid flowing through unsaturated media is forced to traverse through water films, whose expanse of interfacial area is larger than that of a corresponding idealized medium (Fig. 2
). In this study, we used the ratio of interfacial area of the laboratory-measured retention curve to the corresponding retention curve for the idealized medium as a quantitative measure of the additional, tortuous path that the fluid must traverse through the real soil. The objectives of this study were to: (i) develop a modeling approach for the prediction of variable tortuosity based on estimated air–water interfacial areas (aaw); and (ii) test the model's ability to predict saturation-dependent tortuosity for sandy media. In particular, we compared our variable tortuosity model predictions with those based on the empirical pore connectivity and tortuosity models of Mualem (1976a) and Burdine (1953). A comparison of the measured and predicted unsaturated hydraulic conductivity, K, as a function of volumetric moisture content, for 22 coarse-textured sediments (Khaleel et al.,1995) provided a test, albeit indirect, of our proposed variable tortuosity model. While aaw is considered a fundamental variable, very little work has been reported on application of aaw-based models to predict tortuosity and K(). The interfacial friction offers resistance to flow; it is therefore reasonable to expect that the relative permeabilities of immiscible fluids in unsaturated media are functions of aaw (Silverstein and Fort, 2000; Gvirtzman and Roberts, 1991). On the basis of estimated liquid–vapor interfacial configurations for varying matric potentials, Tuller and Or (2002) developed a methodology to calculate liquid saturation and to predict sample-scale unsaturated hydraulic conductivity. They represented the pore space by a bimodal distribution of pore sizes, reflecting two disparate populations of matrix and structural pores.
View larger version (48K):
[in this window]
[in a new window]
|
Fig. 2. (a) Schematic of a real medium and the corresponding idealized medium (image courtesy of D. Or and M. Tuller); (b) matric potential as a function of effective saturation for the real and the idealized media.
|
|
Note that the method presented here considers only the main drainage cycle of the desaturation curve; therefore it does not consider hysteresis in moisture retention. While the interfacial area based variable tortuosity concept was advanced by Saripalli et al. (2002), to our knowledge, its application to predicting K() via aaw-based tortuosity estimates and comparison with actual measurements is new. As detailed below, our approach to deriving the variable tortuosity on the basis of aaw is also different from that of Saripalli et al. (2002).
|
THEORETICAL CONSIDERATIONS
|
|---|
Variable Tortuosity and Unsaturated Hydraulic Conductivity Models
Tortuosity for a porous medium is typically defined as (Le/L)2, where Le is the length of the tortuous flow path and L is that of the straight flow path in the idealized capillary bundle. No methods exist for the measurement or prediction of Le or L. In a direct analogy to (Le/L)2, and based on immiscible fluid interfacial areas, which can be measured and estimated, Saripalli et al. (2002) proposed a new definition for unsaturated media tortuosity, a:
 | [1] |
where aaw is the estimated immiscible fluid (air–water) interfacial area in a real unsaturated medium, and aaw,o is the same variable for the corresponding idealized capillary bundle. The maximum interfacial area corresponds to a condition in which pendular rings have formed (i.e., relatively low water contents such as residual water content). As the water content increases above residual water content, the interfacial area decreases (Fig. 1).
Figure 2a is a schematic of a real soil sample and the corresponding idealized capillary bundle. As shown below, the interfacial area for both real and idealized media is directly proportional to the area under their respective moisture retention curves. While the idealized medium contains many of the same characteristics of the real sample, an important distinction is the absence of tortuosity in the idealized capillary bundle shown in Fig. 2a. This makes the use of such an idealized model an attractive choice to characterize tortuosity that is present in a real soil sample. As suggested by Fig. 2b, for the same water content, the corresponding matric potentials for the real and idealized media are different.
We now relate the a term in Eq. [1] to its analog in the Mualem (1976a) and Burdine (1953) pore-size distribution models. We describe the moisture retention data by either the Brooks–Corey (Brooks and Corey, 1964) or the van Genuchten (van Genuchten, 1980) model, whereas the Mualem (1976a) and Burdine (1953) models are used to predict K().
Mualem (1976a) represented the porous medium as being randomly distributed, interconnected pores with radius r and described by a pore radii distribution f(r), with the smallest pore radius being rmin and the largest being rmax. Mualem introduced two factors: a factor G(,r,q) to account for the partial correlation between pore radii r and q at a water content and a factor T(,r,q) to account for tortuosity. Both T and G are functions of only. Mualem's relative hydraulic conductivity, Kr() model is given by
 | [2] |
where Kr = K()/Ks and Ks is the saturated hydraulic conductivity. The power of 2 of the quadratures in Eq. [2] evolves from the assumption that the pore configuration can be replaced by a pair of capillary elements whose lengths are proportional to their respective radii, r and q. Mualem suggested replacing the correction and tortuosity factors by a single factor Se0.5 where Se = ( – r)/(s – r) = the effective water saturation (0 
Se 1), s = volumetric moisture content at saturation, and r = volumetric moisture content at residual saturation. Using the capillarity equation, the soil moisture retention curve is used to describe the pore radii distribution:
 | [3] |
where h is matric potential, which for notational convenience is considered as being positive (i.e., h denotes tension), 
w is density of water, g is acceleration due to gravity, 
is the fluid surface tension, and 
is the contact angle. For water at 20°C and with being 0°, Eq. [3] is simply
 | [4] |
where h and r are in centimeters and 0.149 is in square centimeters. Combining Eq. [2] and [4] with Se0.5, Mualem obtained
 | [5] |
In contrast to Mualem (1976a), Burdine (1953) presented the functional dependence of tortuosity of the wetting phase in an unsaturated medium as
 | [6] |
where T1.0 is the tortuosity of the wetting phase when Se is 1.0 and TSe is the tortuosity as a function of Se. The TSe term in Eq. [6] is analogous to a in Eq. [1].
The Burdine (1953) and Mualem (1976a) models are special cases of the generalized model (Hoffmann-Riem et al., 1999):
 | [7] |
where the parameters , ß, and 
determine the shape of the conductivity function and are equal to 2, 2, and 1, respectively, for the Burdine model and 0.5, 1, and 2 for the Mualem model. According to Carman (1937), the tortuosity for a saturated granular medium is 
2. With Se interpreted in terms of pore connectivity and tortuosity, Se should always be <1, since 0 < Se < 1. Similar to Se, we let the tortuosity ratio, /a, approach 1 as Se approaches 1. Note that at Se = 1, aaw theoretically is zero; a is calculated only as Se 1.
The Brooks–Corey (BC) capillary pressure desaturation function in terms of Se is (Brooks and Corey, 1964)
 | [8] |
where hb is the air-entry value or bubbling pressure (cm) and is considered to be positive, and 
is a parameter characterizing the pore-size distribution of the soil. The van Genuchten (VG) model in terms of Se is (van Genuchten, 1980)
 | [9] |
where 
and n are empirical fitting parameters ( is in 1/cm, n is dimensionless), and m is approximated as m = 1 – 1/n.
By combining Eq. [8] and [9] with Eq. [7], we obtain the conventional BCB (Eq. [10]), and VGM (Eq. [11]) predictive models. These models use moisture retention measurements and Ks to predict K:
 | [10] |
and
 | [11] |
Furthermore, recognizing the similarity of the interfacial-area-based tortuosity ratio /a (Eq. [1]) with Se in Eq. [7], we obtain, from Eq. [10] and [11], the interfacial-area-based BCB (Eq. [12]) and VGM (Eq. [13]) equations:
 | [12] |
and
 | [13] |
The variable tortuosity ratio /a is now embedded in Eq. [12] and [13], thus allowing a comparison of K() predictions based on the BCB and VGM models.
Interfacial Area for the Real Medium
Relating the change in interfacial free energy to the externally imposed capillary pressure, Leverett (1941) presented an expression, based on thermodynamic principles, for the interfacial area created during drainage between successive saturation states in an initially water-saturated medium:
 | [14] |
where Sw = /s, Aaw is the interfacial area per unit pore volume, and G is the Gibbs free energy. Expressing the interfacial area in Eq. [14] relative to bulk volume of the porous medium and changing the variable from Sw to Se for the integral:
 | [15] |
where Sr = r/s is residual saturation and 
h(Se)
real represents the laboratory-measured MRC for the real soil. The aaw term in Eq. [15] represents the interfacial area of the real medium for a in Eq. [1]; the integral represents the area under the laboratory-measured drainage curve and can be evaluated either analytically or numerically. For the real medium, we obtained the VG parameters s, r, , and n and the BC parameters s, r, hb, and via curve fitting of laboratory-measured moisture retention data; we used numerical integration to obtain the integral aaw between any two values of Se in Eq. [15]. Analytical solutions for the integral are also available for both VG and BC models (Niemet et al., 2002). Both numerical and analytical methods yielded near-identical results. Note that the calculated aaw values (Eq. [15]) on the basis of h vs. Se and rescaled h* vs. Se (i.e., h vs. Se for the VG model and h/hb vs. Se for the BC model) are identical.
Interfacial Area for the Idealized Medium
To determine aaw,o, we followed the formulation developed by Cary (1994). Consider a porous medium having a s and a cross-sectional area A that is idealized as a cylindrical capillary bundle of straight (nontortuous), parallel tubes of different radii (Fig. 2). Each tube is continuous, of the same length, and contains no dead end or stagnant region. All pores with a radius r less than r' are completely filled with water, whereas those equal to or greater than r' are lined with a water film but are otherwise empty (Cary, 1994). The latter pores have a distinct air–water interfacial area. Note that, whereas in the real soil, the flow boundaries are a combination of solid–water and solid–air interfaces, the wetting-phase boundaries of the idealized medium are composed entirely of solid–water interfaces with empty capillaries constituting the nonwetting phase. The incremental interfacial area for the idealized medium is (Cary, 1994; Eq. [3])
 | [16] |
To obtain the total air–water interfacial area (aaw,o) for the bundle of capillary tubes, we integrated the incremental area associated with pores having r > r' (which are all lined with a water film), invoked the capillarity Eq. [3], and changed the variable from Sw to Se for the integral:
 | [17] |
where aaw,o is the area per unit bulk volume of the porous medium. Thus, although Cary used a different approach, the air–water interfacial area for the idealized medium is still related to the integral under the retention curve, h(Se)ideal, for the idealized medium. Equation [17] is equivalent to the thermodynamics-based Eq. [15] containing the h(Se)real curve for the real soil.
Idealized Medium Moisture Retention Curve
Compared with the idealized medium (Fig. 2b), for the same Se, the interfacial area is larger for the real medium. Therefore, as indicated in Fig. 2b, for the same Se, the h value for the idealized medium is smaller than the h value for the real medium. The inverse of interfacial area per bulk volume of porous medium (L–1) can be viewed as a scaling factor, as described below.
Consider the capillary rise, hc, in a tube, with the height of capillary rise given by Eq. [4]. With the idealized medium consisting of a bundle of straight, vertical capillary tubes and with varying tube radii following a Gaussian distribution (Kutilek and Nielsen, 1994), the equilibrium Se–h relation (with h = hc) for the idealized bundle will also follow a Gaussian distribution.
For an idealized medium for sandy soils, Cary (1994) proposed a power function relation (power b = 2) for the radius ri of tubes emptying corresponding to a given Sw:
 | [18] |
where rmax is the radius of the largest tube.
Equation [18], however, is not adequate in representing the h(Se)ideal across the entire range from a complete saturation to residual saturation. The shape of the Se–h curve for the idealized medium is best represented by the VG model (Eq. [9]):
 | [19] |
where the subscript ideal is for the idealized medium. With its right side strictly a function of Se for a constant n, Eq. [19] can be viewed in the context of scaling for similar media (Miller and Miller, 1956; Warrick, 2003). Because tortuosity is present in the soil sample but not in the idealized capillary bundle, there exists for the real medium a length scale real, which is greater than the length scale ideal for the idealized medium. For a given Se, h values for the two media are related so that
 | [20] |
Thus h(Se)real for the real and h(Se)ideal for the idealized media have the same functional form and differ only by a scaling factor. With the dimension of in Eq. [19] in inverse length, for a given h:
 | [21] |
For the samples analyzed in our study, i.e., coarse sandy samples having a narrow pore-size distribution, we used m = 0.5 (Warrick, 2003) and n = 2 in Eq. [19]. Note that unlike the BC model, the VG model does not explicitly contain an air-entry pressure term; it assumes drainage to occur with h = 0. For the VG model for the idealized medium, we therefore used a large value to represent rmax (i.e., h = 0). Using Eq. [4], recognizing that during drainage the pores drain sequentially from the largest to the smallest, and with VG being proportional to rmax, Eq. [19] with n = 2 can also be written analogous to Eq. [18]:
 | [22] |
The retention curve for the idealized medium stays below the curve for the real medium (Fig. 2b). For the BC model for the idealized medium, we related the VG n to the BC (i.e., = n – 1) and set hb equal to 1/ for the VG model. Note that, because a is calculated as a ratio of areas under the rescaled retention curves (i.e., h vs. Se for the VG model and h/hb vs. Se for the BC model), /a calculations are insensitive to the choice of and hb for the idealized medium. Once the h(Se) curve for the idealized medium is defined, the h() curve for the same medium is obtained by setting r and s for the idealized medium equal to those for the real medium.
For the idealized medium, Eq. [17], in conjunction with the VG or BC retention models on the basis of Eq. [19], is used to calculate numerically the integral aaw,o, which represents the interfacial area of the idealized medium for a in Eq. [1]. The integrals aaw (Eq. [15]) and aaw,o (Eq. [17]) are then used in Eq. [1] to obtain a as a function of Se. Again, we let the ratio /a approach 1 as Se approached 1. Equations [1]
through [19] represent a complete description of the two closed-form relations for characterizing /a and the K() relationship across an effective saturation range of 1 to 0 for nonhysteretic soils.
|
MATERIALS AND METHODS
|
|---|
We tested the interfacial-area-based models using the experimental data reported in Khaleel et al. (1995), who obtained moisture retention and unsaturated hydraulic conductivity data in the laboratory for a total of 22 (repacked) sediment samples. The particle-size distribution and the median particle sizes for the 22 samples are given in Khaleel et al. (1995; Table 1). The moisture retention data for the drainage cycle up to –1000 cm of pressure head were measured using Tempe pressure cells; the rest of the drainage curve up to –15 000 cm was measured using the pressure plate extraction method (Klute, 1986). Saturated hydraulic conductivities for the samples were measured in the laboratory using a constant-head permeameter. A variation of the steady-state head control method (Klute and Dirksen, 1986) was used to measure K() in the laboratory (Khaleel and Heller, 2003). Additional details about the characteristics of the sediment samples, experimental procedures, and analysis of data can be found in Khaleel et al. (1995). The fitted BC parameters for the 22 samples are shown in Table 1. The fitted VG parameters (m = 1 – 1/n) for the 22 samples are available in Khaleel et al. (1995; Table 2).
|
RESULTS AND DISCUSSION
|
|---|
To predict K on the basis of the BCB and VGM models, the fitted moisture retention parameters for the 22 samples were used to estimate aaw for each sample, using Eq. [15]. The interfacial area aaw,o for the idealized capillary bundle was based on Eq. [17], combined with Eq. [19] for the idealized-medium MRC. Subsequently, using Eq. [1], a was calculated as a function of effective saturation, /a set to 1 (for Se 1), and substituted in Eq. [12] and [13], to obtain the K() relationship for a given sample.
Clearly, an important aspect of this work is the use of Eq. [1] in defining a based on interfacial area. To our knowledge, no methods are available for the direct measurement of unsaturated media tortuosity as a function of water content. Nonetheless, the new definition of tortuosity is attractive and has physical basis, as the interfacial area for unsaturated media can be quantified using interfacial tracers (e.g., Kim et al., 1997; Saripalli et al., 1997).
Variable Tortuosity and Unsaturated Hydraulic Conductivity Predictions for Individual Samples
Brooks–Corey–Burdine Models
As discussed above, Eq. [12] and [13] allow comparison between the BCB and VGM model predictions using the variable tortuosity ratio of /a. For the BCB model, Fig. 3
shows the /a vs. Se plots for all 22 samples for both a-based tortuosity and Burdine tortuosity (Se2) models. As stated above, the /a estimates in Fig. 3, as well as elsewhere, were obtained by letting /a = 1 as Se approached 1.
As Fig. 3 suggests, for most of the samples, the tortuosity ratio /a (solid line) deviates considerably from the Burdine model tortuosity (Se2) predictions (broken line). The differences between the two curves in Fig. 3 translate to deviations between the K() predictions using the standard BCB model (Eq. [10]) and the a based BCB model (Eq. [12]). Figure 4
presents the standard BCB model as well as the a based BCB model predictions, in addition to K() measurements, for all 22 samples. In general, the a-based BCB model predictions compare well with K() measurements for the entire range of . The majority of the samples in Fig. 4 are relatively coarse textured, with more than 90% of the particle-size distribution representing the sand fraction. The a-based model predictions compare well with K() measurements for nearly all of the coarse-textured samples that contained no fine fraction (i.e., silt and clay). Several samples, however, contain a sizeable fine fraction that ranged from 4 to 27% (Khaleel et al., 1995; Table 1). Sample 0-072 (Fig. 4), for example, suggests a somewhat poor fit of the measured and predicted K; the poor fit might be related to the presence of a relatively high fine fraction (19%) for this particular sample. Nonetheless, the majority of the a-based predictions are well within an order of magnitude of the K() meaurements (Fig. 4). Although not apparent from Fig. 4, in general, predictions based on the standard BCB model (Eq. [10]), at low , deviated somewhat from the measured K. Such a bias in predictions is discussed below for BCB models.
View larger version (31K):
[in this window]
[in a new window]
|
Fig. 4. Comparison of variable tortuosity (a) based Brooks–Corey–Burdine predictions (solid line) with the standard Brooks–Corey–Burdine predictions (broken line) and unsaturated hydraulic conductivity, K, measurements (squares) for 22 samples (the triangles represent the measured saturated hydraulic conductivity, Ks).
|
|
Van Genuchten–Mualem Models
For the identical samples shown in Fig. 3, Fig. 5
shows, for the VGM model, the /a vs. Se plots for all 22 samples for both a-based tortuosity and Mualem (Se0.5) models. Again, similar to BCB models, the tortuosity ratio /a predictions deviate from the Mualem model Se0.5 predictions. Similar to Fig. 4, Fig. 6
presents the standard model (Eq. [11]) and the a-based model (Eq. [13]) predictions as well as the K() measurements for all 22 samples. For the majority of the samples, K() predictions based on the standard model and the a-based VGM model are nearly identical (Fig. 6). Compared with the a-based BCB model predictions (Fig. 4), deviations between the measured and predicted K (Fig. 6) are greater for a number of samples (e.g., Samples 2-1636, 2-1637, and 1-1419). The worst fit is provided by Sample 1-1419 (Fig. 6); predictions are underestimated by as much as three orders of magnitude.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 6. Comparison of variable tortuosity (a) based van Genuchten–Mualem predictions (solid line) with the standard van Genuchten–Mualem predictions (broken line) and unsaturated hydraulic conductivity, K, measurements (squares) for 22 samples (the triangles represent the measured saturated hydraulic conductivity, Ks).
|
|
The standard VGM and BCB model predictions use considerably different formulations. For the Burdine model, = 2, ß = 2, and = 1. For the Mualem model, = 0.5, ß = 1, and = 2. In particular, the parameter (
in Mualem's notation) was estimated as 0.5, an average of some 45 samples. Like our database, Mualem's database consisted of a large number of relatively coarse-textured repacked sediments (Mualem, 1976b); however, while the average was 0.5, fitted values for different samples ranged from about –5 to 5. Schaap and Leij (2000) also found considerable variability in the parameter . Again, at low , the a-based as well as the standard VGM predictions deviate considerably from K() measurements (Fig. 6).
The general agreement in the trend of predictions on the basis of Eq. [12] and [13] suggests that the basic idea of deriving a as a function of the interfacial area is reasonable for the coarse-textured samples considered in our study. Thus, while an independent experimental characterization of the a–Se relationship is not available, the results presented in Fig. 3 through 6 do support the idea of conceptualizing unsaturated media tortuosity as a function of interfacial area. The differences in K() predictions are due to differences in the model formulation. Note that the predictions by the two models apply across an effective saturation range of one to zero. Equation [14], which is based on Leverett's thermodynamic prediction of the air–water interfacial area's dependence on Sw may not be strictly valid for saturations below residual, i.e., Sw < Sr (Leverett, 1941) because of the possible absence of a continuous wetting phase.
Comparison of Measured and Predicted Unsaturated Hydraulic Conductivity
While a qualitative comparison of measured K with predicted K for individual samples looks promising (Fig. 4 and 6), no quantitative measures were provided above on the overall quality of fit in terms of goodness of fit and presence of potential bias in the predictions. Here, we provide the goodness of fit of K() measurements with the a-based K() predictions as well as with the standard predictive models. The goodness of fit for all 22 samples is described by the residual mean squared deviation (RMSD) between the measured and predicted K values; logarithmic values of K are used to avoid bias toward high conductivities in the wet range (Schaap and Leij, 2000). The degree of bias in the data is characterized by the deviation of the fitted linear regression slope to the desired slope of 1:1 as well as by testing the null hypothesis, at the 0.01 probability level, that the mean difference between measured and predicted K values is zero (Conover, 1980). Table 2 presents the summary statistics.
Brooks–Corey–Burdine Models
Figures 7a
and 7b show, respectively, a comparison of the measured and predicted K for all 22 samples using the a-based BCB model (Eq. [12]) and the standard BCB predictive model (Eq. [10]). Note that the standard BCB predictive K values are based on the laboratory-measured saturated conductivity and the moisture retention data.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 7. (a) Variable tortuosity (a) based Brooks–Corey–Burdine model predictions and unsaturated hydraulic conductivity, K, measurements and (b) the standard Brooks–Corey–Burdine model predictions and K measurements (the solid line represents the regression line and the dotted line represents the 1:1 line).
|
|
While the scatter in predictions (Fig. 7a and 7b) and r2 values are of similar magnitude (Table 2), it is apparent that the two data sets behave differently. Unlike the standard BCB model predictions (Fig. 7b), the a-based BCB model predictions (Fig. 7a) compare well with the measured K at low water contents, with the majority of the predictions being within 1 to 1.5 orders of magnitude of measured K. With the slope of the regression line being 0.98 (Table 2), the trend line for the BCB model predictions (Fig. 7a) compares well with the 1:1 line and no significant bias is apparent for either high or low . To the contrary, for the standard BCB model (Fig. 7b), most of the predictions are below the 1:1 line with a regression slope of 1.27 (Table 2), which is an indication that the estimates are biased. Also, as indicated in Table 2, unlike the a-based BCB model predictions (Fig. 7a), the standard BCB model predictions (Fig. 7b) rejected the null hypothesis that the mean difference between the measured and predicted K values is zero, again suggesting that the estimates are biased.
As was shown in Fig. 4, when compared with the coarse-textured samples, the fit between the measured and predicted K is slightly poor for the fine-textured samples. This is supported by observations in the literature that, in general, the standard BCB model is known to produce relatively accurate results for coarse-textured soils, similar to ours, that are characterized by a relatively narrow pore- or particle-size distribution. Results are known to be less accurate for many fine-textured and undisturbed field soils because of the absence of a distinct air-entry value for these soils (van Genuchten et al., 1991). Nonetheless, a poor fit is provided by the standard BCB model, with an increase in bias in predicted K with a decrease in (Fig. 7b).
van Genuchten–Mualem Models
Figures 8a
and 8b show, respectively, scatter plots of the measured and predicted K for all 22 samples using the a-based VGM model (Eq. [13]) and the standard VGM predictive model (Eq. [11]). The RMSD and r2 values (Table 2) are very similar for Fig. 8a and 8b. Compared with the a-based BCB model predictions (Fig. 7a), the a-based VGM model (Fig. 8a) as well as the standard VGM model predictions (Fig. 8b) are noticeably biased across the entire moisture regime. Differences between the measured and predicted K for both models appear to grow with a decrease in , and nearly all K() predictions lie below the 1:1 line. The a-based as well as the standard VGM model predictions (Fig. 8) rejected the null hypothesis that the mean difference between the measured and predicted K values is zero (Table 2).
View larger version (20K):
[in this window]
[in a new window]
|
Fig. 8. (a) Variable tortuosity (a) based van Genuchten–Mualem model predictions and unsaturated hydraulic conductivity, K, measurements and (b) the standard van Genuchten–Mualem model predictions and K measurements (the solid line represents the regression line and the dotted line represents the 1:1 line).
|
|
Clearly, on the basis of comparison of predicted K() for the BCB and VGM models, some overall trends are emerging. The a-based BCB model predictions appear to best fit the measured K() for the 22 samples, with no significant bias in predicted values. A somewhat greater bias is apparent especially for low in both a-based and standard VGM model predictions than is present in a-based BCB model predictions. Results suggest that, when compared with measured K, the a-based Burdine tortuosity model predictions are more accurate than the a-based Mualem tortuosity model predictions for the coarse-textured soils considered in our study.
Our method derives the air–water interfacial area from the area under the Se–h drainage curve, as parameterized using closed-form soil moisture retention functions. As discussed above, interfacial area in unsaturated media can be measured using interfacial tracers (Saripalli et al., 1997); however, our method can be used without this additional measurement. This suggests that the new model has physical basis. In any case, it should be noted that measuring interfacial area using interfacial tracers is much less tedious and time consuming than measuring unsaturated conductivity, especially for low , which for our samples took an average of 5 wk per sample. Thus it is encouraging that the K estimates on the basis of the new model, the laboratory-measured retention data, and saturated conductivity compare well with the measured K; the a-based K() predictions are at least as good as those based on the conventional BCB and VGM models.
|
CONCLUSIONS
|
|---|
While the fluid–fluid interfacial area aaw is generally recognized as a fundamental variable in describing flow of immiscible fluids in unsaturated media, very little work has been reported on the application of aaw-based models to predict tortuosity and unsaturated K. A variable tortuosity for unsaturated media (a) is defined as aaw/aaw,o (ratio of the estimated air–water interfacial area for the real and the corresponding idealized porous media). Unlike the Burdine (Se2) and Mualem (Se0.5) empirical models, the new a-based model parameters have explicit physical significance. For the coarse-textured sediments, the a-based variable tortuosity formulation provides at least as good a fit of the measured and predicted conductivities as the standard BCB and VGM models. Specifically, BCB K() predictions that include the newly developed a-based variable tortuosity compare well with measured data for repacked sands. While results are promising, further testing is needed using undisturbed sediments and other soil types (e.g., gravelly and fine-textured media).
|
ACKNOWLEDGMENTS
|
|---|
The work reported was performed for the U.S. Department of Energy Office of River Protection under Contract DE-AC06-99RL14047 with CH2M Hill Hanford Group, Inc. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof or its contractors or subcontractors. The views and opinions of the authors do not necessarily state or reflect those of the U.S. Government or any agency thereof.
|
REFERENCES
|
|---|
- Anwar, A.H.M.F., M. Bettahar, and U. Matsubayashi. 2000. A method for determining air–water interfacial area in variably saturated porous media. J. Contam. Hydrol. 43:129–146.
- Bradford, S., and F. Leij. 1997. Estimating interfacial areas for multi-fluid soil systems. J. Contam. Hydrol. 27:83–105.
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol. Pap. no. 3. Colorado State Univ., Ft. Collins.
- Brusseau, M.L., N.T. Nelson, and M.S. Costanza-Robinson. 2003. Partitioning tracer tests for characterizing immiscible-fluid saturations and interfacial areas in the vadose zone. Vadose Zone J. 2:138–147.[Abstract/Free Full Text]
- Burdine, N.T. 1953. Relative permeability calculations from pore-size distribution data. Trans. AIME 198:71–78.
- Carman, P.C. 1937. Fluid flow through granular beds. Trans. Inst. Chem. Eng. (London) 15:150–166.
- Cary, J.W., 1994. Estimating the surface area of fluid phase interfaces in porous media. J. Contam. Hydrol. 15:243–248.
- Corey, A.T. 1977. Mechanics of heterogeneous fluids in porous media. Water Resour. Publ., Littleton, CO.
- Conover, W.J. 1980. Practical nonparametric statistics. John Wiley & Sons, New York.
- Costanza-Robinson, M.S., and M.L. Brusseau. 2002. Air–water interfacial areas in unsaturated soils: Evaluation of interfacial domains. Water Resour. Res. 38:1195, doi:10.1029/2001WR000738. (Correction: Water Resour. Res. 39:1091. doi:10.1029/2003WR002018.)[CrossRef]
- Gray, W.G., and S.M. Hassanizadeh. 1991. Unsaturated flow theory including interfacial phenomena. Water Resour. Res. 27:1855–1863.[CrossRef]
- Gvirtzman, H., and P.V. Roberts. 1991. Pore scale spatial analysis of two immiscible fluids in porous media. Water Resour. Res. 27:1165–1176.[CrossRef]
- Hoffmann-Riem, H., M.Th. van Genuchten, and H. Fluhler. 1999. General model of the hydraulic conductivity of unsaturated soils. p. 31–42. In M.Th. van Genuchten, F.J. Leij, and L. Wu (ed.) Proc. Int. Worksh. Characterization and measurement of the hydraulic properties of unsaturated porous media, Riverside, CA. 22–24 Oct. 1997. Univ. of California, Riverside.
- Karkare, M.V., and T. Fort. 1996. Determination of the air–water interfacial area in wet "unsaturated" porous media. Langmuir 12:2041–2044.[CrossRef]
- Khaleel, R., J.F. Relyea, and J.L. Conca. 1995. Evaluation of van Genuchten–Mualem relationships to estimate unsaturated hydraulic conductivity at low water contents. Water Resour. Res. 31:2659–2668.[CrossRef]
- Khaleel, R., and P. R. Heller. 2003. On the hydraulic properties of coarse-textured sediments at intermediate water contents. Water Resour. Res. 39:1233, doi:10.1029/2003WR002387.
- Kim, H., P.S.C. Rao, and M.D. Annable. 1999. Gaseous tracer technique for estimating air–water interfacial areas and interface mobility. Soil Sci. Soc. Am. J. 63:1554–1560.[Abstract/Free Full Text]
- Kim, H., P.S.C. Rao, and M.D. Annable. 1997. Determination of effective air–water interfacial area in partially saturated porous media using surfactant adsorption. Water Resour. Res. 33:2705–2712.[CrossRef]
- Klute, A., 1986. Water retention: Laboratory methods. p. 635–662. In A. Klute (ed.) Methods of soil analysis. Part 1. 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. Agron. Monogr. 9. ASA and SSSA, Madison, WI.
- Kutilek, M., and D.R. Nielsen. 1994. Soil hydrology. Catena Verlag, Cremlingen, Germany.
- Leverett, M.C. 1941. Capillary behavior in porous solids. Trans. Am. Inst. Min. Metall. Pet. Eng. 142:152–169.
- Miller, C.T., M.M. Poirier-McNeill, and A.S. Mayer. 1990. Dissolution of trapped nonaqueous phase liquids: Mass transfer characteristics. Water Resour. Res. 26:2783–2796.[CrossRef]
- Miller, E.E., and R.D. Miller. 1956. Physical theory for capillary flow phenomena. J. Appl. Phys. 27:324–332.[CrossRef][ISI]
- Morrow, N. 1970. Physics and thermodynamics of capillary action in porous media. Ind. Eng. Chem. Res. 62(6):32–56.
- Mualem, Y. 1976a. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.[CrossRef]
- Mualem, Y. 1976b. A catalogue of the hydraulic properties of unsaturated soils. Res. Proj. Rep. no. 442. Technion, Israel Inst. of Technology, Haifa.
- Niemet, M.R., M.L. Rockhold, N. Weisbrod, and J.S. Selker. 2002. Relationships between gas–liquid interfacial surface area, liquid saturation, and light transmission in variably saturated porous media. Water Resour. Res. 38:1135, doi:10.1029/2001WR000785.[CrossRef]
- Niemet, M.R., and J.S. Selker. 2001. A new method for quantification of liquid saturation in 2-D translucent porous media systems using light transmission. Adv. Water Resour. 24:651–666.[CrossRef]
- Peng, S., and M.L. Brusseau. 2005. Impact of soil texture on air–water interfacial areas in unsaturated sandy porous media. Water Resour. Res. 41:W03021, doi:10.1029/2004WR003233.[CrossRef]
- Reeves, P.C., and M.A. Celia. 1996. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour. Res. 32:2345–2358.[CrossRef]
- Saripalli, K.P., H. Kim, P.S.C. Rao, and M.D. Annable. 1997. Measurement of specific fluid–fluid interfacial areas of immiscible fluids in porous media. Environ. Sci. Technol. 31:932–936.
- Saripalli, K.P., R.J. Serne, P.D. Meyer, and B.P. McGrail. 2002. Prediction of diffusion coefficients in porous media using tortuosity factors based on interfacial areas. Ground Water 40:346–352.[Medline]
- Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
- Schaefer, C.E., D.A. Dicarlo, and M.J. Blunt. 2000. Experimental measurement of air–water interfacial area during gravity drainage and secondary imbibition in porous media. Water Resour. Res. 36:885–890.
- Silverstein, D.L., and T. Fort. 2000. Incorporating low hydraulic conductivity in a numerical model for predicting air–water interfacial area in wet unsaturated particulate porous media. Langmuir 16:835–838.[CrossRef]
- Tuller, M., and D. Or. 2002. Unsaturated hydraulic conductivity of structured porous media: A review of liquid configuration-based models. Vadose Zone J. 1:14–37.[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.[ISI]
- van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC Code for quantifying the hydraulic functions of unsaturated soils. EPA/600/2–91/065. USEPA, Washington, DC.
- Warrick, A.W. 2003. Soil water dynamics. Oxford Univ. Press, Oxford, UK.
TOP
ABSTRACT
INTRODUCTION
for the 22 samples in 
