VZJ sign up for citetrack
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Liu, H.-H.
Right arrow Articles by Bodvarsson, G. S.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Liu, H.-H.
Right arrow Articles by Bodvarsson, G. S.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Liu, H.-H.
Right arrow Articles by Bodvarsson, G. S.
Related Collections
Right arrow Vadose Zone Processes and Chemical Transport
Right arrow Unstable Flow/Fingering
Right arrow Fractured Rock
Vadose Zone Journal 2:259-269 (2003)
© 2003 Soil Science Society of America

The Active Fracture Model

Its Relation to Fractal Flow Patterns and an Evaluation Using Field Observations

Hui-Hai Liu*, Guoxiang Zhang and Gudmundur S. Bodvarsson

Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA
* Corresponding author (hhliu{at}lbl.gov)

Received 18 December 2002.



    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
The active fracture model (AFM) (Liu et al., 1998) has been widely used in modeling flow and transport in the unsaturated zone of Yucca Mountain, Nevada, a proposed repository of high-level nuclear wastes. This study presents an in-depth evaluation of the AFM, based on both theoretical arguments and field observations. We first argue that flow patterns observed from different unsaturated systems (including the unsaturated zone of Yucca Mountain) may be fractals. We derive an interesting relation between the AFM and the fractal flow behavior, indicating that the AFM essentially captures this important flow behavior at a subgrid scale. Finally, the validity of the AFM is demonstrated by the favorable comparison between simulation results based on the AFM and 14C age and fracture coating data collected from the unsaturated zone at Yucca Mountain. These data sets independently provide important insight into flow and transport processes at the Yucca Mountain site. Potential future improvements of the AFM include expanding it to consider film flow and multifractal concepts.

Abbreviations: AFM, active fracture model • CFu, mostly nonwelded and altered Crater Flat (undifferentiated) Group • CHn, mostly nonwelded and sometimes altered Calico Hills Formation • DLA, diffusion-limited aggregation • PTn, mainly nonwelded Paintbrush Group • TCw, welded Tiva Canyon Tuff • TSw, welded Topopah Spring Tuff (TSw)


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
FLOW AND TRANSPORT in unsaturated fractured rocks are generally complicated because of the complexity of fracture–matrix interaction mechanisms, distinct differences in hydraulic properties between fractures and the matrix, and nonlinearity involved in unsaturated flow. Recently, the investigation into using the unsaturated zone at Yucca Mountain, Nevada, as a storage facility for the geological disposal of high-level nuclear wastes has generated intensive research interests in modeling flow and transport in unsaturated fractured rocks (e.g., Bodvarsson and Tsang, 1999; Robinson et al., 1997; Pruess et al., 1999; Liu et al., 1998; Ho, 1997). To reliably assess the performance of the repository, it is desirable to accurately predict flow and transport processes within the mountain. Modeling flow and transport in unsaturated fractured rocks is also of interest in other areas including environmental contamination in arid and semiarid regions.

Several approaches are available in the literature for modeling flow and transport in unsaturated fractured rocks. When classified according to the manner in which fracture networks are treated in the model structure, these approaches mainly fall into one of the two categories: the continuum approach and the discrete fracture network approach. Excellent reviews on these approaches, which have been developed and used in different fields (including oil reservoir engineering, groundwater hydrology, geothermal engineering, and soil physics), can be found in Bear et al. (1993), National Research Council (1996), and Pruess et al. (1999).

In the continuum approach, fractures are considered to be sufficiently ubiquitous and distributed in such a manner that they can be meaningfully described statistically (Bear et al., 1993). The role of individual fractures in fractured media is considered to be similar to that of individual pores in porous media. Therefore, one can describe average fracture properties as macroscopic and those associated with individual fractures as microscopic. In the continuum approach, connected fractures and rock matrix are viewed as two or more overlapped interacting continua. In other words, at a "point", two or more continua are considered to coexist. In this case, the continuum mechanics formulations, such as those used for porous media, can be used to describe flow and transport in each continuum. Coupling of processes between different continua is determined by their interaction mechanisms at a subgrid scale.

The fracture network approach involves the generation, by computer simulation, of synthetic fracture networks, and the subsequent modeling of flow and transport in these networks. The approach has been extensively used for single-phase flow and transport, with deterministic, stochastic, artificial, or site-specific fracture networks (e.g., Bear et al., 1993; National Research Council, 1996), and has recently been applied to unsaturated conditions (Rasmussen, 1991; Kwicklis and Healey, 1993; Karasaki et al., 1993; Zimmerman and Bodvarsson, 1996; Liu and Bodvarsson, 2001; Liu et al., 2003a).

