HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Tuller, M. Articles by Or, D. Search for Related Content PubMed Articles by Tuller, M. Articles by Or, D. GeoRef GeoRef Citation Agricola Articles by Tuller, M. Articles by Or, D. Related Collections Hydraulic Conductivity Structure and Properties Preferential Flow

REVIEWS & ANALYSES

Unsaturated Hydraulic Conductivity of Structured Porous Media

A Review of Liquid Configuration–Based Models

^a Department of Plant, Soil & Entomological Sciences, University of Idaho, Moscow, Idaho 83844-2339
^b Department of Plants, Soils & Biometeorology, Utah State University, Logan, Utah 84322-4820
* Corresponding author (mtuller{at}uidaho.edu)

Received 27 December 2001.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Common approaches for modeling hydraulic functions of unsaturated structured porous media (SPM) rely on macroscopic continuum representation, where parameterization schemes and constitutive relationships originally developed for homogeneous porous media are extended to represent hydraulic behavior of dual (or multi) continuum SPM. Such models often result in inconsistencies due to lack of consideration of structural pore space geometry and the neglect of underlying physical processes governing liquid retention and flow under unsaturated conditions. We review a new framework that considers equilibrium liquid configurations in dual continuum pore space as the basis for calculation of liquid saturation and subsequent introduction of hydrodynamic considerations. The SPM pore space is represented by a bimodal distribution of pore sizes, reflecting two disparate populations of matrix and structural pores. Three steady-state and laminar flow regimes are considered to derive unsaturated hydraulic conductivity functions: (i) flow in completely filled pore spaces, (ii) corner flow in partially filled pores and grooves, and (iii) film flow on solid surfaces. Two key assumptions are used in deriving the average cross-sectional flow velocities in these regimes: (i) that equilibrium liquid–vapor interfaces remain stable under slow laminar flows and (ii) that flow pathways are parallel. Liquid–vapor interfacial configurations for different matric potentials are calculated and statistically upscaled to derive sample-scale saturated and unsaturated hydraulic conductivity from velocity expressions weighted by the appropriate liquid-occupied cross-sectional areas, neglecting three-dimensional (3-D) network effects. Similarly, the hydraulic functions for matrix and structural pores are derived separately and later combined by weighting the individual contributions by the porosities of the associated pore spaces. A parameter estimation scheme was developed to calculate liquid saturation and to predict sample-scale unsaturated hydraulic conductivity. Model evaluation using measured data for homogeneous porous media, fractured welded tuff, and macroporous and aggregated soils shows favorable agreement (within the limitations of model assumptions). Effects of nonequilibrium conditions between matrix and structural pore domains on the hydraulic conductivity and approximate consideration of 3-D network effects are discussed.

Abbreviations: BCC, bundle of cylindrical capillaries • CPA, critical path analysis • FPM, fractured porous media • SPM, structured porous media • TCwt, Tiva Canyon welded tuff • VGM, van Genuchten–Mualem [model] • 3-D, three-dimensional • 2-D, two-dimensional

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
FLOW AND TRANSPORT in partially saturated SPM takes place in two or more vastly different pore domains, each with distinct pore-size distribution, shape, porosity, and topology. Flow regimes in matrix and structural pore spaces coexist and interact within the same medium. Rapid flow in structural pore space tends to bypass the bulk matrix and penetrates to large depths, leading to preferential transport of chemicals and contamination of water resources. Therefore, development of physically based hydraulic functions for SPM is an essential first step towards quantitative description of preferential flow phenomena in fractured rock and aggregated and macroporous soils.

Structured porous media consist of interconnected networks of matrix and structural pores forming two (or more) distinct pore spaces. Typically, spaces of sizes in the order of 10^-4 to 10^-2 m are associated with structural pores, whereas the porous matrix contains smaller pore sizes in the range of 10^-7 to 10^-5 m. The resulting large disparity in hydraulic behavior between structural pores and matrix in unsaturated SPM presents practical and theoretical challenges to modeling of total system response. In contrast to relatively well-developed theory and experimental data for flow behavior in saturated SPM, theory and measurements for unsaturated SPM are limited. Most modeling approaches employ macroscopic continuum representation where parameterization schemes and constitutive relationships developed for homogeneous porous media are extended to SPM (Mohanty, 1999; Durner, 1994). Such approaches often fail to capture localized and time-varying flow phenomena, and may result in inconsistencies due to oversimplified pore space representation as a bundle of cylindrical capillaries (BCC), and the lack of a physical basis for their effective hydraulic parameters. Liquid retention due to adsorptive forces and flow in thin liquid films are often neglected, leading to poor predictive capabilities of unsaturated hydraulic conductivity at low saturation.

In the absence of consensus regarding the spatial scales, flow behavior, and types of pore space heterogeneity amenable to continuum representation, there is a clear need for physically sound constitutive relationships between matric potential, liquid retention, and hydraulic conductivity. Recent progress in modeling hydraulic behavior of partially saturated homogeneous porous media (Tuller et al., 1999; Or and Tuller, 1999; Tuller and Or, 2001b) provides the basis for assembling such relationships for arbitrary pore space geometry using fundamental solid–liquid physical interactions and basic flow behavior. This approach could bridge some gaps between volume-averaged macroscopic representation of flow properties and discrete representation of structural pores in SPM, and offer a framework where the total system response is derived from fundamental physical processes (instead of nonphysical extension of standard porous media parameters). Nevertheless, the approaches reviewed in this study are based on certain assumptions leading to existence of steady-state flow regimes in such SPM. Evidence shows that under certain circumstances flow patterns in structural pore spaces tend to be spatially irregular (fingering) and temporally intermittent. Such phenomena are not amenable to simple representation and often involve extreme sensitivity to initial and boundary conditions, which cannot lead to constitutive relationships in the standard sense. These flow regimes are subject of ongoing research (Or and Ghezzehei, unpublished data, 2002; Su et al., 1999).

The primary objective of this study was to review constitutive hydraulic functions for SPM based on equilibrium liquid configurations as a function of matric potential (µ), considering steady-state flow regimes within simplified geometrical representation of media pore space. The specific objectives were to review: (i) new conceptual models for pore space geometry of fractured, aggregated, and macroporous media; (ii) derivation of closed-form functions for liquid retention in matrix and structural pore spaces (dual-continuum), based on liquid configuration as a function of µ; (iii) the use of equilibrium liquid configurations to introduce hydrodynamic considerations leading to derivation of unsaturated hydraulic conductivity functions for SPM (ignoring network effects); and (iv) simple tests of the proposed models using available datasets.

The review is organized as follows. In Theoretical Considerations we introduce pore-scale geometrical representations for matrix and structural pore spaces, review recent developments in modeling hydraulic properties of homogeneous porous media (matrix pores) (Tuller and Or, 2001b), and illustrate the derivation of functions for liquid saturation and unsaturated hydraulic conductivity for two-dimensional (2-D) cross sections of structural pores, considering corner, film, and parallel-plate flow regimes. Subsequently, we introduce a statistical upscaling scheme for representation of sample-scale response of the SPM using a statistical distribution of elements and combining flow processes in the two pore domains (matrix and structural), assuming hydraulic equilibrium (i.e., no interactions). In the Results and Discussion section, we apply the proposed model using experimental data for homogeneous porous media, fractured welded tuff, and aggregated and macroporous soils. We examine aspects of nonequilibrium between flow domains and their potential effect on the constitutive relationships. Finally, we discuss approximations for inclusion of 3-D network effects in the proposed 2-D scheme.

	THEORETICAL CONSIDERATIONS

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Pore Geometry (Unit Pore Space Element)
Matrix Pore Space
Even cursory inspection of soil or rock micrographs reveals that natural pore spaces do not resemble cylindrical capillaries (Fig. 1). This is because many natural porous media are formed by aggregation of primary particles and various mineral surfaces, with the resulting pore space more realistically described by angular or slit-shaped pores rather than by cylindrical capillaries (Li and Wardlaw, 1986; Mason and Morrow, 1991).

View larger version (149K): [in this window] [in a new window]   Fig. 1. (a) Thin section of Devonian Sandstone (Roberts and Schwartz, 1985) revealing angular pore space. (b) Scanning electron micrograph (SEM) of calcium-saturated montmorillonite.

View larger version (26K): [in this window] [in a new window]   Fig. 2. Unit cell for pore space representation comprising an angular central pore attached to slit-shaped spaces.

