Published online 20 November 2007
Published in Vadose Zone J 6:925-934 (2007)
DOI: 10.2136/vzj2006.0141
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: HANFORD SITE
Saturation-Dependent Hydraulic Conductivity Anisotropy for Multifluid Systems in Porous Media
Z. F. Zhang,
M. Oostrom* and
A.L. Ward
Hydrology Group, Environmental Technology Division, Pacific Northwest National Lab., P.O. Box 999, Richland, WA 99352
* Corresponding author (mart.oostrom{at}pnl.gov).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 22 September 2006.
 |
ABSTRACT
|
|---|
The hydraulic conductivity of unsaturated anisotropic soils has recently been described with a tensorial connectivity–tortuosity (TCT) concept. We extend this concept to unsaturated porous media with two or three immiscible fluids. Mathematical expressions to describe the conductivity of each fluid in anisotropic porous media under unsaturated condition are derived in the form of symmetric second order tensors. The theory is applicable to the generalized hydraulic conductivity model and compatible types of saturation–pressure formulation. The extended model shows that the anisotropic coefficient of any one of the fluids is independent of the saturation of other fluids. Synthetic Miller-similar soils generated to be anisotropic were defined by allowing the saturated hydraulic conductivity to have different correlation ranges for different directions of flow. The extended TCT concept was tested using synthetic soils with four levels of heterogeneity and four levels of anisotropy. Numerical experiments of infiltration of two liquid phases, water and the nonaqueous phase liquid (NAPL) carbon tetrachloride, were performed to test the extended model. The results show that, similar to water in a two-fluid (air–water) system, NAPL retention curves in a three-fluid (air–NAPL–water) system were independent of flow direction but dependent on soil heterogeneity, whereas the connectivity–tortuosity coefficients are functions of both soil heterogeneity and anisotropy. The extended TCT model accurately describes unsaturated hydraulic functions of anisotropic soils and can be combined into commonly used relative permeability functions for use in multifluid flow and transport numerical simulations.
Abbreviations: DNAPL, dense nonaqueous phase liquid • NAPL, nonaqueous phase liquid • TCT, tensorial connectivity–tortuosity
|
INTRODUCTION
|
|---|
Nonaqueous phase liquids (NAPLs), such as chlorinated solvents and petroleum hydrocarbons, are found in the subsurface at many contaminated sites. Nonaqueous phase liquids in the subsurface are often long-term contaminants, as their solubilities in water and air are usually relatively low but exceed quality standards. Plutonium recovery operations at the USDOE's Hanford Site in Richland, WA, resulted in organic and aqueous wastes that were released at several locations. The organic wastes consisted of carbon tetrachloride (CCl4) mixed with lard oil, tributyl phosphate, and dibutyl butyl phosphonate. The three major disposal facilities, the 216-Z-9 trench, the 216-Z-1A tile field, and the 216-Z-18 crib, received a total of about 13,400,000 L of liquid waste containing 363,000 to 580,000 L of dense NAPL (DNAPL). Two remediation technologies have been applied. Between 1992 and 2000, about 76,500 kg (48,100 L) of CCl4 was removed using a soil vapor extraction system in the vadose zone. In addition, a pump-and-treat system for the unconfined aquifer removed 4570 kg (2870 L) of CCl4 from groundwater between 1996 and 2000. Last and Rohay (1993) conducted a rough order-of-magnitude estimate of the discharge inventory and computed that 21% of the CCl4 may have been lost to the atmosphere, 12% retained in the vadose zone (dissolved in water, sorbed to the solid phase, and as a component of the gas phase), and 2% dissolved in the saturated zone. The remaining 65% remains unaccounted for in the system.
One of the greatest challenges facing any modeling efforts of DNAPL migration at the Hanford Site and elsewhere is to adequately understand and describe the subsurface flow processes of multiple immiscible fluids, for example, air, NAPL, and water. In past assessments, the hydrogeologic units in Hanford were generally assumed to be either homogeneous and isotropic or homogeneous and anisotropic in the saturated hydraulic conductivity only. In unsaturated porous media, it is often assumed that the anisotropy is the same as at complete saturation (Oostrom et al., 2004; 2006). In reality, these units display complex inter- and intrasedimentary structures at various scales. The effects of these complex structures are generally thought to enhance lateral spreading and impede downward migration (Ward et al., 2006; Stewart et al., 2006). Thus, the effect of these small-scale structures needs to be more thoroughly understood and properly accounted for in the assignment of physical properties to the larger modeled units.
Several researchers have found that for air–water systems, soil anisotropy is dependent on fluid saturation (Mualem, 1984; Yeh et al., 1985; Stephens and Heermann, 1988; McCord et al., 1991; Ursino et al., 2001; Zhang et al., 2003). Zhang et al. (2003) proposed a TCT concept to describe the hydraulic conductivity of anisotropic unsaturated soil. The TCT concept states that in the hydraulic conductivity of unsaturated anisotropic soil, the anisotropy is expressed not merely by a proportionality to the hydraulic conductivity at saturation but also by three connectivity–tortuosity coefficients, Li = (L1, L2, L3), corresponding to the three principal directions. In other words, the TCT concept implies saturation-dependent anisotropy. Raats et al. (2004) presented a mathematical formalization of a connectivity–tortuosity tensor, assuming that the principal axes coincide with those of the hydraulic conductivity tensor at saturation. The hydraulic conductivity tensor of such unsaturated anisotropic soils was shown to be the product of a saturation-dependent scalar variable, the connectivity–tortuosity tensor, and a hydraulic conductivity tensor at saturation.
The purpose of this contribution is to extend the TCT concept of Zhang et al. (2003) and Raats et al. (2004) to saturation-dependent anisotropies for each of the fluid phases in unsaturated porous media with multifluids and to test the concept using synthetic soils with different levels of heterogeneity and anisotropy through numerical simulation of NAPL infiltration.
|
Materials and Methods
|
|---|
Multifluid Flow Governing Equations and Constitutive Relations
The governing equation of mass conservation of an immiscible fluid j can be written as
 | [1] |