Both continuum and fracture network approaches have advantages and disadvantages. While the fracture network approach is useful as a tool for concept evaluation or model-based process studies, it has several limitations. First, the approach requires geometric parameters that may strongly impact flow and transport, such as fracture apertures and conductivity, but typically cannot be well constrained from field observations (Pruess et al., 1999). Second, it is difficult to separate the conductive fracture geometry from the nonconductive fracture geometry (National Research Council, 1996). Third, flow and transport models based on the approach can be complex and computationally intensive for realistic fracture densities (National Research Council, 1996). Fourth, so far the studies based on the fracture network approach have rarely considered fracture–matrix interaction (flow and transport between fractures and the matrix), because of computational intensity for unsaturated flow and transport (Pruess et al., 1999). The fracture–matrix interaction has important effects on flow and transport processes in unsaturated fractured rocks. Because the continuum approach is relatively simple and straightforward to implement, it is preferred for most applications encountered in practice (National Research Council, 1996). For example, because the number of fractures is on the order of 109 at Yucca Mountain (Doughty, 1999), it is practically impossible to construct and calibrate a discrete fracture network site-scale model with so many fractures, considering data availability and computational feasibility. Therefore, the dual-continuum approach, in which fractures and the matrix are treated as two overlapping and interactive continua, has been used as the baseline approach for modeling flow and transport within Yucca Mountain (Bodvarsson and Tsang, 1999; Liu, 2000). Note that a similar approach has also been used to model unsaturated flow and transport in structured soils (Gerke and van Genuchten, 1993).

A traditional continuum approach assumes uniformly distributed flow patterns at a subgrid scale and therefore cannot be used for representing gravity-driven fingering flow and transport in fracture networks, resulting from subsurface heterogeneities and nonlinearity involved in unsaturated flow (Glass et al., 1996; Liu et al., 1998; Su et al., 1999; Pruess et al., 1999). In an effort to incorporate this flow behavior into the continuum approach, Liu et al. (1998) developed an AFM that assumes only a portion of fractures in a connected unsaturated fracture network contribute to liquid water flow. Because the AFM can provide a large range of flow behaviors in unsaturated fractures (including fast flow behavior), it has been used extensively for modeling large-scale flow and transport in the unsaturated zone of Yucca Mountain (e.g., Moridis and Hu, 2000; Wu et al., 2000; Wu et al., 2002; Liu et al., 2003b; Zhou et al., 2003; Buscheck et al., 2002; Xu et al., 2003).

While previous studies have shown that simulation results based on the AFM are generally consistent with field observations from the unsaturated zone of Yucca Mountain (Liu et al., 1998, 2003b), a more in-depth evaluation of the AFM is highly desirable, because accurately modeling complex, large-scale flow and transport within unsaturated fractured rocks is of interest to many research areas. Therefore, the major objective of this work is to provide a further evaluation of the AFM, on the basis of both theoretical arguments and field observations. In this paper, we first show that flow patterns in unsaturated systems (including the unsaturated zone of Yucca Mountain) may be characterized by fractals. Then we will show that the AFM is approximately consistent with fractal flow patterns at the subgrid scale. Finally, we demonstrate that simulation results based on the AFM can represent different data sets collected from the Yucca Mountain site. These data sets contain important information about large-scale flow and transport processes at the site.


    ACTIVE FRACTURE MODEL
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
The AFM was developed within the context of the dual-continuum approach (Liu et al., 1998). While the details of the AFM can be found in Liu et al. (1998), a brief introduction of the AFM is provided here for convenience. The active fracture concept is based on the reasoning that, because of the fingering flow, only a portion of fractures in a connected, unsaturated fracture network contribute to liquid water flow, while other fractures are simply bypassed. The portion of the connected fractures that actively conduct water are called active fractures. We hypothesize that the number of active fractures in the unsaturated zone of Yucca Mountain is small compared with the total number of connected fractures such that active fractures, rather than total connected fractures, must be used in numerical models. We further hypothesize that the number of active fractures within a gridblock is large; hence, the continuum approach is still valid for describing fracture flow. These hypotheses are consistent with the consideration that fractures conducting water in the unsaturated zone of Yucca Mountain are many and highly dispersed (Liu, 2000).

To use the active fracture concept for modeling flow and transport in fractures, we treat active fractures as a portion of the "homogeneous" fracture continuum for a given gridblock. Note the differences between the active fracture model and the conventional, capillary equilibrium–based, fracture-water distribution model. The latter assumes that liquid water occupies first fractures with small apertures, and then fractures with relatively large apertures, as water potential (or water saturation) increases. In contrast, the active fracture model presumes gravity-dominated, nonequilibrium, preferential liquid water flow in fractures, which is expected to be similar to fingering flow in unsaturated porous media. A liquid finger can bypass a large portion of a porous medium, which does not necessarily correspond to large pores. Note that the above discussion is valid for large-scale flow processes and not inconsistent with possible validity of a capillary equilibrium–based water distribution concept at relatively small scales corresponding to an individual flow path or a single flow finger.

Flow and transport conditions and fractured rock properties should determine the fraction of active fractures in a connected fracture network, fa. An expression for fa must satisfy the following conditions: all connected fractures are active (fa = 1) if the system is fully liquid saturated, all fractures are inactive (fa = 0) if the system is at residual saturation, and fa should be related to water flux in fractures. It is generally believed that more fractures are conducive to a larger water flux. This water flux in fractures is considered to be mainly dependent on fracture saturation because fracture water flow is gravity dominated. A simple expression for fa that meets these conditions and includes only one parameter is a power function of effective water saturation in connected fractures, Se:

[1]
where {gamma} is a positive constant depending on properties of the corresponding fracture network, and the effective water saturation in connected fractures is given by

[2]
where Sf is the water saturation of all connected fractures and Sr is the residual fracture saturation. In this study, Eq. [1] is used to determine the fraction of active fractures. As will be discussed below, this simple equation is actually the first-order approximation for fractal flow behavior at a subgrid scale, although it was initially developed as an empirical relation (Liu et al., 1998).