Structural Pore Space
In this review we show examples for fractured, aggregated, and macroporous media. Due to its abundance and higher complexity we will illustrate derivations for fractured porous media (FPM), and comment on their application to aggregated and macroporous media. Or and Tuller (unpublished data, 2002) illustrated that a 2-D fracture network can be represented as an assembly of basic fracture elements, each comprised of two parallel surfaces separated by a certain aperture size (B). Each surface contains a single groove (or pit) representing surface roughness (Or and Tuller, 2000). The fracture element length (B), and groove depth (B), are assumed to be proportional to fracture aperture size (B), as depicted in Fig. 3a. Although such geometrical scaling imposes constraints on the nature of surface roughness, it simplifies the model and facilitates derivation of closed-form expressions for FPM hydraulic conductivity. Nevertheless, the proposed geometry shows reasonable versatility as illustrated in Fig. 3b through 3e, using several combinations of scale parameters for mated and unmated fracture surfaces. Moreover, scaling constraints can be relaxed as more information on surface roughness becomes available (at the expense of requiring a numerical scheme for evaluation of the hydraulic functions).

View larger version (31K): [in this window] [in a new window]   Fig. 3. Definition sketch for (a) a unit fracture element representing a partially saturated fracture with liquid retained in crevices and adsorbed liquid films; (b–e) various combinations of scale parameters and surface mating based on the same geometrical definitions as in (a).

Equilibrium Liquid Configuration within Basic Unit Elements
Equilibrium liquid–vapor interfacial configurations evolve with changes in matric potential and determine the liquid-occupied cross-sectional areas (i.e., saturation) within the unit elements representing matrix and structural pore spaces discussed above. The desaturation of initially saturated elements by a gradual decrease of matric potential is not a continuous process, but rather involves spontaneous liquid displacement (formation of separated interfaces) at certain critical potentials as determined by liquid properties and geometry. Before proceeding with identification of these critical potentials at the transition between various element filling stages, we first address a few basic relationships between adsorbed liquid film thickness and liquid–vapor interfacial curvature as functions of matric potential.

Tuller et al. (1999) developed a model for liquid configuration and liquid saturation in partially saturated homogeneous porous media that serves as the basis for subsequent derivations. The framework is based on two complementary elements: (i) a unitary approach for explicit consideration of the individual contributions of adsorptive and capillary forces to the matric (chemical) potential and (ii) its implementation within the pore space geometry elements introduced in the previous section. Modern interface science formalism and the concept of the disjoining pressure (Derjaguin et al., 1987) provided the physicochemical basis for incorporation of adsorption phenomena into the augmented Young-Laplace equation (Derjaguin, 1957; Philip, 1977; Novy et al., 1989; Blunt et al., 1995). Simplifications of the rigorous Young-Laplace equation have been instrumental in the development of practical closed-form expressions for both saturation and interfacial area as a function of matric (chemical) potential at the pore scale.

While adsorptive and capillary forces are considered simultaneously in the augmented Young-Laplace equation, their contribution to liquid retention is treated separately in the simplified approach, termed the shifted Young-Laplace equation (Tuller et al., 1999). For a given geometry and matric potential, we first calculate the thickness of adsorbed films h(µ) lining all exposed surfaces within the unit elements, and then shift the radius of interface curvature r(µ) obtained from the classical Young-Laplace equation by the film thickness h(µ). Figure 4 illustrates the application of the SYL to a unit cell representing matrix pores.

View larger version (59K): [in this window] [in a new window]   Fig. 4. Conceptual sketch for the simplified approach (SYL). The radius of interface curvature r(µ) is shifted by the thickness h(µ) of the adsorbed liquid film. Note that r(µ) and h(µ) are functions of the matric potential µ.

[1]

where A_svl is the Hamaker constant for solid–vapor interactions through the intervening liquid, and is the density of the liquid. The matric potential is expressed in terms of energy per unit mass (J kg^-1) and may be converted to pressure (energy/volume) (Pa) by multiplication with water density (). Equation [1] represents the simplest form of liquid adsorption on solid surfaces considering van der Waals forces only (Derjaguin et al., 1987). The radius of liquid–vapor interface curvature (r) for unsaturated conditions is dependent on matric potential (µ) according to the Young–Laplace relationship:

[2]

where is the surface tension of the liquid. An interesting feature of this uniform capillary radius of curvature is that all central pores with the same shape retain the same amount of liquid-filled cross-sectional area regardless of their size.

Both here and in subsequent derivations, we assume thermodynamic equilibrium between structural and matrix pore domains (an assumption that was verified experimentally in some fractured porous media systems [Wang and Narasimhan, 1993]).

Spontaneous Liquid Displacement (Snap-Off) During Drainage and Imbibition
Liquid configuration (and saturation) within unit element cross sections evolve with changes in matric potential. Imbibition and drainage processes of unit elements are not continuous processes, but rather show "jumps" at certain critical potentials where liquid interfaces spontaneously and abruptly attain new configurations. Such spontaneous displacement of liquid is termed snap-off and depends on the size and shape of the unit elements. In the following we will discuss critical snap-off potentials for unit matrix and fracture elements. Critical potentials derived for the central pore in matrix unit cells are also employed for triangular, circular, and polygonal pore cross sections representing macropores and structural pores in aggregated media.

In the transition from adsorption- to capillary-dominated imbibition of matrix unit cells and vice versa (drainage) (Fig. 5), four cell-filling stages and two major snap-off mechanisms are considered. At relatively dry conditions (low matric potentials), thin liquid films are adsorbed on all surfaces of the unit cell, and liquid accumulates in corners of the central pore as a result of capillary forces (Fig. 5a). As the matric potential gradually increases to a certain critical value, the two separate films within the slits become unstable and bridge the gap. This is followed by a spontaneous filling of the slits (assuming the films are connected with the bulk liquid phase) (Fig. 5b). This critical potential µ_M-1i for slit snap-off at imbibition is given as (Iwamatsu and Horii, 1996):

[3]

where A_llv is the Hamaker constant for liquid–liquid interactions through the intervening vapor phase. The conditions for slit snap-off during imbibition including time scales for onset of film instability were recently discussed by Restagno et al. (2000). A further increase in matric potential leads to a further increase of film thickness and to an increase of the radius of capillary interface curvature in the central pore, until the capillary corner menisci make contact to form an inscribed circle (Lenormand et al., 1983) (Fig. 5c). At this critical potential, liquid fills up the central pore spontaneously (pore snap-off) (Fig. 5d). The critical potential µ_M-2i is calculated according to Mason and Morrow (1986) and Tuller et al. (1999):

[4]

with A as the cross-sectional area of the central pore, P as the pore perimeter (P = nL, with n as the number of corners of the central pore), and F_n as a pore shape dependent angularity factor given as (Tuller et al., 1999)

[5]

where is the corner angle of the regular polygonal central pore. This defines the curvature coefficient C_ni:

[6]

View larger version (49K): [in this window] [in a new window]   Fig. 5. A sketch illustrating liquid–vapor interfacial configurations in a unit cell with triangular central pore during transition from adsorption to capillary-dominated imbibition. (a) Partially saturated unit cell with adsorbed films and liquid accumulated in corners of the central pore, (b) liquid configuration after slit snap-off, (c) liquid saturation at pore snap-off when capillary menisci form an inscribed circle, and (d) completely liquid saturated unit cell after pore snap-off.

[7]

with C_nd as the drainage radius of curvature coefficient,

[8]

There is ambiguity regarding the critical matric potential for slit snap-off during drainage. It is possible to define a critical slit spacing (L*) by (Eq. [9]) that separates slit sizes according to capillary drainage (Derjaguin and Churaev, 1976) vs. adsorption-dominated (Eq. [3]) drainage. Capillary based slit snap-off would be applied to slit spacing greater than L* and adsorption-based snap-off for spacing smaller or equal to L*.

[9]

We noticed that while the introduction of such critical slit spacing has a minor effect on the resulting saturation behavior, it greatly complicates calculations dealing with a population of unit cells (see Upscaling to Sample Scale—Statistical Pore-Size Distribution below). For simplicity we choose to apply Eq. [3] uniformly for drainage conditions as well.

