|
|
||||||||
Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA
* Corresponding author (hhliu{at}lbl.gov)
Received 18 December 2002.
| ABSTRACT |
|---|
|
|
|---|
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 |
|---|
|
|
|---|
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 fracturematrix interaction (flow and transport between fractures and the matrix), because of computational intensity for unsaturated flow and transport (Pruess et al., 1999). The fracturematrix 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 |
|---|
|
|
|---|
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 equilibriumbased, 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 equilibriumbased 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] |
is a positive constant depending on properties of the corresponding fracture network, and the effective water saturation in connected fractures is given by
![]() | [2] |
Note that only the active fracture continuum, a portion of the total fracture continuum, contributes to flow and transport in fractures and fracturematrix 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,
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] |
(Pa-1), n, and m = 1 - 1/n are van Genuchten parameters. In the active fracture model, however, the van Genuchten capillarypressure 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
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] |
![]() | [8] |
Combining Eq. [7] and [8] yields
![]() | [9] |
In general, relative permeability (kr) is affected by
in a complicated manner for a given Se. A larger
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
value corresponds to a smaller value of fa. Because the former effect is dominant, a larger
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 fracturematrix 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 fracturematrix 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 fracturematrix interaction, flow, and transport in the fracture continuum. The
factor may be interpreted as a measure of the "activity" of connected fractures. Generally speaking, a smaller
value corresponds to a larger number of active fractures in a connected fracture network. For example,
= 0 results in fa = 1 in Eq. [1], corresponding to all connected fractures being active. On the other hand,
= 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 |
|---|
|
|
|---|
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] |
|
![]() | [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 invasionpercolation 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].
|
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] |
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
and the fractal dimension, while
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
= 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 
Equation [20] implies that the AFM parameter
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
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
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
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] |
![]() | [22] |
![]() | [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 BrooksCorey (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,
is a constant under certain conditions for porous media. Consequently, we expect that
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 |
|---|
|
|
|---|
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
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
values associated with the TSw formation, where the repository is proposed to be located. The value of the AFM parameter
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
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
indicates that 14C data are useful for validating the AFM and for constraining the
values for the TSw unit. For
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
value generally corresponds to a larger travel time (for the matrix), because of the smaller degree of matrix diffusion resulting from a smaller fracturematrix 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.
|
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
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
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
. 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
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
,
= 0.4 may roughly represent the upper limit for the actual
values. For
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.
|
= 0.2 to 0.4 approximately match the observed 14C age data for the two boreholes simultaneously. The similar range of
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
= 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 |
|---|
|
|
|---|
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 |
|---|
| REFERENCES |
|---|
|
|
|---|
This article has been cited by other articles:
![]() |
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] |
||||