Note that only the active fracture continuum, a portion of the total fracture continuum, contributes to flow and transport in fractures and fracture–matrix interaction. Fracture hydraulic properties should thus be defined for active fractures. The effective water saturation of active fractures, Sae, is related to the effective water saturation in connected fractures, Se, by

[3]

Because Sae <= 1, {gamma} should be smaller than or equal to one. The effective water saturation of active fractures is related to the actual water saturation in active fractures, Sa, by

[4]

If all connected fractures are considered to be active in conducting water, as assumed in previous studies, the water capillary pressure for the fracture continuum may be described by the well-known van Genuchten relation (van Genuchten, 1980):

[5]
where {alpha} (Pa-1), n, and m = 1 - 1/n are van Genuchten parameters. In the active fracture model, however, the van Genuchten capillary–pressure relation is considered to be relevant for the active fracture continuum, rather than for the whole fracture continuum. The capillary pressure for active fractures is determined by replacing Se in Eq. [5] with Sae:

[6]

Equation [6] rather than [5] should be used to simulate water flow in the fracture continuum. For a given effective water saturation in connected fractures, a larger {gamma} value corresponds to a larger effective water saturation in active fractures, and therefore to a lower absolute value for capillary pressure.

The liquid-phase relative permeability for the active fracture continuum, kar, is directly determined by the effective water saturation of active fractures. However, because only a portion of the fractures are active, the relative permeability of the entire fracture continuum, kr, should be the relative permeability of active fractures multiplied by fa, or

[7]
where kar can be given by the following van Genuchten permeability relation:

[8]

Combining Eq. [7] and [8] yields

[9]

In general, relative permeability (kr) is affected by {gamma} in a complicated manner for a given Se. A larger {gamma} value, resulting in a higher effective water saturation in active fractures (Sae), gives rise to a larger value of kar. On the other hand, a larger {gamma} value corresponds to a smaller value of fa. Because the former effect is dominant, a larger {gamma} value gives a larger relative permeability for a given effective water saturation of the fracture continuum.

Note that the van Genuchten relationships were developed for porous media rather than fracture networks. Recently, Liu and Bodvarsson (2001) developed a new constitutive-relationship model for unsaturated flow in fracture networks, based mainly on numerical experiments. They found that the van Genuchten model is approximately valid for low fracture saturations corresponding to the ambient conditions in the unsaturated zone of Yucca Mountain. Therefore, the van Genuchten relationships are still employed because it is the ambient flow process that is of interest here. However, we emphasize that the key hypothesis of the AFM is Eq. [1]. This equation can be combined with any appropriate constitutive-relationship model to simulate unsaturated flow processes in fractured porous media. Consequently, the focus of this work is on evaluation of Eq. [1].

In an unsaturated fracture network, the ratio of the interface area contributing to flow and transport between fractures and the matrix to the total interface area determined geometrically from the fracture network, is called the fracture–matrix interface area reduction factor. In the AFM, this reduction factor results from three considerations. First, the average interface area between mobile water in an active fracture and the surrounding matrix is smaller than the geometric interface area. Second, the number of active fractures is smaller than that of connected fractures. (Conventionally, all connected fractures are assumed to contribute to fracture–matrix interaction.) Third, average active fracture spacing is much larger than that for connected fractures. Under the quasi-steady-state condition, flow and transport between fractures and surrounding matrix is inversely proportional to the corresponding fracture spacing. Based on these considerations, Liu et al. (1998) derived an expression for the reduction factor:

[10]

In summary, the active fracture model uses a combination of the volume-averaged method and a simple filter to deal with fracture flow and transport. Inactive fractures are filtered out in modeling fracture–matrix interaction, flow, and transport in the fracture continuum. The {gamma} factor may be interpreted as a measure of the "activity" of connected fractures. Generally speaking, a smaller {gamma} value corresponds to a larger number of active fractures in a connected fracture network. For example, {gamma} = 0 results in fa = 1 in Eq. [1], corresponding to all connected fractures being active. On the other hand, {gamma} = 1 corresponds to zero fracture capillary pressure (Eq. [6]), indicating that all active fractures are saturated. In the latter case, the fraction of active fractures is very small for small percolation fluxes because relatively high fracture permeabilities allow most of the water to flow through only a few fractures.


    THE ACTIVE FRACTURE MODEL AND FRACTAL FLOW BEHAVIOR
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
A flow system exhibits so-called fractal flow behavior when the corresponding flow patterns can be characterized by fractals. In this section, a brief discussion of fractal dimension (used for characterizing fractal patterns) is presented, followed by a discussion of evidence for fractal flow behavior in unsaturated systems and relations between the AFM and this fractal flow behavior.

Fractal Dimension
Fractal dimension, df, is generally a noninteger and less than the corresponding Euclidean (topological) dimension of a space, D. Different kinds of definitions of fractal dimension exist (e.g., similarity dimension, Hausdorff dimension, and box dimension), although they provide very close fractal dimension values for practical applications (Feder, 1988). The most straightforward definition is the so-called box dimension, based on a simple "box-counting" procedure. This dimension is determined from Eq. [11] (below) by counting the number (N) of "boxes" (e.g., line segment, square, and cubic for one-, two-, and three-dimensional problems, respectively), needed to cover a spatial pattern, as a function of the box size (l) (e.g., Feder, 1988):