In unit fracture elements we are confronted with slightly different conditions. We distinguish between three element filling stages (Fig. 6) in the transition from complete saturation (high matric potential) to dry conditions (low matric potential). Starting with a completely saturated fracture element cross section and lowering the matric potential gradually, we will reach a certain potential threshold value µ_F-2d, where the fracture spontaneously empties and two separate liquid–vapor interfaces will form on the opposite fracture surfaces. The critical potential µ_F-2d at this point of separation is derived from capillarity considerations:

[10]

where B is the fracture aperture. The presence of roughness on the fracture surface ensures that some liquid is retained (by capillary forces) in pits and surface grooves. The radius of interface curvature of a meniscus anchored at the edges of the pit is simply B/2 (or -/µ_F-2d) at the separation potential. For certain pit depths (parameterized by ), such as,

[11]

the radius of interface curvature at fracture evacuation results in menisci that are tangent to the pit surfaces, which greatly simplifies subsequent calculations. In cases where the inequality given in Eq. [11] is not satisfied, we have to introduce a second potential threshold (µ_F-1d) for mathematical tractability, which marks the starting point for recession of capillary menisci into the surface pit. The critical potential µ_F-1d is obtained from simple geometrical considerations:

[12]

where is the pit angle (Fig. 3). Hence, for a given geometry, we first evaluate Eq. [11] and calculate the relevant critical potentials. The drainage potentials derived above also hold for imbibition of the unit fracture elements, unless the fracture aperture is very small and adsorptive forces start dominating the system. As with considerations for matrix elements, introduction of a critical aperture size separating capillary and adsorption dominated drainage and imbibition processes will greatly enhance mathematical complexity with only negligible effects on saturation behavior. We therefore apply Eq. [10] and [12] uniformly for drainage and imbibition.

View larger version (16K): [in this window] [in a new window]   Fig. 6. A sketch illustrating liquid configurations and critical potentials during fracture drainage. (a) Three-step transition for geometries where the capillary meniscus is first anchored at the pit edges after interface separation and then recedes into the surface pit (note that the second transition was introduced for mathematical tractability). (b) Two-step transition for geometries where the capillary meniscus immediately tangents the pit walls after interface separation.

For potentials µ_M-2 > µ µ_M-1 in matrix elements, where the central pore is partially saturated and the slits are completely liquid-filled (Fig. 5b), degree of saturation S_wM-2(µ) is calculated according to Tuller et al. (1999):

[13]

where n is the number of corners of the central pore and A_n is an area factor for the central pore defined as (Tuller et al., 1999):

[14]

For potentials smaller than µ_M-1 (partially saturated pore and slits), S_wM-3(µ) is given as:

[15]

Similar algebraic expressions were derived for liquid–vapor interfacial area, an important medium property for interfacial exchange processes (e.g., oxygen exchange and bioremediation) (Tuller et al., 1999).

For a given fracture element geometry, we first evaluate Eq. [11] and calculate the relevant critical potentials (Or and Tuller, unpublished data, 2002). If the geometry requires introduction of µ_F-1, the relative saturation curve is obtained by employing the following expressions. As mentioned above, for all potentials µ µ_F-2, the unit element is completely saturated and the relative saturation S_wF-1(µ) is simply 1. For all potentials µ_F-2 > µ µ_F-1, where capillary menisci are anchored at the pit edges, relative saturation is given as:

[16]

where is a dimensionless fracture length scaling parameter. Note that the denominator in Eq. [16] is the cross-sectional area of the unit fracture element without solid shell. Finally, for all potentials µ < µ_F-1 where menisci start receding into the pit corners and additional film covered area is exposed, relative saturation is calculated as:

[17]

where F is an angularity factor defined as:

[18]

For all cases where capillary menisci are tangent to the pit surfaces immediately after interface separation (i.e., where Eq. [11] is satisfied), S_wF-1 = 1 for µ µ_F-2, and Eq. [17] is employed for all potentials µ < µ_F-2.

The calculated equilibrium liquid configurations in unit element cross sections provide the basis for introduction of hydrodynamic considerations within the unit element geometries of matrix and structural pores.

Hydrodynamic Considerations for a Unit Cell
The detailed picture of liquid–vapor interface configuration and the underlying assumption that equilibrium liquid–vapor interfaces remain relatively stable under slow laminar flow conditions enable introduction of hydrodynamic aspects. Several other studies (e.g., Ransohoff and Radke, 1988; Blunt and Scher, 1995; Dillard and Blunt, 2000) have used the equilibrium radius of interface curvature r(µ) to define boundaries of flowing liquid. Calculations of restraining capillary forces in small open channel capillaries also show that break-up or significant deformation of the liquid–vapor interface requires relatively large viscous forces (Romero and Yost, 1996), especially at lower potentials. We use the equilibrium liquid–vapor interfaces as fixed boundaries to define laminar pore-scale flow regimes. In subsequent derivations we ignore network effects (pore topology and connectivity) and we assume parallel flow paths in direction perpendicular to the unit cell cross section—similar to derivations used in the BCC model (Fatt and Dykstra, 1951; Millington and Quirk, 1961; Mualem, 1976a). On the basis of filling stage, we can identify four primary flow regimes in matrix pore elements (Fig. 7a; numbers following refer to the figure). For complete cell saturation, liquid flows in (1) ducts (central pore) and between (2) parallel plates (slits); For partial saturation, and after pore and slit snap-off, flow occurs in (3) thin adsorbed films lining all flat parts of the cell, and in (4) corners of the central pore (Fig. 8a).

View larger version (28K): [in this window] [in a new window]   Fig. 7. Conceptual sketch of considered flow regimes under various cell-filling stages in (a) matrix, and (b) fracture unit elements.

View larger version (57K): [in this window] [in a new window]   Fig. 8. Three-dimensional representation of corner and film flow regimes within (a) a partially saturated matrix unit cell and (b) a unit fracture element.

With boundaries defined by the equilibrium liquid–vapor interfaces, we can formulate Navier–Stokes equations for each flow regime and solve for average flow velocity. Average hydraulic conductivity for the entire unit cell is obtained by assuming a linear relationship between flux density and hydraulic head gradient and weighting the individual contributions of each flow regime over the associated liquid-occupied cross-sectional area. For solutions of the Navier–Stokes relationships for all flow regimes considered in matrix and structural pore spaces, we refer interested readers to Tuller and Or (2001b) and Or and Tuller (unpublished data, 2002).

