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 Structure and Properties Preferential Flow Hydraulic Conductivity

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 hydrau