[11]
where L refers to the size of the entire spatial domain under consideration. Figure 1 shows a box-counting procedure for a spatial pattern with df = 1.6, in a two-dimensional domain with size L (Yamamoto et al., 1993).



View larger version (14K):
[in this window]
[in a new window]
 
Fig. 1. Demonstration of the "box" counting procedure for several box sizes.

 
Obviously, if a spatial pattern is uniformly distributed in space, the fractal dimension will be identical to the corresponding Euclidean dimension. In this case, the box number, N*, and the box size l have the following relation:

[12]

A fractal pattern exhibits similarity at different scales. When df < D, the corresponding pattern does not fill the whole space, but only part of it (Fig. 1).

Fractal Flow Behavior in Unsaturated Systems
Fractals have been shown to provide a common language for describing many different natural and social phenomena (Mandelbrot, 1982). While a vast literature exists on the validity of the fractal concept for a great number of fields, fractals have been found to be useful for representing many spatial distributions in subsurface hydrology, including soil particle-size distribution, roughness of fracture surface, distribution of permeability in heterogeneous formations, and large-scale solute dispersion processes (e.g., Tyler and Wheatcraft, 1990; Rieu and Sposito, 1991a, b; Perfect et al., 1996; Neuman, 1990, 1994; Molz and Boman, 1993; Molz et al., 1997). In particular, recent studies have suggested that complex flow (or solute transport) patterns in unsaturated systems can be characterized by fractals (below).

Flury and Flühler (1995) first indicated that solute leaching patterns, observed from three field plots consisting of an unsaturated loamy soil, could be well represented by a diffusion-limited aggregation (DLA) model (Witten and Sander, 1981), although the relation between DLA parameters and soil hydraulic properties is still an unresolved issue. It has been documented that DLA generates fractal patterns (e.g., Feder, 1988; Flury and Flühler, 1995). The observation of Flury and Flühler (1995) is further supported by a study of Persson et al. (2001), who used dye-infiltration data to investigate field pathways of water and solutes under unsaturated conditions. Persson et al. (2001) showed that field observations are well described by the DLA model. Furthermore, they demonstrated that observed mean power spectrum for dye penetration of a field plot displays a typical power-law relationship, another important indication of fractal flow behavior. Glass (1993) first showed that unsaturated flow in a single vertical fracture is characterized by gravity-driven fingers, and the resulting flow patterns could be modeled by an invasion–percolation approach (Wilkinson and Willemsen, 1983). Again, percolation-based models generate fractal clustering patterns (Stauffer and Aharony, 1991). Closely related to flow and transport processes in unsaturated systems, flow patterns observed from multiphase systems are also related to fractals. For example, viscous fingering in porous media has been shown experimentally to be fractal (Feder, 1988). The problem of viscous fingering in porous media is of central importance in oil recovery. Smith and Zhang (2001) also reported that DNAPL fingering in water saturated porous media, observed from sandbox experiments, is fractal.

Detailed experimental studies on unsaturated flow patterns in large-scale fracture networks are scarce in the literature because of the technical difficulties in observing these patterns in the field. Fortunately, some geochemical data sets closely related to flow patterns in large-scale fracture networks have recently become available from the unsaturated zone of Yucca Mountain. For example, a spatial distribution of fractures with mineral coatings was determined along an underground tunnel, using a detailed-line-survey method (Paces et al., 1996; Fabryka-Martin et al., 2000). The total survey length along the tunnel is several thousand meters, broken up by a number of disconnected 30-m survey intervals. As will be discussed below, fracture coating is roughly a signature of water flow paths. Therefore, this data set can be used to infer potential fractal flow behavior at a large fracture-network scale.

While the detailed line survey was performed for several geological units, we focus on the one that is densely fractured (average fracture spacing is about 0.2 m) and has the largest number of 30-m survey intervals. In this unit, about 130 coated fractures were found over the approximate 1000-m survey length. The intersections of coated fractures with the survey line (along the underground tunnel) form a set of points. If the corresponding flow pattern is fractal, the point set should be fractal too. A box-counting method is used to detect the fractal pattern from this point set. Figure 2 shows number of boxes (line segments), N, covering at least one intersection point as a function of box size (length of a segment), l. Note that the largest box size corresponds to the length of survey intervals (30 m). The smallest box size used here is 3 m. If all fractures are coated, all boxes of this size or larger should at least cover one of the intersection points because the average fracture spacing is much smaller than this size. This corresponds to a dimension of one (the topological dimension) based on Eq. [11].



View larger version (9K):
[in this window]
[in a new window]
 
Fig. 2. Relation between the number of boxes (N), covering at least one intersection between coated fractures and the survey line, as a function of box size l. The solid line corresponds to a power function (with a power of -0.5) used to fit the data points.

 
Figure 2 shows that the data points can be fitted by a power function with a power of -0.5, indicating the existence of a fractal pattern for coated fractures (or flow paths). The power value suggests a fractal dimension of 0.5 that is not an integer but smaller than the corresponding topological dimension (one). To the best of our knowledge, the fractal flow behavior in large-scale unsaturated fracture networks, observed from field data, has not been previously reported in the literature. While more studies may be needed to further confirm the findings reported here, it is generally reasonable to expect the existence of fractal flow behavior in large-scale unsaturated fracture networks, given the evidence shown here and the fact that complex unsaturated flow patterns in several different natural systems are fractals (and fractals are valid for describing many different processes).