where
 | [2] |
and t is time (T), 
is the porosity, 
is the density (M L–3), S is the actual liquid saturation, k is the intrinsic permeability (L2), krj is the relative permeability of phase j, µ is the viscosity (M L–1 T–1), P is the pressure (M L–1 T–2), and gz is the gravitational vector (L T–2). An important aspect of multifluid flow is the mathematical description of the relative permeability (kr), fluid saturation (S), and pressure (P) constitutive relations. The two most widely used multifluid kr–S –P models are the extended van Genuchten (1980) S-P relation in combination with kr–S equations derived from Mualem's (1976) theory, as described by Parker et al. (1987), and the Brooks–Corey (Brooks and Corey, 1964) S-P relation in combination with the kr–S equations derived from Burdine's (1953) pore-size distribution model. Here, although we have used the extended van Genuchten (1980) S–P relation and the kr–S equations derived from Mualem's (1976) theory, the theory is applicable to the combination of any type S–P formulation and the generalized hydraulic conductivity model given below. To allow for straightforward comparisons between this paper and related papers in the soil physics literature, we use water-equivalent pressure heads instead of capillary pressures. The water-equivalent pressure head, hj (L), of a fluid j and capillary pressure head, hij (L), of a nonwetting fluid i and a wetting fluid j are defined as, respectively,
 | [3a] |
and
 | [3b] |
Assuming that fluid wettability follows the sequence water > NAPL > air (Leverett, 1941), the total-liquid saturation (i.e., the sum of water and NAPL saturation) is a function of the air–NAPL capillary pressure, and the aqueous-phase saturation is a function of the NAPL–water capillary pressure. Using the subscript w for water, o for NAPL, and a for air, the extended van Genuchten (1980) relations for effective water, 
w, and total-liquid (subscript t) saturation, t, can be written as, respectively,
 | [4] |
and
 | [5] |
where
 | [6] |
 | [7] |
and ßij are interfacial-tension dependent scaling factors, 
(L–1) and n are curve-shape parameters, m = 1– 1/n, and Sr is the irreducible water saturation. For this application, ßij = 
aw/ij is used (Lenhard et al., 2002), where (M L–2 T–2) is the interfacial or surface tension.
General expressions for the water, oil, and air relative permeability, based on theoretical considerations by Burdine (1953) and Mualem (1976), can be written as
 | [8a] |
 | [8b] |
 | [8c] |
where
 | [9] |
and L is the connectivity–tortuosity coefficient, and 
and 
are constants. With = 1 and = 2, Eq. [8a] is equivalent to the Mualem (1976) relationship for water relative permeability. Substituting the extended van Genuchten (1980) relation (Eq. [4]) into Eq. [8] yields
 | [10a] |
 | [10b] |
 | [10c] |
Multifluid Tensorial Connectivity-Tortuosity Theory
Mualem (1976) found that a connectivity–tortuosity coefficient of L = 0.5 was an optimal value for a data set of 45 disturbed and undisturbed vertically oriented samples for his theoretical model, while a value of L = 2 is often used for the Burdine (1953) model. These values of L have typically been used in describing the permeability of multifluid system (Parker et al., 1987; Parker and Lenhard, 1987; Oostrom and Lenhard, 1998; Lenhard et al., 2004). For air–water systems in anisotropic soils, the parameter L may be directional due to the difference in flow path connectivity–tortuosity in different directions (Zhang et al., 2003; Raats et al., 2004). However, directionally dependent hydraulic property measurements are commonly not conducted, and therefore, experimentally obtained information regarding the directional dependence of hydraulic properties is usually not available.
The TCT concept of Zhang et al. (2003) and Raats et al. (2004) can be extended to multifluid systems. For a three-fluid anisotropic system, fluid permeabilities, defined as the product of soil intrinsic permeability and fluid relative permeability, for the three immiscible fluids in the ith principal direction can be written, respectively, as
 | [11a] |
 | [11b] |
 | [11c] |