Average Hydraulic Conductivity for Flow in Parallel Plates, Ducts, Films, and Corners
Linear relationships between steady-state flux densities obtained from the Navier–Stokes solutions and hydraulic gradient (similar to Darcy's Law) are used to derive permeabilities for parallel plates, ducts, isosceles triangles, films, and corners. Similar analogy of pore-scale flow processes with Darcy's Law was instrumental to the computation of unsaturated hydraulic conductivity in studies by Childs and Collis-George (1950), Fatt and Dykstra (1951), Burdine (1953), Gardner (1958), Mualem (1976a), and others. Darcy's Law is given as:

[19]

where Q is the volumetric discharge rate, A is the cross-sectional area occupied by the liquid, K is the hydraulic conductivity, and g is the acceleration of gravity. Assuming a unit pressure head gradient, and equating the right-hand side of Eq. [19] with Navier–Stokes solutions for average flow velocity, we receive expressions for hydraulic conductivity K for each individual flow regime considered. Closed-form solutions of the Navier–Stokes relationships for average flow velocity () between parallel plates and in full ducts are provided in many fluid mechanics textbooks (e.g., Spurk, 1997). Hydraulic conductivity expressions resulting from equating with Darcy's Law are given as:

[20]

[21]

[22]

[23]

where L is the length of the central pore in the matrix unit cell as depicted in Fig. 2, ₀ is the viscosity of bulk liquid, is the dimensionless slit-spacing parameter, B_S is a dimensionless shape parameter for square ducts (Tuller and Or, 2001b), and A_n is the area factor defined in Eq. [14].

To obtain average flow velocities in corners bounded by a liquid–vapor interface, and flow in isosceles triangular ducts with a fluid–fluid boundary (Fig. 7 and 8), we apply a finite difference–based numerical scheme to solve the Navier–Stokes relationships (Or and Tuller, unpublished data, 2002; Ransohoff and Radke, 1988). The resulting conductivity functions are given as:

[24]

[25]

where r(µ) is the capillary radius of interface curvature defined in Eq. [2], B is the fracture element aperture size, is the groove depth scaling parameter (Fig. 3), and and are dimensionless flow resistance parameters obtained from parameterization of the numerical Navier–Stokes solutions (Or and Tuller, unpublished data, 2002; Ransohoff and Radke, 1988).

On the basis of experimental and theoretical evidence showing the presence of a thin liquid layer with modified viscosity close to the solid surface, we consider two scenarios for flow in adsorbed liquid films (Fig. 8). For flow in very thin films, we substitute an exponential viscosity function relating viscosity to distance from the solid surface (Low, 1979; Or and Wraith, 1999) into the Navier–Stokes equations for plane flow, and for thicker films we assume that surface forces do not significantly alter the viscosity of bulk liquid, and therefore solve the plane flow equation with constant liquid viscosity. The resulting conductivity expressions are given as:

[26]

[27]

where h(µ) is the film thickness given in Eq. [1], and B(µ) is a constant dependent on film thickness and potential. Calculations for film flow on rough fracture surfaces (Or and Tuller, 2000) compare reasonably well with experimental data reported by Tokunaga and Wan (1997).

Note that Eq. [20] through [27] contain a fluidity term (g/₀) that could be factored out to obtain expressions for permeability (with dimensions of length squared, except Eq. [27], which involves nonlinear interactions).

Saturated and Unsaturated Hydraulic Conductivity for Unit Matrix and Fracture Elements
Saturated and unsaturated hydraulic conductivities for the entire unit elements are calculated by weighting the conductivities of single elements discussed above over the liquid occupied cross-sectional areas and dividing by the total cross-sectional element area (pore plus solid shell). Note that a key assumption here is that flow into the 2-D cross-sectional area is taking place in concurrent and parallel pathways.

The total cross-sectional area for a matrix unit cell A_MT is simply calculated by dividing the pore cross-sectional area by the porosity _M of the matrix pore space:

[28]

where A_n is the area factor of the central pore (Eq. [14]).

To calculate the unsaturated hydraulic conductivity K(µ) we consider the spontaneous redistribution (snap-off) of liquid within the slits and the central pore for drainage and imbibition (see above). The saturated hydraulic conductivity (completely filled unit matrix cell, µ µ_M-2) is given as:

[29]

where K_S and K_D are the parallel plate and duct hydraulic conductivities defined above. Unsaturated hydraulic conductivity for potentials µ_M-2 > µ µ_M-1 (partially saturated central pore and liquid-filled slits) is given as:

[30]

with K_F(µ) and K_C(µ) as the hydraulic conductivities for films and corners (Eq. [26], [27], and [24]). Note that the thin film solution (Eq. [27]) is employed for films thinner than 10 nm. Unsaturated conductivity for potentials µ < µ_M-1 (partially saturated slits and pores) is given as:

[31]

The total cross-sectional fracture element area A_FT is obtained by dividing the fracture area by the porosity of the fracture domain _F:

[32]

The resulting expressions for hydraulic conductivity corresponding to various unit element-filling stages are given as:

[33]

Separated Interfaces—Capillary Menisci anchored at pit edges: µ_F-2 > µ µ_F-1

[34]

Separated Interfaces—Capillary menisci tangent to pit surfaces: µ < µ_F-1

[35]

where k_s, k_dI, K_F(µ), and K_C(µ) are the average hydraulic conductivities given in Eq. [20] through [27], and is a groove connectivity factor with values ranging from zero to one (Or and Tuller, 2000) that accounts for partial connectivity among neighboring pits or grooves in the direction of flow. The parameter ensures that isolated pits are not considered as part of the k_dI and K_C(µ) contributions in Eq. [33] to [35].

For a fracture geometry that satisfies the inequality in Eq. [11], we use Eq. [33] for µ µ_F-2, and Eq. [35] for all matric potentials µ < µ_F-2. Note that employing K_C(µ) (Eq. [24]) in Eq. [34] leads to a slight underestimation of K_F2(µ) within the narrow matric potential range from µ_F-2 to µ_F-1 that may be neglected for all practical purposes.

Relative hydraulic conductivity is obtained by dividing Eq. [29] through [31] and [33] through [35] by the saturated hydraulic conductivities given in Eq. [29] and [33], respectively.

Upscaling to Sample Scale—Statistical Pore-Size Distribution
To express liquid saturation and unsaturated hydraulic conductivity for a sample of structured porous medium, we consider a sample cross section made of a population of unit matrix and structural element sizes, represented by two statistical probability distributions showing distinct disparity in pore sizes between the two domains (Fig. 9).

View larger version (97K): [in this window] [in a new window]   Fig. 9. Images of (a) fractured rock, (b) macroporous and, (c) aggregated soils with conceptual sketches for dual continuum pore space representation. Note the pore size disparity between the two domains.

In the following, we illustrate derivations for matrix pores and fractures as depicted in Fig. 9a. Or and Tuller (1999) developed a statistical upscaling scheme for the pore-scale liquid retention in matrix pores (homogeneous porous media) discussed above in Liquid Saturation Relationships for Unit Pore Elements. They represented the central pore size L of the unit cell (Fig. 2) by a gamma distribution (Rice, 1995) that resembles the positive skewness of the commonly used lognormal distribution, and at the same time facilitates the derivation of closed-form sample-scale expressions for liquid retention and interfacial area. The gamma density function for the central pore size, f(L), is dependent on two parameters, and :

[36]

The parameter is limited to integer values. Calculations involving expectations of f(L) are greatly simplified by the choice of = 2, which provides a balance between flexibility and tractability ( = 2 was used in this study). The moment-generating function of the gamma distribution (Rice, 1995) is used to obtain expressions for the mean m(L) and variance v(L) of L given as:

[37]

[38]

The range of admissible L values for the assumed gamma distribution is limited to values between L_min and L_max, representing the smallest and largest central pore sizes, respectively. An indirect proportional relationship to central pore size L is used to express slit length as ß(L_max - L). Such relationship facilitates the representation of clayey soils by unit cells with relatively small central pores and large slits (on average), whereas sandy soils would be represented by larger central pores attached to shorter slits (i.e., less internal surface area). The relationships between pore length distribution and slit spacing is expressed as L. Constraints on the values of , ß, and f(L) are imposed based on measurable medium properties such as porosity, specific surface area, and measured liquid retention characteristic.

A conceptual sketch of the proposed upscaling scheme is depicted in Fig. 10. A population of gamma-distributed triangular pores (Fig. 10a) is represented for illustrative purposes by six bins (L₁–L₆). The physical model predicts the shapes of liquid–vapor interfaces for each pore size and matric potential µ, leading to different stages of pore filling, according to matric potential and pore geometry (Fig. 10b). The fraction of pores at each of the several filling stages is determined from the statistical distribution of pore length, f(L), and is expressed as the expected value of a certain range of pore spaces to be completely or partially liquid-filled. The total sample scale saturation is calculated as the weighted sum of different pore filling stages.

View larger version (40K): [in this window] [in a new window]   Fig. 10. Definition sketch showing the applied upscaling scheme. (a) Gamma distribution of central pore lengths with = 2, and six hypothetical bins (note the inverse relationship between ß and L). (b) Cell-filling stages for a population of unit cells (represented by L1 to L6) defined at three matric potential values, µ1 to µ3 (dry to wet).

Sample-Scale Liquid Saturation within the Matrix and Fracture Domains
Application of the upscaling scheme to the pore-scale expressions derived above led to closed-form expressions for sample-scale saturation as a function of matric potential (Or and Tuller, 1999). General closed-form expressions were derived for all matrix elements with regular polygon-shaped central pores and for fracture elements considering both drainage and imbibition conditions. Note that expressions derived for matrix pores need only slight modifications (e.g., neglecting slits and employing the appropriate statistical distribution) to be applicable for the structural domain in macroporous (circular pores) and aggregated (triangular pore spaces) media.

To allow separation of liquid retained due to adsorptive and due to capillary forces in the matrix domain, we express liquid saturation as a function of matric potential as the sum of five terms related to the expected values of different unit cell filling stages:

[39]

with:

[40]

[41]

[42]

[43]

[44]

The term S_wM1(µ) (Eq. [40]) is the expected value of completely filled central pores. The expected value operation amounts to integration between the smallest central pore dimension L_min (the lower limit of integration depicted in Fig. 11), and a certain pore size denoted by L₁(µ). The upper limit L₁(µ) is determined from the radius of capillary interface curvature at the onset of central pore snap-off by rearranging Eq. [7] for drainage (or Eq. [4] when considering imbibition):

[45]

View larger version (31K): [in this window] [in a new window]   Fig. 11. (a) Limits of integration for determining the expected values of liquid saturation for pores in various stages of filling. (b) Integration limits expressed in terms of pore length L as a function of potential µ.

The second term S_wM2(µ) (Eq. [41]) describes the contribution to sample scale saturation due to liquid retained in full slits. It is implicitly assumed that slits remain full for all slit spacing smaller than L₂(µ), which defines slit spacing for the onset of spontaneous slit snap-off for a given potential µ. Expressions for L₂(µ) are derived by rearranging Eq. [3]:

[46]

The third term S_wM3(µ) accounts for adsorbed liquid films within the slits after slit snap-off, with integration taken between L₂(µ) and the largest pore length L_max. The fourth term S_wM4(µ) represents the contribution of adsorbed films within the central pore after pore snap-off, with L₁(µ) and L_max being the limits of integration. The last term S_wM5(µ) calculates the contribution (in terms of its expected value) of liquid retained in corners of the central pore due to capillary forces. The limits of integration are defined between L₁(µ), the onset of pore snap-off, and L_max. Figure 11 is a definition sketch for the limits of integration as related to different pore filling stages. Generalized analytical solutions for the integrals in Eq. [40] through [44] are derived in Tuller and Or (2001b).

Sample-scale expressions for liquid saturation in the fracture domain are derived in a similar fashion by taking expectations or integrating pore scale expressions (Eq. [16]–[17]) multiplied with the gamma distribution (Eq. [36]) over portions of the fracture population associated with the different filling stages discussed above in Equilibrium Liquid Configuration within Basic Unit Elements.

The integration limits separating the fracture population are obtained by rearranging Eq. [10] and [12] and solving for the fracture apertures that separate the fracture filling stages at a certain potential, as depicted in Fig. 12. Note that the expression for the completely saturated fraction of the fracture element population is simply obtained by integrating the gamma distribution with B_min and B₁(µ) as integration limits, since the relative saturation of a saturated fracture element is 1. Individual contributions of all population groups are finally summed up for the entire matric potential range under consideration to receive the sample-scale saturation curves.

View larger version (17K): [in this window] [in a new window]   Fig. 12. Critical aperture sizes determining expected fracture-filling stages at different matric potentials.

Sample-Scale Hydraulic Conductivities for Matrix and Fracture Domains
Sample-scale expressions for unsaturated hydraulic conductivity are derived in the same fashion. The unsaturated hydraulic conductivity function for the matrix domain is expressed as the sum of five terms subjected to identical integration limits (cell filling stages) as derived in the last section:

[48]

The separation into five terms is necessary for separating contributions due to film flow and corner and duct flow. The contributions are given as follows:

[49]

[50]

[51]

[52]

[53]

where f(L) is the gamma density function with = 2 (Eq. [36]).

The first term K_M1(µ) is the expected value of hydraulic conductivity for completely filled central pores (duct flow) integrated between the minimum pore dimension L_min and the upper limit L₁(µ) (Eq. [45]). The second term K_M2(µ) describes the contribution of full slits (parallel plate flow) with integration limits L_min and L₂(µ) (Eq. [46]). Flow in films adsorbed on slit walls is accounted for in K_M3(µ) with the integration limits L₂(µ) and the maximum pore dimension L_max. The fourth term K_M4(µ) represents the contribution of films aligning the flat parts of the central pore (between the curved corner menisci). The limits of integration are defined between L₁(µ) and L_max. Finally, K_M5(µ) calculates the contribution of corner flow within the central pore from pore snap-off, L₁(µ), to L_max. For integration limits see also Fig. 11. The contribution of film flow to the overall unsaturated hydraulic conductivity is represented by the sum of terms K_M2(µ), K_M3(µ), and K_M4(µ).

The expected value of sample-scale saturated hydraulic conductivity K_Msat is given as:

[54]

The relative hydraulic conductivity is simply the quotient of K_M(µ) and K_Msat. Closed-form expressions for Eq. [49] through [53] were derived in Tuller and Or (2001b).

To avoid excessive duplication, we provide an example for the derivation of the upscaled hydraulic conductivity for the completely saturated fraction of the fracture population that ranges from the smallest aperture B_min to B₁(µ) (see Fig. 12):

with subscript F standing for "fracture". Note that the saturated hydraulic conductivity (all fractures completely saturated) is readily calculated employing Eq. [55] by changing the upper integration limit to the maximum aperture B_max. The same procedure is applied to Eq. [34] and [35] with integration limits B₁(µ) - B₂, and B₂ - B_max, respectively. The upscaled expressions are then added to yield the composite response of the whole fracture population at a certain matric potential.

For geometrical configurations with pit scaling parameters [ cos(/2)/2tan(/2)] (Eq. [11]), the upscaled expressions can be solved analytically as shown in Or and Tuller (2002). Note that these analytical solutions cover a wide variety of different geometrical configurations. Only cases with [ < cos(/2)/2tan(/2)] require numerical evaluation of the upscaled expressions.

Composite Sample-Scale Response of Matrix and Fracture Domains
The individual contributions of the matrix and fracture domains to liquid saturation are added and weighed by the porosities of the individual domains to obtain the composite saturation curve for the FPM:

[57]

where _M and _F are the matrix and fracture porosities. The composite hydraulic conductivity curve is obtained by simple addition of individual contributions:

[58]

Note that the matrix and fracture conductivity functions are already weighted by their respective porosities (i.e., the expressions are divided by total cross-sectional areas, including associated solid shells).

	MODEL APPLICATION

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Estimation of Model Parameters
The analytical sample-scale expressions for matrix and structural pore domains contain a number of free model parameters that can be either obtained from direct measurements or via fitting to measured properties (e.g., saturated hydraulic conductivity).

Model Parameters for the Matrix Domain
The proposed matrix sample-scale model for liquid saturation and hydraulic conductivity contains four free parameters; the dimensionless slit-length parameter ß, the gamma distribution parameter , the matric potential µ_d at the onset of drainage (air-entry value) to determine the largest pore length L_max, and the distribution overlap parameter that relates the dimensionless slit-spacing parameter to the largest L_max and the mean pore length m(L) (Eq. [37]) (Or and Tuller, 1999).

The dimensionless scaling parameter is determined from the adjoint Gamma distributed slit-spacing L. The extent of overlap between central pore length (L) distribution and its adjoint distribution for slit spacing (L) (Fig. 10a) plays an important role in realistic representation of media pore space. The largest slit aperture (or the leading edge of slit-spacing distribution) is determined by a distribution overlap parameter (a fitting parameter) that relates slit-spacing ratio () to largest (L_max) and mean m(L) pore lengths according to:

[59]

According to these relationships, increasing the parameter results in a decrease in overlap and smaller slit spacing. The dimensionless slit-length parameter ß is a fitting parameter and is highly dependent on measured media specific surface area.

The potential at the onset of drainage µ_d (air-entry value) is often attributed to the largest pore size present in the porous medium, which we denote as L_max. Equation [45] is used to estimate the appropriate pore size L₁(µ_d). A small correction is introduced by calculating film thickness h(µ_d) lining pore walls according to Eq. [1]. The resulting maximum pore size is given as:

[60]

The minimum pore length (L_min) was set to an arbitrary value of 10^-9 m in this study (clearly more research is needed to better bracket this parameter). These two pore sizes define the limits of the gamma distribution. L_max (or the potential µ_d at the onset of drainage) is estimated as one of the fitting parameters.

These model parameters are estimated by fitting the sample scale expression for liquid saturation (Eq. [39]) to measured liquid retention data (drainage or imbibition branch), while imposing constraints on medium specific surface area and porosity. The objective function employed in this process assures that measured medium porosity is preserved, and the calculated medium specific surface area is within 90% of the independently measured specific surface area (Or and Tuller, 1999). The medium effective Hamaker constant A_svl for solid–liquid interactions through the intervening vapor phase, used to calculate thickness of adsorbed liquid films (Eq. [1]), is remarkably constant across different soil types. It was estimated beforehand by fitting Eq. [1] to the dry end of the measured retention curve as discussed in Or and Tuller (1999). The resulting model parameters are used to calculate the sample-scale saturation curve (Eq. [39]), to predict sample-scale liquid–vapor interfacial area (Or and Tuller, 1999) and to predict sample-scale unsaturated (Eq. [48]) and saturated (Eq. [54]) hydraulic conductivities.

Model Parameters for the Structural (Fracture) Domain
In parallel we determine parameters for sample-scale expressions for liquid saturation and unsaturated hydraulic conductivity of the structural pore space, as will be illustrated here for fractured porous media. If the aperture distribution is known from measurements, we are left with four free model parameters for the fracture domain hydraulic functions. These parameters, related to fracture geometry, are the dimensionless fracture-lengths and pit-depth scaling parameters and , the pit connectivity factor , and the pit angle . For cases where information about fracture geometry is available, the free parameters can be further reduced (e.g., average pit angle , pit depth B, and pit spacing B could be obtained from image analyses of fracture cross sections). A conceptual flow chart of the proposed parameter estimation scheme is depicted in Fig. 13.

View larger version (38K): [in this window] [in a new window]   Fig. 13. A conceptual flow chart of the parameter estimation scheme. Assumed matrix pore geometry and measured matrix liquid retention data are used as input parameters to estimate free matrix model parameters while imposing surface area and porosity constraints, resulting in liquid saturation and hydraulic conductivity. Assumed fracture geometry and measured aperture-size distribution in combination with fracture porosity and saturated fracture permeability are used to determine fracture model parameters for calculating the continuous liquid saturation curve and predicting fracture hydraulic conductivity as a function of matric potential. Note that the individual contributions of the matrix and fracture domains are superimposed and weighted by matrix and fracture porosities to receive the composite response of the fractured porous medium.

The TCwt dataset (Wang and Narasimhan, 1993) contains information about matrix porosity, saturated matrix and fracture permeabilities, van Genuchten (1980) parameters _VG and n_VG for the matrix liquid saturation–matric potential relationship (water characteristic curve), and fracture spacing and effective aperture for vertical fractures. Fracture porosity for a unit volume was calculated by dividing the effective aperture through aperture spacing. The shape of the aperture distribution is approximated with the gamma distribution (Eq. [36]) with = 2 and B_min and B_max set to values of 1 x 10^-9 m and 5 x 10^-4 m, respectively. The second gamma distribution shape parameter is calculated based on the assumption that the critical aperture size is equal to the mean of the distribution m = . Reported and derived model input parameters for TCwt are listed in Table 3.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

Hydraulic Functions for Homogeneous Porous Media
A summary of sample-scale model parameters and predicted values for K_sat and medium specific surface area for the soils under investigation is given in Table 5.

View larger version (24K): [in this window] [in a new window]   Fig. 14. (a) Calculated liquid saturation for a sandy loam. (b) Predicted relative hydraulic conductivity. (Note that 1 J kg-1 10-2 bar.)

The saturated hydraulic conductivity (Eq. [54]) is readily predicted using the estimated model parameters. The predicted K_sat value of 0.032 m d^-1 (Table 5) is reasonably close to the measured value of 0.09 m d^-1 given in Table 2.

In addition, the estimated model parameters may be used to predict the liquid–vapor interfacial area, an important medium property for interfacial exchange processes (e.g., oxygen transport), as a function of matric potential. Closed-form expressions for liquid–vapor interfacial area, previously discussed in Or and Tuller (1999), were generalized, and may be now obtained with a general equation presented in Tuller and Or (2001b).

In Fig. 15 we show results for Gilat loam and compare the proposed model with the VGM model. The relatively poor fit of the saturation curve to measured data in the intermediate potential range (Fig. 15a) is attributed to the limited flexibility of the statistical gamma distribution employed for upscaling from pore- to sample scale, and the surface area constraint imposed. For a detailed discussion on the gamma distribution flexibility issues, readers are referred to Tuller and Or (2001b). The VGM model is in good agreement with measurements at a saturation range from 0.4 to 0.75, but some difficulties are encountered in matching data points close to saturation and at the dry end. When comparing the predicted relative conductivity curves, the proposed model is in favorable agreement with the measurements, capable of capturing the tail at lower matric potentials as evidenced by the measurements (Fig. 15b). In contrast, the VGM model underestimates the hydraulic conductivity close to saturation and was unable to capture the data at the dry end of the curve. The observed deviations in K(µ) prediction at low matric potentials are attributed to the rapid and monotonously decreasing functional form of the van Genuchten–Mualem model that is based solely on capillary response. The standard VGM model does not consider film flow, which contributes to hydraulic conductance at different rates, thereby altering the slope of K(µ) relationships at lower potentials. The transition from corner- and duct-dominated flow to film-dominated flow occurs under moderately wet conditions at matric potentials around -80 to -100 J kg^-1 (-0.8 to -1.0 bars).

View larger version (19K): [in this window] [in a new window]   Fig. 15. (a) Calculated liquid saturation for Gilat loam. (b) Predicted relative hydraulic conductivity.

Hydraulic Functions for Fractured Rock
The lack of complete and definitive datasets for model testing introduces undesired degrees of freedom into the evaluation. We thus view the use of the TCwt data (Wang and Narasimhan, 1993) as an illustrative example rather than a test of the model in the strict sense. Input parameters used in the scheme illustrated in Fig. 13 are given in Table 3. Figure 16 depicts the resulting water characteristic curves for matrix and the fracture domains. (Note the fracture domain is seen on the bottom left corner of Fig. 16.) As expected, the matrix domain dominates the saturation–matric potential relationships. In contrast, the permeability function (Fig. 16) is dominated by the fracture domain at low matric potentials (close to complete saturation). The transition between fracture to matrix permeability occurs at a potential of about 50 J kg^-1, and a second transition occurs at potentials of 2000 J kg^-1, where matrix film flow provides the dominant contribution to the medium permeability. The resulting unsaturated permeability curve contains three "humps", one for each of the processes. Note that the VGM model employed for the matrix domain does not show the hump at the transition from corner and duct flow–dominated to film flow–dominated regions of the permeability curve, because of its functional form (see also the previous section). It is interesting to note that fracture film flow on separated surfaces provides only a marginal contribution to transport processes, probably because of the presence of only a few and mostly small aperture-size fractures. Finally, all model-fitting parameters are summarized in Table 6.

View larger version (20K): [in this window] [in a new window]   Fig. 16. Calculated saturation and predicted permeability curves for the Tiva Canyon welded tuff (TCwt) unit. Note the corner and film flow contributions within the matrix and fracture domains. (Note that 1 J kg-1 10-2 bar.)

View larger version (23K): [in this window] [in a new window]   Fig. 17. Permeability curves for Tiva Canyon welded tuff (TCwt) under nonequilibrium conditions.

View larger version (26K): [in this window] [in a new window]   Fig. 18. (a) Fitted liquid saturation for silt loam soil with biological macropores. (b) Predicted relative hydraulic conductivity. (Note that 1 J kg-1 10-2 bar.) (Data source: Mohanty et al., 1997.)

View larger version (25K): [in this window] [in a new window]   Fig. 19. (a) Fitted liquid saturation for aggregated loam soil. (b) Predicted relative hydraulic conductivity. (Smettem and Kirkby, 1990.)

For known structural pore-size (aperture) distribution, we envision the application of an alternative theoretical approach based on concepts of CPA from percolation theory (Ambegaokar et al., 1971; Friedman and Seaton, 1998; Banavar and Johnson, 1987). The implementation of CPA in this context is based on the following argument. Given a broad structural pore-size distribution forming a 3-D network, we begin by removing all the structural pores, and then replacing the pore segments in order of decreasing size, back to their original location (from largest to smallest size). The pore size that completes the first conductive pathway across the network is labeled as the "critical" pore size (aperture). According to CPA, all pore sizes larger than the critical size are essentially in series (all flow must pass through the critical size), and all pore sizes smaller than the critical size could be in parallel but are much less conductive, thus providing a limited contribution to the overall hydraulic conductivity. Consequently, the hydraulic conductivity of the structural pore network can be represented by the hydraulic conductivity of the critical structural pore (aperture). The critical pore size is determined by finding the cumulative fraction of pore sizes larger than the critical size (L_cr) that equals the percolation threshold of the network (pc). The percolation threshold is the minimal fraction of pore sizes that span a conductive pathway, and its value depends mainly on the dimensionality of the network (d = 1, 2, or 3) and on the coordination number, Z. For simple cubic lattices pc = 0.2488. The value of the coordination number Z is difficult to determine a priori; however, evidence suggests that for structural pore networks, Z values close to 3 are common (Doyen, 1988; Sahimi, 1995; Perret et al., 1999). Hence, for Z = 3 in a 3-D pore network, the value of pc = 0.5, and the critical pore size (aperture) is equal to the mean value of the pore-size distribution. The value of L_cr can be used to estimate the saturated hydraulic conductivity of a unit structural pore element to represent the entire structural pore domain of the medium. In summary, we propose that we constrain the estimates of structural pore-size distribution such that the calculated planar (2-D) saturated hydraulic conductivity will match the 3-D estimated from the critical pore size identified by CPA.

	SUMMARY AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
The large disparity in hydraulic behavior between structural and matrix pore domains in unsaturated SPM presents a challenge to developing constitutive hydraulic relationships that are required for modeling total system response. This study reviews recent developments in modeling equilibrium liquid–vapor interfacial configurations in various pore spaces that provide approximate boundary conditions for introduction of hydrodynamic behavior. Media pore space is represented by a bimodal distribution of the two disparate populations of matrix and structural pores. The model considers a cross section in the SPM where liquid configurations take place in the plane and steady unit gradient flow takes place in the perpendicular direction, resulting in parallel flow pathways (no lateral interactions are considered). In this review we illustrate the derivation of hydraulic conductivity for the matrix domain and show detailed derivations of hydraulic functions for fractures as an example for structural pore space.

Unlike the negligible role of structural domains in liquid retention behavior, these large voids tend to dominate hydraulic conductivity of the SPM near saturation. The composite unsaturated hydraulic conductivity function contains multiple "humps" representing at least three processes: (i) flow in structural pores (dominates under wet conditions), (ii) matrix capillary flow (controls conductivity for intermediate wetness range), and (iii) film flow (dry conditions). The biggest constraint to model application is the lack of definitive datasets.

The 2-D analysis presented in this review requires some measure of structural pore-size distribution, which is likely to be derived from image analysis, because of the limited information contained in the structural domain liquid retention. However, considering evidence suggesting that <20% of the fractures in fractured porous media are water-conducting (Chiles and de Marsily, 1993), there are some doubts on the direct use of imagery information. We are thus left with limited options with respect to reliable input data for the model, and the most robust input variable is a direct measurement of SPM saturated permeability. Critical path analysis offers an additional avenue for considerations of potential 3-D effects using the structural pore-size distribution (rather than the detailed 3-D network). The results reduce the problem to calculation of the saturated permeability of a critical unit fracture element.

The primary results of this study include (i) review of a new geometrical model for SPM pore space representation; (ii) derivation of closed-form expressions for SPM liquid saturation, considering individual contributions of matrix and fracture domains; (iii) derivation of physically based functions for prediction of SPM unsaturated hydraulic conductivity; (iv) illustration of potential effects of non equilibrium conditions between matrix and fractures on the unsaturated hydraulic conductivity function in fractured porous media; and (v) introduction of potential strategies for inclusion of 3-D network effects via measured saturated permeability or CPA. Future work will focus on detailed experimental characterization of SPM properties to be used as input parameters for testing and refining the proposed model. Additionally, intermittent flow regimes and spatially nonuniform flow patterns associated with structural pore spaces must be considered and characterized.

	APPENDIX

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Notation
Dimensionless slit-spacing parameter for matrix pore geometry (Or and Tuller, 1999)

ß Dimensionless slit-length parameter for matrix pore geometry (Or and Tuller, 1999)

Angle of matrix central pore corners and fracture surface pits and grooves (°)

Groove connectivity factor (Or and Tuller, 2000)

Dimensionless flow resistance parameter for corner flow (Or and Tuller, 2000)

₀ Viscosity of bulk liquid (kg m^-1 s^-1)

Slit-spacing distribution overlap parameter for matrix pores (Or and Tuller, 1999)

µ Matric potential (J kg^-1)

µ_d Critical matric potential at the onset of matrix pore drainage (J kg^-1) (air-entry value) (Or and Tuller, 1999)

µ_F-1 Critical matric potential for the onset of meniscus recession into surface grooves (J kg^-1)

µ_F-2 Critical matric potential for the onset of spontaneous fracture element liquid displacement (J kg^-1)

µ_M-1, µ_M-2 Critical matric potentials for the onset of spontaneous matrix slit and pore liquid displacement (J kg^-1)

Dimensionless scaling parameter for fracture element length

Gamma distribution parameter (set to = 2 in this study)

Density of the liquid (kg m^-3)

Surface tension at the liquid–vapor interface (N m^-1)

Dimensionless flow resistance parameter for flow in isosceles triangular ducts

Dimensionless pit-depth scaling parameter for surface grooves

_F Porosity of the fracture pore space

_M Porosity of the matrix pore space

Gamma distribution parameter

A_n Central matrix pore area factor

A_MT, A_FT Total cross-sectional matrix and fracture element areas (m²)

A_svl Hamaker constant for solid–vapor interactions through intervening liquid (J)

A_llv Hamaker constant for liquid–liquid interactions through intervening vapor (J)

B(µ) Thin film function to account for modified liquid viscosity close to solid surfaces (Or and Tuller, 2000)

B Fracture aperture size (m)

B₁(µ) Critical aperture size that separates completely filled and partially filled fractures within the fracture population (integration limit for sample-scale expressions) (m)

B₂(µ) Critical aperture size that separates fractures with filled and partially filled surface grooves within the fracture population (integration limit for sample-scale expressions) (m)

B_maxMaximum aperture (m)

B_minMinimum aperture (m)

B_S Dimensionless shape parameter for square ducts (Tuller and Or, 2001b)

C_n Drainage and imbibition radius of curvature coefficients

f(B) Fracture aperture distribution

f(L) Matrix central pore-size distribution

F Fracture element pit or groove angularity factor

F_n Angularity factor of the central matrix pore

g Acceleration of gravity (m s^-2)

h(µ) Film thickness as a function of matric potential (m)

K_C(µ) Corner hydraulic conductivity (m s^-1)

K_D Duct hydraulic conductivity (m s^-1)

k_d Term used to calculate duct hydraulic conductivity (m^-1 s^-1)

K_F(µ) Film hydraulic conductivity (m s^-1)

K_S Parallel plate hydraulic conductivity (m s^-1)

k_s Term used to calculate parallel plate hydraulic conductivity (m^-1 s^-1)

K_Fsat, K_Msat Saturated fracture and matrix element hydraulic conductivities (m s^-1)

K_F(µ), K_M(µ) Unsaturated hydraulic conductivity for a unit fracture and matrix elements (m s^-1)

K_FPM(µ) FPM unsaturated hydraulic conductivity (m s^-1)

L Dimension of the central matrix pore (m)

L₁(µ), L₂(µ) Critical central matrix pore sizes that separate pore populations in different stages of filling (integration limit for sample-scale expressions) (m)

L_max Maximal dimension of the central matrix pore (m)

L_min Minimal dimension of the central matrix pore (m)

n Number of corners of the central matrix pore

P Perimeter of the central matrix pore (m)

p Hydraulic pressure head (m)

r(µ) Radius of interface curvature as a function of matric potential (m)

S_wF, S_wM Relative liquid saturation for the fracture and matrix domains

S_FPM(µ) FPM relative liquid saturation

z Spatial coordinate along flow path (m)

	ACKNOWLEDGMENTS

 
The authors gratefully acknowledge the partial support of the NSF under grant EAR-9805409, the USDA-NRI under grant 2001-35107-10181, the Utah Agricultural Experimental Station (UAES), the Idaho Agricultural Experimental Station (IAES), and the collaboration supported by the regional technical committee W-188. Approved as IAES journal paper 2706.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORETICAL CONSIDERATIONS MODEL APPLICATION RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 

• Ambegaokar, V., N.I. Halperin, and J.S. Langer. 1971. Hopping conductivity in disordered systems. Phys. Rev. B 4:2612–2620.
• Banavar, J.R., and D.L. Johnson. 1987. Characteristic pores sizes and transport in porous media. Phys. Rev. B 35:7283–7286.
• Blunt, M., and H. Scher. 1995. Pore-level modeling of wetting. Phys. Rev. E 52:6387–6403.
• Blunt, M., D. Zhou, and D. Fenwick. 1995. Three-phase flow and gravity drainage in porous media. Transp. Porous Media 20:77–103.
• Burdine, N.T. 1953. Relative permeability calculation size distribution data. Petrol. Trans. Am. Inst. Min. Metal. Petrol. Eng. 198:71–78.
• Childs, E.C., and G.N. Collis-George. 1950. The permeability of porous materials. R. Soc. London, Proc., A 201:392–405.
• Chiles, J.P., and G. de Marsily. 1993. Stochastic models of fracture systems and their use in flow and transport modeling. p. 169–236. In J. Bear et al. (ed.) Flow and contaminant transport in fractured rock. Academic Press Inc., San Diego, CA.
• Derjaguin, B.V. 1957. Correct form of the equation of capillary condensation in porous bodies, Proc. 2nd int. Congr. Surf. Action 2:153–159.
• Derjaguin, B.V., and N.V. Churaev. 1976. Polymolecular adsorption and capillary condensation in narrow slit pores. J. Colloid Interface Sci. 54:157–175.
• Derjaguin, B.V., N.V. Churaev, and V.M. Muller. 1987. Surface forces. Plenum Publishing, New York, NY.
• Dillard, L.A., and M.J. Blunt. 2000. Development of a pore network simulation model to study nonaqueous phase liquid dissolution. Water Resour. Res. 36:439–454.
• Doyen, P.M. 1988. Permeability, conductivity and pore shape evolution of sandstone. J. Geophys. Res. 93:7729–7740.
• Dullien, F.A.L., F.S.Y. Lai, and I.F. Macdonald. 1986. Hydraulic continuity of residual wetting phase in porous media. J. Colloid Interface Sci. 109:201–218.
• Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30:211–223.
• Fatt, I., and H. Dykstra. 1951. Relative permeability studies. Petrol. Trans. AIME 192:249–255.
• Friedman, S.P., and N.A. Seaton. 1998. Critical path analysis of the relationship between permeability and electrical conductivity of three dimensional pore networks. Water Resour. Res. 34:1703–1710.
• Gardner, W.R. 1958. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85:228–232.
• Iwamatsu, M., and K. Horii. 1996. Capillary condensation and adhesion of two wetter surfaces. J. Colloid Interface Sci. 182:400–406.
• Lenormand, R., C. Zarcone, and A. Sarr. 1983. Mechanisms of the displacement of one fluid by another in a network of capillary ducts. J. Fluid Mech. 135:337–353.
• Li, Y., and N.C. Wardlaw. 1986. Mechanisms of nonwetting phase trapping during imbibition at slow rates. J. Colloid Interface Sci. 109:473–486.
• Low, P.F. 1979. Nature and properties of water in montmorillonite–water systems. Soil Sci. Soc. Am. J. 43:651–658.[Abstract/Free Full Text]
• Mason, G., and N.R. Morrow. 1986. Meniscus curvatures in capillaries of uniform cross-section. J. Chem. Soc. Faraday Trans. I. 80:2375–2393.
• Mason, G., and N.R. Morrow. 1991. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci. 141:262–274.
• Millington, R.J., and J.P. Quirk. 1961. Permeability of porous media. Nature 183:387–388.
• Mohanty, B.P. 1999. Scaling hydraulic properties of macroporous soil. Water Resour. Res. 35:1927–1931.
• Morrow, N.R., and X. Xie. 1998. Surface energy and imbibition into triangular pores. In M.Th. van Genuchten et al. (ed.) Proc. Int. Workshop on the Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. 1998. University of California, Riverside, CA.
• Mualem, Y. 1976a. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
• Mualem, Y. 1976b. A catalogue of the hydraulic properties of unsaturated soils. Research Project Report No. 442. Technion, Israel Institute of Technology, Haifa, Israel.
• Novy, R.A., P.G. Toledo, H.T. Davis, and L.E. Scriven. 1989. Capillary dispersion in porous media at low wetting phase saturations. Chem. Eng. Sci. 44:1785–1797.
• Or, D., and M. Tuller. 1999. Liquid retention and interfacial area in variably saturated porous media: Upscaling from pore to sample scale model. Water Resour. Res. 35:3591–3605.
• Or, D., and M. Tuller. 2000. Flow in unsaturated fractured porous media: Hydraulic conductivity of rough surfaces. Water Resour. Res. 36:1165–1177.
• Or, D., and J.M. Wraith. 1999. Temperature effects on soil bulk dielectric permittivity measured by time domain reflectometry: A physical model. Water Resour. Res. 35:371.
• Pachepsky, Ya.A., R.A. Shcherbakov, G. Varallyay, and K. Rajkai. 1984. On obtaining soil hydraulic conductivity curves from water retention curves. (Original in Russian.) Pochvovedenie 10:60–72.
• Perret, J., S.O. Prasher, A. Kantzas, and C. Langford. 1999. Three-dimensional quantification of macropore networks in undisturbed soil cores. Soil Sci. Soc. Am. J. 63:1530–1543.[Abstract/Free Full Text]
• Philip, J.R. 1977. Unitary approach to capillary condensation and adsorption. J. Chem. Phys. 66:5069–5075.
• Ransohoff, T.C., and C.J. Radke. 1988. Laminar flow of a wetting liquid along the corners of a predominantly gas-occupied noncircular pore. J. Colloid Interface Sci. 121:392–401.
• Restagno, F., L. Bocquet, and T. Biben. 2000. Metastability and nucleation in capillary condensation. Phys. Rev. Lett. 84:2433–2436.[Medline]
• Rice, J.A. 1995. Mathematical statistics and data analysis. 2nd ed. Duxbury Press, Belmont, CA.
• Roberts, J.N., and L.M. Schwartz. 1985. Grain consolidation and electrical conductivity in porous media. Phys. Rev. B 31:5990–5997.
• Romero, L.A., and F.G. Yost. 1996. Flow in an open channel capillary. J. Fluid Mech. 322:109–129.
• Sahimi, M. 1995. Flow and transport in porous media and fractured rock: From classical methods to modern approaches. John Wiley & Sons, New York, NY.
• Smettem, K.R.J., and C. Kirkby. 1990. Measuring the hydraulic properties of a stable aggregated soil. J. Hydrol. 117:1–13.
• Spurk, J.H. 1997. Fluid mechanics. Springer-Verlag, Berlin, Heidelberg, Germany.
• Su, G.W., J.T. Geller, K. Pruess, and F. Wen. 1999. Experimental studies of water seepage and intermittent flow in unsaturated, rough-walled fractures. Water Resour. Res. 35:1019–1037.
• Tokunaga, T.K., and J. Wan. 1997. Water film flow along fracture surfaces of porous rock. Water Resour. Res. 33:1287–1295.
• Tuller, M., and D. Or. 2001a. Hydraulic functions for macroporous soils. p. 37–40. In ASAE Proceedings of the 2nd International Symposium on Preferential Flow. Honolulu, HI. ASAE, St. Joseph, MI.
• Tuller, M., and D. Or. 2001b. Hydraulic conductivity of variably saturated porous media—Film and corner flow in angular pore space. Water Resour. Res. 37:1257–1276.
• Tuller, M., D. Or, and L.M. Dudley. 1999. Adsorption and capillary condensation in porous media: Liquid retention and interfacial configurations in angular pores. Water Resour. Res. 35:1949–1964.
• van Genuchten, M.T. 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]
• Wang, J.S.Y., and T.N. Narasimhan. 1993. Unsaturated flow in fractured porous media. p. 325– 394. In J. Bear et al. (ed.) Flow and contaminant transport in fractured rock. Academic Press Inc., San Diego, CA.

This article has been cited by other articles:

				R. Khaleel and K. P. Saripalli An Air-Water Interfacial Area Based Variable Tortuosity Model for Unsaturated Sands Vadose Zone J., May 26, 2006; 5(2): 764 - 776. [Abstract] [Full Text] [PDF]

				E. Priesack and W. Durner Closed-Form Expression for the Multi-Modal Unsaturated Conductivity Function Vadose Zone J., January 26, 2006; 5(1): 121 - 124. [Abstract] [Full Text] [PDF]

				S. Bitterlich, W. Durner, S. C. Iden, and P. Knabner Inverse Estimation of the Unsaturated Soil Hydraulic Properties from Column Outflow Experiments Using Free-Form Parameterizations Vadose Zone J., August 1, 2004; 3(3): 971 - 981. [Abstract] [Full Text] [PDF]

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Tuller, M. Articles by Or, D. Search for Related Content PubMed Articles by Tuller, M. Articles by Or, D. GeoRef GeoRef Citation Agricola Articles by Tuller, M. Articles by Or, D. Related Collections Hydraulic Conductivity Structure and Properties Preferential Flow