A Relation between the AFM and Fractal Flow Patterns
We have shown that flow processes in unsaturated systems (including fracture networks) may be fractal. This has many important implications regarding the development of numerical models for large-scale unsaturated systems, because fractal patterns are generally related to fingering and/or preferential flow paths. In this section, we will demonstrate the AFM's consistency with fractal flow behavior in unsaturated fracture networks.

Consider Fig. 1a to be a gridblock containing a fracture network and the corresponding flow pattern in the fracture network to be fractal. In this case, only a portion of the medium within a gridblock contributes to water flow (Fig. 1). This is conceptually consistent with the AFM (Liu et al., 1998). Note that in Fig. 1, a box is shadowed if it covers one or more fractures (or fracture segments) that conduct water. For simplicity, further consider that fractures are randomly distributed in space, and thus the dimension for water-saturation distribution is the corresponding Euclidean dimension when all the connected fractures actively conduct water. As will be obvious from the derivations to follow, this assumption does not alter our general conclusion regarding the connection between the AFM and fractal flow patterns.

Combining Eq. [11] and [12] yields

[13]

The average water saturation (Se) for the whole gridblock (Fig. 1a) is determined to be

[14]
where V is the total water volume (excluding residual water) in fractures within the gridblock (Fig. 1a), and {phi} is fracture porosity. Similarly, the average water saturation (Sb) for shadowed boxes with a size of l is

[15]

From Fig. 1, it is obvious that there exists a box size l1 < L satisfying

[16]

Based on Eq. [13] through [16], the average saturation for shadowed boxes with size l1, Sb1, can be expressed by

[17]

Because a fractal is similar at different scales, the procedure to derive Eq. [17] from a gridblock with size L can be applied to shadowed boxes with the smaller size of l1. In this case, for a given box size smaller than l1, the number of shadowed boxes will be counted as an average number for those within the (previously shadowed) boxes with a size of l1. Again, we can find a box size l2 < l1 to obtain a saturation relation:

[18]

The procedure to obtain Eq. [18] can be continued until it reaches an iteration level, n*, at which all the shadowed boxes with a size of ln cover active fractures only. The resultant average saturation for these shadowed boxes is

[19]

By definition of active fractures, Sbn should be equivalent to the effective saturation of active fractures. It is remarkable that Eq. [19] is similar to Eq. [3], obtained from a key hypothesis of the AFM that the fraction of active fractures in an unsaturated fracture network is a power function of the average effective saturation of the network. Comparing these two equations yields

[20]

Equation [20] provides the first theoretical relation between the AFM parameter {gamma} and the fractal dimension, while {gamma} was initially developed as an empirical parameter (Liu et al., 1998). Thus, the AFM essentially captures fractal flow behavior at the gridblock scale (df < D), whereas traditional continuum approaches assume a uniform flow pattern (or effective-saturation distribution) at that scale (corresponding to df = D or {gamma} = 0). In other words, the AFM can be used for simulating fractal flow behavior (in an unsaturated fracture network) that cannot be handled by traditional continuum approaches.

Preliminary Assessment of Saturation Dependency of the AFM Parameter {gamma}
Equation [20] implies that the AFM parameter {gamma} is not a constant, as assumed by Liu et al. (1998), but a function of saturation (or other flow conditions), because both iteration level n* and df may be dependent on water saturation for a given fracture network. The relations among n*, df, and the relevant flow conditions (such as saturation) are not totally clear at this point. However, a constant {gamma} is a reasonable treatment at least for a limited range of water saturations (or flow conditions), which is the case for the unsaturated zone of Yucca Mountain where the fracture saturation is believed to be typically <10% under ambient conditions. On the other hand, experimental evidence seems to indicate that {gamma} is a weak function of saturation at least for porous media under certain conditions (which will be discussed below). It is obvious from the derivation of Eq. [20] that the AFM concept and Eq. [20] can be applied to porous media also, as long as fingering flow patterns in them are fractals. With this in mind, results from porous media may be used to conceptually evaluate the relation between {gamma} and water saturation for unsaturated fracture networks.

On the basis of laboratory experimental observations reported in the literature, Wang et al. (1998) derived a relation between flow conditions and a parameter, F, defined as the ratio of horizontal cross-sectional area occupied by gravity fingers to the total cross-sectional area. F corresponds to fa, which is defined as the portion of active fractures in a fracture network (Eq. [1]). Wang et al. (1998) related F to the ratio of average water flux through the whole cross-sectional area, q, to the saturated hydraulic conductivity of the porous medium, Ks, as follows:

[21]
for q/Ks = 0.4 - 1.0. By definition, we can relate the average water flux within fingers (qF) to q by

[22]
and also relate the average water saturation of fingers, SF, to the average water saturation for the whole cross-sectional area, Se, by

[23]

It is expected that flow within a gravitational finger is gravity dominated. In this case, we can approximately have

[24]

Equation [24] uses the Brooks–Corey (1964) model for describing the relationship between relative permeability (kr) and saturation, with ß a constant. Combining Eq. [21] to [24] yields

[25]

Comparing the above equation with Eq. [1] gives

[26]