If the values of Li are equal in all the directions, Eq. [11] reduce to traditional model of constant anisotropy or isotropy. In the form of symmetric second order tensors, Eq. [11] can be expressed as
 | [12a] |
 | [12b] |
 | [12c] |
where k is the intrinsic-permeability tensor and Tf the connectivity–tortuosity tensor with f being either w, o, or a. These three tensors are defined as
 | [13] |
 | [14] |
 | [15] |
The superscripts 1, 2, and 3 denote the three principal directions. According to Eq. [12], the fluid permeability tensor of an unsaturated anisotropic soil, kf, for each of the three fluids is given as the product of the intrinsic permeability tensor, k, the symmetric connectivity–tortuosity tensor, Tf, and a scalar variable denoting the fluid distribution over the available pore space. Based on the above formulation, the anisotropy coefficient of hydraulic conductivity for fluid phase f between the ith and jth direction, Afij (i = 1, 2, 3; j = 1, 2, 3; i
j) can be written as
 | [16] |
Equation [16] shows that Afij is only dependent on the directional intrinsic permeability and the saturation of fluid of interest. In this paper, we will use superscripts n and p for flow directions normal and parallel to the imposed strata, respectively.
Synthetic Soil Development
The same 16 synthetic soils as those in Zhang et al. (2003) were used in this study to test the saturation-dependent anisotropy in NAPL conductivity. In this section, only the main aspects of the geostatistical, synthetic-soil generation method are discussed. For details of the method, the reader is referred to Zhang et al. (2003, p. 315–316). The van Genuchten (1980) model was used to represent the soil hydraulic properties of a sandy soil with base properties = 4.0 m–1, n = 2.0, = 0.40, Sr = 0.0, and L = 0.5. The saturated hydraulic conductivity, Ks, of this sand was assumed to be 10–4 m s–1, corresponding to an intrinsic permeability, k, of 1.02 x 10–11 m2. The soils were constructed by assigning different levels of heterogeneity and anisotropy to Ks and , while leaving the other variables unchanged. The main statistics of these two variables for the 16 synthetic soils are shown in Table 1
Both soil heterogeneity and anisotropy characteristics are based on the assumption that the soils are "Miller similar" (Miller and Miller, 1956) at the macroscopic scale. On the basis of this assumption, the entire domain may be characterized by a reference state {
h*, K*, h*(
), K*()}, with being the characteristic length, and a unique scaling relationship between the hydraulic functions at different spatial locations, s (Roth, 1995). To generate h(,s) fields, h(s) was treated as a stationary random function with unique probability density and autocovariance functions. An exponential covariance model with correlation lengths along the three principal directions was assumed for the autocovariance function, generated using the sequential Gaussian simulation (SGSIM) program from the GSLIB library (Deutsch and Journel, 1992) with different correlation lengths (
) in the three principal directions. The characteristic length scale was originally based on the mean or median (d50) particle diameter (Miller and Miller, 1956). However, the characteristic length scale can also be based on the average air-entry pressure, which is the pressure at which a nonwetting fluid (typically air) would enter a liquid-saturated porous medium (Rockhold et al., 1996), or the inverse of the van Genuchten parameter (Zhang et al., 2003).
Soil heterogeneity was described by the variance of f, denoted f2, where f is the logarithm of the Miller and Miller (1956) scaling factor. The geometric (denoted by a superscript g) means of the saturated hydraulic conductivity, Ksg, and of the inverse of the capillary length, g, as well as the variances of Y = ln(Ks), denoted Y2, and of A = ln(), denoted A2, were then calculated using the generated Ks and (Table 1). Soil anisotropy was described by the ratio (R) of the correlation lengths parallel (p) and normal (n) to soil strata. To produce 16 different soils, four levels of soil heterogeneity (f2 = 0.1, 0.25, 0.5, and 1.0) and four levels of soil anisotropy (R = 1:1, 10:1, 50:1, and 100:1) were generated on a vertical 1.0 m2 domain uniformly discretized with a grid spacing of 0.02 m. Each of the 2500 soil elements is homogeneous and isotropic but has properties different from other elements.
Numerical Simulations
Using the generated synthetic soils, multistep numerical experiments of gravity-induced flow were performed using the Water-Oil mode of the STOMP simulator (White and Oostrom, 2006). In this mode, the governing equations are those for water and oil component mass conservation, transported over the aqueous phase and NAPL. The aqueous phase comprises liquid water with dissolved oil, while the NAPL comprises a single organic liquid. Equilibrium dissolution is assumed to occur in the presence of NAPL. The Water-Oil mode, which does not include an active air phase, was used to simplify the analyses and presentation of the results. The air phase pressure for this particular mode is kept constant at 101,325 Pa. For the set of simulations in this paper, the van Genuchten (1980) and Mualem (1976) models were used to describe fluid retention curve and fluid conductivity. The NAPL used for the numerical experiments has properties representative of the DNAPL mixture disposed at the Hanford Site at the 216-Z9 trench (Oostrom et al., 2004, 2006). The site NAPL's main component was the chlorinated solvent carbon tetrachloride, with a density of 1426 kg m–3 and dynamic viscosity of 0.00111 Pa s. The assumed surface tensions for the air–water and air–NAPL were 72 dyne cm–1 and 36 dyne cm–1, respectively, while a value of 36 dyne cm–1 was used for the NAPL–water surface tension.
Constant pressure boundary conditions for both water and NAPL were imposed for each step on the top and bottom boundaries, while zero flux conditions were imposed on the other boundaries. By fixing the water and NAPL pressures for each step, capillary pressures, defined as the difference between the nonwetting and wetting fluid pressure at a fluid interface, are constant. Since the fluid pressure values at the top and bottom boundaries were kept the same for each of the two immiscible liquids during each drainage step, gravity was the only driving force when steady-state conditions were obtained.
Each numerical experiment, repeated at four different initial water pressure heads of –2.0, –1.75, –1.5, and –1.25 m, started at the condition at which the total saturation was 1.0. The NAPL in the system was subsequently drained gradually by lowering the NAPL pressure heads in steps ranging from 1 to 5 cm, while keeping the water pressure heads constant, until almost all the NAPL was drained from the domain. The range in NAPL pressure heads is from 0 m to a value slightly higher than the imposed initial water pressure heads. A schematic of the imposed water and NAPL boundary conditions (constant water pressure head but variable NAPL pressure head with time) is shown in Fig. 1
. Simulation of each step continued until steady-state NAPL flow was obtained. The real time for the system to reach steady-state conditions after each boundary condition modification was strongly dependent on the average total liquid saturation of the system and ranged from a few years to several hundred years. After determining the NAPL hydraulic conductivity values in the direction normal to the soil strata, the two-dimensional domain was rotated 90°, similar to the approach in Zhang et al. (2003). The numerical simulations were then repeated to determine the NAPL conductivity in the direction parallel to the soil strata. The total number of simulations equaled 16 (number of soils) x 4 (system water pressures) x 2 (flow directions) = 128.
View larger version (8K):
[in this window]
[in a new window]
|
FIG. 1. Schematic of the water and nonaqueous phase liquid (NAPL) pressure head variation with time at the top and bottom boundaries of the simulation domain. For each simulation set, the water pressure head was kept constant (–2.0, –1.75, –1.5, and –1.25 m) while the NAPL pressure was lowered in steps from zero to a value slightly larger than water pressure head.
|
|
At the completion of each step in a simulation, the mean values of the water pressure, NAPL pressure, water saturation, and NAPL saturation over the computational domain were computed. The ho(St) data pairs were used to fit the effective value of and the effective value of n, while fixing the porosity and the irreducible water saturation. The intrinsic permeability (kp) and connectivity–tortuosity coefficient (Lop) in the direction parallel to the soil strata were then fitted using the kop(ho) data with the already fitted and n treated as constants. Similarly, the parameters in the direction normal to soil strata, kn and Lon, were fitted to the kon(ho) data.
Recent research (Luckner et al., 1989; Vogel et al., 2000; Schaap and Leij, 2000; Schaap et al., 2001; Schaap and van Genuchten, 2005) has suggested several shortcomings of van Genuchten (1980) functions near saturation, most notably the lack of second-order continuity of the soil water retention function at saturation and the inability of the hydraulic conductivity function to account for macroporosity. To avoid these concerns, the obtained data under conditions where the total-liquid saturation was larger than 0.75 were not considered in the analyses of hydraulic conductivities.
|
Results and Discussion
|
|---|
Fluid Retention Curves
Although both water and NAPL contents varied during the numerical simulations, the change of water content was relatively small (
0.05 m3 m–3) because the water pressure head was controlled between –2.0 and –1.25 m. Since the change in NAPL content was relatively large (0.30 m3 m–3), the relationship between the air–NAPL capillary pressure head and total-liquid saturation, as given by Eq. [5], is discussed. The total-liquid retention for each of the 16 soils was obtained by relating the spatially averaged air–NAPL capillary pressure head and total-fluid saturation in both directions. An example of the results for R = 50 is shown in Fig. 2
for various heterogeneity levels. Each plot in this figure includes the measurements for the four different water pressures (–2.0, –1.75, –1.5, and –1.25 m). The results clearly show that the relationship between the capillary pressure head at the air–NAPL interface and total-fluid content was unique for each of the heterogeneous soils, regardless of the initial water saturation, as the measurements under four different water pressures fell to one single curve (Fig. 2).
When the imposed NAPL pressures at the top and boundaries were slightly lower than the gas pressure (i,.e., about –0.01 m water pressure head), the average NAPL pressure over the whole computational domain were often positive. Detailed examination of the spatial distribution of NAPL pressure heads under these conditions shows that positive NAPL pressures generally appeared in relatively finer grained porous materials, while negative pressures occurred in relatively coarser grained materials. Due to the nonlinear effects of soil heterogeneity on NAPL pressure, the average NAPL pressure was positive for all the cases, although the boundary pressure is slightly less than zero. Under similar NAPL pressure head conditions (i.e., about –0.01 m), soil heterogeneity had little effect on total liquid saturation (Fig. 2a–c) except when soil heterogeneity was high and flow was normal to soil strata (Fig. 2d). For the latter case, the relatively finer materials of the soil were saturated while the coarser region was not; hence, on average, the domain was not fully saturated, as shown by the three circles in the upper-right corner of Fig. 2d. Therefore, although all the data are plotted in Fig. 2, the curves were fitted to the data for which the NAPL pressure head were less than zero.
The measurements of fluid retention in two directions also show that the NAPL retention curves obtained when flow was parallel to soil strata were the same as those when flow was normal to soil strata, indicating that the retention curves are independent of flow direction. Similar results were reported by Zhang et al. (2003) for air–water systems and in Mayes et al. (2003), who measured directional soil water retentions of strongly stratified soils using undisturbed soil cores and found that soil water retentions were independent of bedding orientation.
As soil heterogeneity varied, the total liquid retention curve also varied slightly, as shown by different fitted parameters and n (Fig. 2). Compared with the values of and n for air–water systems in Zhang et al. (2003), the values of n for their systems and the air–NAPL–water systems analyzed in this paper are almost identical. The values of for the two systems show similar trends, although with different intercepts (Fig. 3
). Further examinations of the parameter ratios of and n of the air–water system (aw and naw) to those of the air–NAPL–water system (ao and nao) are shown in Fig. 4
. The ratio of ao/aw varied between 1.93 and 2.37, with an average of 2.28. These ratios are very close to the ratio of 2.0 of the interfacial tension of an air–water system (aw = 72 dyne cm–1) to that of an air–NAPL interface (ow = 36 dyne cm–1). These results agree with the findings of Lenhard et al. (2002). The ratio of nao/naw varied between 1.00 and 1.04, with an average of 1.01, indicating that the n parameter is independent of the fluid component.
NAPL Conductivities
After the retention parameters, and n, were determined, the tortuosity coefficient, Lon and intrinsic permeability, kn, were then fitted using the measured NAPL conductivities in the direction normal to soil strata under four different water pressure heads and multiple NAPL pressure heads, whereas Lop and kp were fitted using the data measured in the direction parallel to soil strata. Results for the case with hw = –2.0 m are shown in Fig. 5
. For the isotropic case (R = 1), regardless of soil heterogeneity, the NAPL conductivity curves in both directions were nearly identical, as were the tortuosity–connectivity coefficients. For the anisotropic cases (R = 10, 50, and 100), the NAPL permeabilities in the direction parallel to soil strata, kop, were larger than those in the direction normal to soil strata, kon, for all levels of heterogeneity. The NAPL permeability values in both directions converge at a total liquid saturation of approximately 0.9, suggesting that a minimum anisotropy ratio is obtained near that saturation value. In contrast to the permeability trends, the values of Lop were smaller than Lon for the same soil. As shown in Fig. 6
, the NAPL anisotropic coefficients, defined as kop/kon, increase with decreasing NAPL saturation. Generally, depending on the level of soil anisotropy, the coefficient increased from a value slightly larger than unity at So = 0.6 to seven or more at So = 0.1. However, when compared with the water anisotropy coefficient, reported in Zhang et al. (2003), the saturation dependency of the NAPL anisotropy coefficient is less than the water anisotropy coefficient. This is reflected by smaller differences between the direction tortuosity coefficients for the air–NAPL–water systems than those for the air–water systems.
View larger version (31K):
[in this window]
[in a new window]
|
FIG. 5. Nonaqueous phase liquid (NAPL) permeability as a function of total fluid saturation for anisotropic heterogeneous soils 9–12 in Table 1 and for hw = –2.0 m. Lp, connectivity–tortuosity coefficient in the direction parallel to strata; Ln, connectivity–tortuosity coefficient in the direction normal to strata; R, ratio of correlation lengths of the direction parallel and normal to strata; Y = ln(Ks) with Ks being the saturated hydraulic conductivity; Y2, variance of Y.
|
|
View larger version (27K):
[in this window]
[in a new window]
|
FIG. 6. Nonaqueous phase liquid (NAPL) anisotropy coefficients as a function of NAPL saturation of soils 9–12 in Table 1 for imposed water pressure head of –2.0 m. Lines: computed relations using Eq. [16] with permeabilities listed in Table 2 and L factors values in plots. Circles: ratio of NAPL permeability values for simulations with flow direction parallel to those normal to soil strata.
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 2. Geometric means of intrinsic permeability (kg) and fitted permeability in the direction normal (kn) and parallel (kp) to soil strata.
|
|
The fitted intrinsic permeability, k, can be thought of as the permeability value at full saturation obtained by extrapolation, which may be different from the true value. The fitted value is useful for systems that are not fully saturated, such as the vadose zone in arid or semiarid areas. The fitted intrinsic permeabilities in the directions normal, kn, and parallel, kp, to soil strata for all the cases are summarized in Table 2. As shown by the kp/kg and kn/kg ratios, both the fitted kn and kp for all but one value were slightly larger (by a factor between 1.01 and 3.67) than the geometric mean, kg, of the intrinsic permeability. Generally, both the kp/kg and kn/kg ratios increased as the degree of soil heterogeneity was increased. However, there was no clear effect of correlation length ratio, R, on the kp/kg and kn/kg ratios. The anisotropic coefficient of the fitted intrinsic permeabilities, kp/kn, were 
1 for 11 cases and <1 for 5 cases and varied by a factor no more than 1.5 overall. These results suggest that the fitted intrinsic permeability value for the two directions were very similar and that the geometric mean of a heterogeneous soil may be a good approximation of the fitted kn or kp. The reason for similar fitted kp and kn values may be that, as is shown in Fig. 5b through 5d, at a total liquid saturation of approximately 0.9, the permeability values in both directions converge, indicating that the soils had the smallest anisotropy coefficient at that saturation. For total saturations either larger or smaller than this value, the soil anisotropy coefficients are larger. This phenomenon has also been observed in stratified soils by Mualem (1984). However, under full or near saturated condition, using a fitted k may introduce relative large errors.
Connectivity–Tortuosity Coefficient
The relationships between the connectivity–tortuosity coefficient and the correlation length ratio are depicted in Fig. 7
. The value of L decreased linearly with increasing value of ln(R*), where R* is the ratio of the correlation length of f in the direction parallel to flow to that normal to flow. It is noted that R* is defined relative to flow direction while R was defined relative to the direction of soil strata. Hence, R* = R if flow is parallel to soil strata, and R* = 1/R if flow is normal to soil strata. The results show that L is inversely proportional to ln(R*), with values of R* ranging from 0.01 to 100. Larger values of L indicate less connected pores and/or more tortuous flow paths. Therefore, L was smaller when flow was parallel to soil strata [ln(R*) > 0] than when flow was normal to the soil strata [ln(R*) < 0]. The effect of anisotropy on L is stronger when the soil is more stratified, meaning that the difference between Lp and Ln becomes larger as a soil becomes more stratified. In Fig. 8
, slopes and intercepts of the L vs. ln(R*) relation are plotted. The figure shows that the absolute value of the slope of the regression line increases as the soil became more heterogeneous, suggesting that the effect of soil stratification on L increases when a soil is more heterogeneous. Hence, the value of L is a function of both soil heterogeneity and soil anisotropy.
View larger version (23K):
[in this window]
[in a new window]
|
FIG. 7. Relationships between L vs. ln(R*) of the soils with different levels of heterogeneity. L, tortuosity coefficient; R*, ratio of correlation length in the direction parallel and normal to flow; Y2, variance of Y, where Y = ln(Ks). with Ks being the saturated hydraulic conductivity.
|
|
View larger version (10K):
[in this window]
[in a new window]
|
FIG. 8. Comparisons of the (a) intercepts and (b) slopes of the relationship between L vs. ln(R*) for an air–water (Zhang et al., 2003) and air–NAPL–water system. Y2, variance of Y, where Y = ln(Ks), with Ks being the saturated hydraulic conductivity.
|
|
Although the trends relating L to ln(R*) for an air–NAPL–water systems are in agreement with those found for air–water systems in Zhang et al. (2003), differences are observed in the actual values of Lao and Law (Fig. 8). In general, the values of Lao were smaller than that of Law, meaning that for the same soil heterogeneity, the effect of soil anisotropy was smaller for an air–NAPL–water system than that for an air–water system. These results indicate that tortuosity–connectivity coefficients for water and NAPL are different. This result is consistent with the findings of Tuli et al. (2005), who measured the water permeability and air permeability for several disturbed and undisturbed soil cores and found that the L parameter for the water permeability and air permeability were different. A possible reason is that in an air–water system, water moves in the relatively smaller pores, while in an air–NAPL–water system, NAPL flows in intermediate-size pores. Using a fluid-dependent tortuosity–connectivity coefficient is different from the general practice of using the same parameter value of L for all fluids.
|
Conclusions
|
|---|
The TCT concept of Zhang et al. (2003) and Raats et al. (2004) was extended to porous media with multiple immiscible fluids. Mathematical expressions to describe the conductivity of each fluid in anisotropic porous media under unsaturated conditions were derived in the form of symmetric second order tensors. The theory is applicable to the generalized hydraulic conductivity model and compatible types of saturation–pressure formulation. The extended model shows that the anisotropic coefficient of a fluid is independent of the saturation of other fluids except the fluid of interest. The extended TCT concept was tested using synthetic soils with different levels of heterogeneity and anisotropy and numerical experiments of infiltration of two immiscible liquid phases in synthetic coarse-textured soils. When the connectivity–tortuosity coefficients in different directions are equal, the extended model reduces to the traditional model of constant anisotropy or isotropy.
The numerical experiments showed that for different water and NAPL contents, the total-fluid retention curves obtained when flow was parallel to soil strata were the same as those when flow was normal to soil strata. As soil heterogeneity varied, the total liquid retention curve also varied slightly, reflected in the different effective van Genuchten parameters and n. A comparison of NAPL retention for air–NAPL–water systems with water retention data for air–water systems in the same synthetic soils (Zhang et al., 2003) shows that effective value of parameter n is the same for the two systems (air–water vs. air–NAPL–water), while that of needs to be scaled by the ratio of the surface interfacial tension.
Similar to the air–water model presented by Zhang et al. (2003), the extended model shows that the anisotropy of the NAPL phase is also saturation dependent. For anisotropic soils, the anisotropic coefficients for NAPL increased with decreasing NAPL saturation. However, this increase was generally smaller than that for water in air–water systems.
The connectivity–tortuosity coefficient is a function of both soil heterogeneity and soil anisotropy. The effects of soil stratification on this coefficient increase as a soil becomes more heterogeneous. These effects were not as strong for the air–NAPL–system as those for the air–water system in Zhang et al. (2003). These results indicate that the tortuosity–connectivity coefficient for water and NAPL are different. Further research is warranted to examine the saturation-dependent anisotropy of laboratory porous media and natural soils containing multiple immiscible fluids. In particular, detailed controlled column and intermediate-scale flow cells experiments in heterogeneous systems for air–water and air–NAPL–water systems are needed to provide data sets for comparisons with the Zhang et al. (2003) and the model forwarded in this paper. Successful comparisons between theory and laboratory experiments will support the use of the theory to simulate NAPL infiltration and redistribution at the heterogeneous and anisotropic Hanford Site.
|
ACKNOWLEDGMENTS
|
|---|
Funding for this research was provided by Hanford's Remediation and Closure Science Project. The Pacific Northwest National Laboratory (PNNL) is operated for the U.S. Department of Energy by Battelle under Contract DE-AC05-76RL01830. The fluid properties were obtained in the Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy's Office of Biological and Environmental Research and located at PNNL.
|
REFERENCES
|
|---|
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrology Paper No. 3. Colorado State Univ., Fort Collins.
- Burdine, N.T. 1953. Relative permeability calculations form pore size distribution data. Trans. Am. Inst. Min. Metall. Pet. Eng. 198:71–78.
- Deutsch, C.V., and A.G. Journel. 1992. GSLIB, geostatistical software library and user's guide. Oxford Univ. Press, Oxford, UK.
- Last, G.V., and V.J. Rohay. 1993. Refined conceptual model for the volatile organic compounds: Arid integrated demonstration and 200 West Area carbon tetrachloride expedited response action. PNL-8597. Pacific Northwest Laboratory, Richland, WA.
- Lenhard, R.J., M. Oostrom, and J.H. Dane. 2002. Prediction of capillary pressure-relative permeability relations. p. 1591–1619. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis: Part 4. Physical methods. SSSA, Madison, WI.
- Lenhard, R.J., M. Oostrom, and J.H. Dane. 2004. A constitutive model for air–NAPL–water flow in the vadose zone accounting for immobile, non-occluded (residual) NAPL in strongly water-wet media. J. Contam. Hydrol. 73:283–304.[CrossRef][Medline]
- Leverett, M.C. 1941. Capillary behavior in porous solids. Trans. Am. Inst. Min. Metall. Pet. Eng. 142:152–169.
- Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res. 25:2187–2193.
- Mayes, M.A., P.M. Jardine, T.L. Mehlhorn, B.N. Bjornstad, J.L. Ladd, J.M. Zachara. 2003. Transport of multiple tracers in variably saturated humid region structured soils and semi-arid region laminated sediments. J. Hydrol. 275:141–161.[CrossRef]
- McCord, D., B. Stephens, and J.L. Wilson. 1991. Hysteresis and state dependent anisotropy in modeling unsaturated hillslope hydrologic processes. Water Resour. Res. 27:1501–1518.[CrossRef]
- Miller, E.E., and R.D. Miller. 1956. Physical theory for capillary flow phenomena. J. Appl. Phys. 27:324–332.[CrossRef][Web of Science]
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.[CrossRef]
- Mualem, Y. 1984. Anisotropy of unsaturated soils. Soil Sci. Soc. Am. J. 48:505–509.[Web of Science]
- Oostrom, M., and R.J. Lenhard. 1998. Comparison of relative permeability–saturation–pressure parametric models for infiltration and redistribution of a light monaqueous-phase liquid in sandy porous media. Adv. Water Resour. 21:145–157.[CrossRef]
- Oostrom, M., M.L. Rockhold, P.D. Thorne, G.V. Last, and M.J. Truex. 2004. Three-dimensional modeling of DNAPL in the subsurface of the 216-Z-9 trench at the Hanford Site. PNNL-14895. Pacific Northwest National Laboratory, Richland, WA.
- Oostrom, M., M.L. Rockhold, P.D. Thorne, G.V. Last, and M.J. Truex. 2006. Carbon tetrachloride flow and transport in the subsurface of the 216-Z-9 trench at the Hanford Site: Heterogeneous model development and soil vapor extraction modeling. PNNL-15914. Pacific Northwest National Laboratory, Richland, WA.
- Parker, J.C., R.J. Lenhard, and T. Kuppusamy. 1987. A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res. 23:618–624.[CrossRef]
- Parker, J.C., and R.J. Lenhard. 1987. A model for hysteretic constitutive relations governing multiphase flow: 1. Saturation–pressure relations. Water Resour. Res. 23:2187–2196.
- Raats, P.A.C., Z.F. Zhang, A.L. Ward, and G.W. Gee. 2004. The relative connectivity–tortuosity tensor for conduction of water in anisotropic unsaturated soils. Vadose Zone J. 3:1471–1478.[Abstract/Free Full Text]
- Rockhold, M.L., R.E. Rossi, and R.G. Hills. 1996. Application of similar media scaling and conditional simulation for modeling water flow and tritium transport at the Las Cruces trench site. Water Resour. Res. 32:595–609.[CrossRef]
- Roth, K. 1995. Steady state flow in an unsaturated, two-dimensional, macroscopically homogeneous, Miller-similar medium. Water Resour. Res. 31:2127–2140.[CrossRef]
- 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]
- Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.[CrossRef]
- Schaap, M.G., and M.Th. van Genuchten. 2005. A modified Mualem–van Genuchten formulation for improved description of the hydraulic conductivity near saturation. Vadose Zone J. 5:27–34.[Abstract/Free Full Text]
- Stephens, D.B., and S. Heermann. 1988. Dependence of anisotropy on saturation in a stratified sand. Water Resour. Res. 24:770–778.
- Stewart, M.L., A.L. Ward, and D.R. Rector. 2006. A study of pore geometry effects on anisotropy in hydraulic permeability using the lattice–Boltzmann method. Adv. Water Resour. 29:1328–1340.[CrossRef]
- Tuli, A., J.W. Hopmans, D.E. Rolston, and P. Moldrup. 2005. Comparison of air and water permeability between disturbed and undisturbed soils. Soil Sci. Soc. Am. J. 69:1361–1371.[Abstract/Free Full Text]
- Ursino, N., T. Gimmi, and H. Fluhler. 2001. Combined effects of heterogeneity, anisotropy, and saturation on steady state flow and transport: A laboratory sand tank experiment. Water Resour. Res. 37:201–208.[CrossRef]
- 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.[Web of Science]
- Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2000. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. Adv. Water Resour. 24:133–144.[CrossRef]
- Ward, A.L., Z.F. Zhang, and G.W. Gee. 2006. Upscaling unsaturated hydraulic parameters for flow through heterogeneous anisotropic sediments. Adv. Water Resour. 29:268–280.[Medline]
- White, M.D., and M. Oostrom. 2006. STOMP subsurface transport over multiple phases, version 4.0, user's guide. PNNL-15782. Pacific Northwest National Laboratory, Richland, WA.
- Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985. Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media with variable . Water Resour. Res. 21:457–464.
- Zhang, Z.F., A.L. Ward, and G.W. Gee. 2003. A tensorial connectivity–tortuosity concept to describe the unsaturated hydraulic properties of anisotropic soils. Vadose Zone J. 2:313–321.[Abstract/Free Full Text]
This article has been cited by other articles:
|
|
|
|
 
G. W. Gee, M. Oostrom, M. D. Freshley, M. L. Rockhold, and J. M. Zachara
Hanford Site Vadose Zone Studies: An Overview
Vadose Zone J.,
November 20, 2007;
6(4):
899 - 905.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