Therefore, {gamma} is a constant under certain conditions for porous media. Consequently, we expect that {gamma} may be a weak function of saturation for unsaturated fracture networks if fingering flow patterns in a porous medium are considered to be analogous to flow patterns in the networks. However, we acknowledge that this connection needs more investigation in future studies.


    MODEL EVALUATION USING FIELD OBSERVATIONS
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
So far, we have shown that fractal flow behavior may be common in different unsaturated systems (including unsaturated fracture networks) and also that the AFM is, in general, consistent with this behavior. In this section, we will further evaluate the AFM using different data sets collected from the Yucca Mountain site. Considering the complexity of large-scale unsaturated flow and transport in fractured porous media, the direct use of actual field data for evaluating a model is critical. Previously, it has been shown that simulation results based on the AFM match the matrix saturation and water potential data collected from long boreholes at the Yucca Mountain site (e.g., Liu et al., 1998). These results are also consistent with flow and tracer transport data observed from a field test (on a scale of about 20 m) performed in a densely fractured unit within Yucca Mountain (Liu et al., 2003b). In this study, we focus on evaluating the AFM with 14C age data and fracture coating data that provide important information regarding large-scale unsaturated flow processes in fractures and fracture–matrix interaction.

Geological Setting
Geological formations of the unsaturated zone of Yucca Mountain have been grouped into stratigraphic units (on the basis of welding) by Montazer and Wilson (1984). The stratigraphic units consist of the following, in descending order from the land surface: welded Tiva Canyon Tuff (TCw), mainly nonwelded Paintbrush Group (PTn), welded Topopah Spring Tuff (TSw), mostly nonwelded and sometimes altered Calico Hills Formation (CHn), and the mostly nonwelded and altered Crater Flat (undifferentiated) Group (CFu). The nonwelded zones near the water table in the CHn and CFu can be subject to zeolitic alteration that reduces the matrix permeability by orders of magnitude (Bandurraga and Bodvarsson, 1999). Furthermore, 35 hydrogeologic layers are defined to correspond to geologic formations and coincident hydrogeologic unit boundaries.

Welded formations are highly fractured, and nonwelded formations have relatively small fracture densities. Liquid water flow occurs mainly in the tuff matrix for nonwelded formations and in fractures for welded formations. Conceptual models of flow and transport in the unsaturated zone of Yucca Mountain can be found in Bodvarsson et al. (2000), Flint et al. (2001), and Liu (2000). For each of the 35 hydrogeologic layers, hydraulic properties include permeability and van Genuchten parameters for both fracture and matrix continua, with AFM parameter {gamma} as an additional parameter. These properties are estimated by model calibration, based on small-scale property measurements, matrix saturation and potential data, and pneumatic data. The methodology of the model calibration has been reported in Bandurraga and Bodvarsson (1999), Liu et al. (1998), and Zhou et al. (2003).

Carbon-14 Data
Carbon-14 data have been collected from perched water, pore water, and gas samples from the unsaturated zone of Yucca Mountain (Yang, 2002; Fabryka-Martin et al., 2000). Water travel times from ground surface to perched water bodies were dominated by PTn, where flow occurs mainly in the rock matrix, and thus simulated water travel time to the groundwater table is insensitive to the AFM parameters. As a result, carbon-14 data collected from perched water are not used for validating the AFM. Pore-water carbon-14 data from various boreholes at Yucca Mountain may not be representative of the pore-water residence time because of probable contamination by atmospheric 14CO2 during drilling, resulting in apparently younger residence times (Yang, 2002; Fabryka-Martin et al., 2000). Carbon-14 data from gas samples are considered to be most representative of in situ conditions (Yang, 2002). Atmospheric 14CO2 during drilling is not expected to contaminate those gas samples collected several years after drilling. Gas samples were collected from two kinds of boreholes, open surface-based boreholes and instrumented (closed) surface-based boreholes. Because the data from the latter boreholes (USW SD-12 and USW UZ-1) are the most reliable indicators of in situ conditions (Fabryka-Martin et al., 2000), 14C residence ages (Fabryka-Martin et al., 2000) calculated using the data from these two boreholes are used for validating the AFM.

Gas-phase 14C ages are assumed to be representative of ages of the in situ pore water. The rationale for this assumption is provided in detail by Yang (2002). This assumption presumes rapid exchange of gas-phase CO2 (in hours to days) with dissolved CO2 and HCO-3 in pore water. Furthermore, the amount of C in an aqueous-phase reservoir relative to C in the CO2 gas-phase reservoir is a hundred times greater. Consequently, the aqueous phase will dominate the gaseous phase when exchange occurs, indicating that the assumption is reasonable (Yang 2002).

Fracture Coating Data
The process of unsaturated-zone mineral deposition is initiated during infiltration where meteoric water interacts with materials in the soil, after which a portion may then enter the bedrock fracture network. Fracture coating is generally a signature of water flow paths. Thus, the coating data are useful for validating the AFM that describes water flow in fractures.

Fracture coating data were collected in an underground tunnel in the unsaturated zone of Yucca Mountain. Observed spatial distribution of fractures with coatings is used to estimate the portion of active fractures. For a given survey interval along the tunnel, a frequency of coated fractures can be estimated for a geologic unit based on the total number of coated fractures. The ratio of the coated fracture frequency to the total fracture frequency provides an estimate of portion of the active fracture for the given geologic unit. The estimated average portion of active fractures for the TSw is about 7.2%. Note that fracture coatings may not precisely represent all the active flow paths in the unsaturated zone of Yucca Mountain (Liu et al., 1998). Nevertheless, these values give us at least an order-of-magnitude estimate of the portion of active fractures that is about 10%.

Mineral-coating growth rate data imply that the unsaturated zone fracture network has maintained a large degree of hydrologic stability with time and that fracture flow paths in the deep unsaturated zone are buffered from climate-induced variations in precipitation and infiltration (Fabryka-Martin et al., 2000). If the AFM accurately represents water flow processes, modeling results based on the AFM should be consistent with this important observation.

Comparison between Simulation Results and Field Observations
One-dimensional dual-permeability numerical models are developed for boreholes USW SD-12 and USW UZ-1. Calibrated rock properties (Liu and Ahlers, 2003) are used, except for {gamma} values associated with the TSw formation, where the repository is proposed to be located. The value of the AFM parameter {gamma} for the TSw formation is varied for different simulations to check the sensitivity of simulated water travel times to this parameter within the formation. The top boundary condition corresponds to the present-day infiltration rate for flow simulations and a constant tracer concentration for transport simulations. The initial condition for solute transport is zero concentration within the fractured rocks. Previous studies indicate that dispersion processes have an insignificant effect on overall solute transport behavior in unsaturated fractured rocks (Bodvarsson et al., 2000; Liu et al., 2003b), and therefore they are ignored here. An effective diffusion coefficient value of 1.97 x 10-10 m2 s-1 is employed, which is equal to the average value of coefficients for tritiated water measured from tuff matrix samples. TOUGH2 and T2R3D codes (Pruess, 1991; Wu et al., 1996) are used for simulating steady-state water flow and tracer transport processes. Liu et al. (2000) reported good agreement between simulation results using the above two codes for solute transport in the unsaturated zone of Yucca Mountain and those obtained using a particle tracking method without numerical dispersion, indicating that the effects of numerical dispersion are insignificant for the problem under consideration. Simulated water travel times (or ages) for rock matrix are compared with 14C ages. A simulated water travel time at a location is determined as the time when the matrix concentration reaches 50% of the top-boundary concentration. It represents the average travel time for water particles from the ground surface to the location under consideration.

Figure 3 shows simulated water travel times (ages) for different {gamma} values of the TSw formation. Note that a major path for tracer transport into the tuff matrix is from fractures to the tuff matrix, which is different from tracer transport processes within a single continuum. The considerable sensitivity of simulated results to {gamma} indicates that 14C data are useful for validating the AFM and for constraining the {gamma} values for the TSw unit. For {gamma} values ranging from 0.2 to 0.4, simulated results approximately match the observations for the two boreholes simultaneously, although a better match is obtained for USW SD-12 than that for USW UZ-14 (as a result of subsurface heterogeneity). In our current model, the heterogeneity within each geological layer is not considered. A larger {gamma} value generally corresponds to a larger travel time (for the matrix), because of the smaller degree of matrix diffusion resulting from a smaller fracture–matrix interfacial area available for mass transport between fractures and the matrix (Eq. [10]). Note that the observed 14C ages and simulated water travel times result from a combination of solute transport in fractures and the relatively slow matrix diffusion process. The spatial variability of the degree of matrix diffusion is the reason why the simulated water travel times and the observed ages in Fig. 3 are not always monotonically increasing with depth in the TSw.



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 3. Comparisons between simulated water travel times (ages) for rock matrix at boreholes (a) USW UZ-1 and (b) USW SD-12 and the corresponding 14C ages for several {gamma} values.

 
A sharp change in simulated water travel times occurs at an elevation of about 1100 m for two boreholes (Fig. 3). This is because the upper portion of the TSw unit has relatively small fracture density values and therefore corresponds to a smaller degree of matrix diffusion for a given {gamma} value. For the borehole USW UZ-1, simulated water travel time is generally longer than the observation for a given elevation. This may be owing to a subsurface heterogeneity that gives larger fracture densities (resulting in a larger degree of matrix diffusion) at the borehole location than those used in the numerical model. Layer-averaged fracture properties are used in the site-scale model of the unsaturated zone of Yucca Mountain. In general, a comparison between simulated water travel times and observed 14C ages indicates that the AFM with {gamma} values for the TSw between 0.2 to 0.4 can reasonably represent the data.

To check the consistency of the AFM with the coating data, a one-dimensional model for borehole USW SD-12 was used. The model is the same as that described above. USW SD-12 was chosen because it is located near the middle of the underground tunnel where coating data were collected. Two infiltration rates (Flint et al., 1996), present-day mean infiltration rate (3.4 mm yr-1) and glacial maximum infiltration rate (17.3 mm yr-1), were used for simulations. Again, uniform {gamma} distributions for the TSw formation were employed. The latter infiltration rate was about five times as large as the former rate and represents the maximum infiltration rate in past climates.

Figure 4 shows the simulated average portion of active fractures, fa, for the TSw formation as a function of infiltration rate and {gamma}. The average portion is calculated from Eq. [1] using the average effective saturation for the TSw formation. The calculated fa values are about 10 to 40% for {gamma} values ranging from 0.4 to 0.2, which are similar to those used for matching the 14C data. Recall that the estimate of the active fracture portion from fracture coating data in the TSw was about 10%. Note that not all the active flow paths are associated with coatings and that coatings can also disappear due to dissolution. For example, Wang et al. (1999) found a flow feature under ambient conditions from the unsaturated zone of Yucca Mountain. This flow feature does not have coatings. Therefore, coated fractures may represent the lower limit of the number of active fractures at the Yucca Mountain site. Since the number of active fractures increases with {gamma}, {gamma} = 0.4 may roughly represent the upper limit for the actual {gamma} values. For {gamma} values <0.4, the calculated fa values do not change significantly for the two infiltration rates (Fig. 4), which is consistent with the observation of the stability of flow paths over time.



View larger version (44K):
[in this window]
[in a new window]
 
Fig. 4. Simulated average portion of active fracture for the TSw formation as a function of infiltration rate and {gamma}.

 
In summary, we have shown in this section that the AFM-based simulation results for {gamma} = 0.2 to 0.4 approximately match the observed 14C age data for the two boreholes simultaneously. The similar range of {gamma} values results in fa values (for TSw) ranging from 10 to 40%, consistent with the portion of active fracture (10%) estimated from the fracture coating data. This portion (10%) is believed to represent the lower limit for fa, as discussed above. The insensitivity of fa values for {gamma} = 0.2 to 0.4 to infiltration rates is also consistent with the stability of flow paths over time, as observed in the unsaturated zone of Yucca Mountain. All these support the validity of the AFM.


    CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 
Accurately modeling flow and transport processes in natural unsaturated systems is very difficult, largely because these processes are generally localized and associated with fingering and preferential flow paths. While the continuum approach is still the most feasible choice for many large-scale problems encountered in the real world, the traditional continuum approach has generally failed in capturing localized flow behavior, because of the use of volume averaging at a subgrid scale.

Two schools of thought exist regarding how to address the problem mentioned above. A number of researchers suggest using modeling approaches that are completely different from the continuum approach, such as the DLA and percolation-based approaches (e.g., Ewing and Berkowitz, 2001; Glass 1993; Flury and Flühler 1995). These approaches have been successfully used to represent relatively small-scale observations. However, the use of these discrete approaches is very limited for large-scale applications. Furthermore, completely satisfactory theories underlying these approaches are still missing, and some critical steps in these approaches are somewhat arbitrary (Meakin and Tolman, 1989; Ewing and Berkowitz, 2001).

The other school of thought is to modify the continuum approach by incorporating key features of the discrete approaches (or small-scale observations), considering that the continuum approach is computationally efficient and robust, and often preferred for large-scale problems. Obviously, the AFM is a product of this school of thought. Note that although these discrete approaches have different physical origins, they are connected by the same class of flow patterns: fractals. Fractal flow behavior is also supported by field observations from different unsaturated systems, including the unsaturated zone of Yucca Mountain. Because of the relative simplicity of fractal-based characterizations, we believe that the key small-scale features in unsaturated systems can be successfully incorporated into large-scale continuum approaches. This is partially supported by the consistency between simulation results based on the AFM and a variety of field observations from the unsaturated zone of Yucca Mountain.

While we have demonstrated that the AFM is generally consistent with fractal flow behavior and can capture the major flow and transport features in the unsaturated zone of Yucca Mountain, more studies are needed to further improve it. For example, the AFM assumes a homogeneous distribution of flow field within the active-fracture continuum, whereas in reality the flow field may be very heterogeneous. One way to resolve this issue may be based on the concept of the multifractal (e.g., Feder 1988; Liu and Molz 1997). Another issue is that film flow is not explicitly considered in the current version of the AFM, while the film flow may be a potentially important mechanism of unsaturated flow in fractured rocks (Tokunaga and Wan, 1997). Note that the practical significance of film flow in unsaturated fractured rocks is still an issue of current debate (Tokunaga and Wan, 1997; Liu et al., 1998; Pruess et al., 1999; Or and Tuller, 2000; National Research Council, 2001). Currently, field tests in the unsaturated zone of Yucca Mountain are under way, with a major objective of further evaluating and improving the AFM. In the future, we will report on further evaluation of the AFM with data observed from these tests.


    ACKNOWLEDGMENTS
 
We are indebted to Q. Zhou and D. Hawkes at Lawrence Berkeley National Laboratory for their careful review of a preliminary version of this manuscript. We also thank two anonymous reviewers for their constructive comments. This work was supported by the Director, Office of Civilian Radioactive Waste Management, U.S. Department of Energy, through Memorandum Purchase Order EA9013MC5X between Bechtel SAIC Company, LLC, and the Ernest Orlando Lawrence Berkeley National Laboratory (Berkeley Lab). The support is provided to Berkeley Lab through the U.S. Department of Energy Contract no. DE-AC03-76SF00098.


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 ACTIVE FRACTURE MODEL
 THE ACTIVE FRACTURE MODEL...
 MODEL EVALUATION USING FIELD...
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
Vadose Zone JHome page
J. T. Birkholzer and Y. Zhang
The Impact of Fracture-Matrix Interaction on Thermal-Hydrological Conditions in Heated Fractured Rock
Vadose Zone J., May 26, 2006; 5(2): 657 - 672.
[Abstract] [Full Text] [PDF]