VZJ sign up for etocs
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 Similar articles in ISI Web of Science
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 ISI Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
Related Collections
Right arrow Flow
Right arrow Fractured Rock
Published in Vadose Zone Journal 3:233-246 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH PAPERS

A Modeling Study of Flow Diversion and Focusing in Unsaturated Fractured Rocks

Lehua Pan*, Yu-Shu Wu and Keni Zhang

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

Received 5 March 2003.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
 RESULTS AND DISCUSSIONS
 CONCLUSION
 REFERENCES
 
This study presents a systematic approach to analyze the flow diversion and flow focusing caused by the natural flow-barrier system in the unsaturated zone (UZ) of Yucca Mountain, Nevada, under ambient steady-state flow conditions. An existing analytical solution for analyzing capillary barrier in porous media has been extended to apply to the fractured porous rock. The new analytical solutions are used to identify the critical layers and to provide the guidance for generation of a proper three-dimensional (3-D), site-scale numerical grid. A large-scale 3-D numerical model (with more than a million grid blocks) has been developed with site-specific data to analyze the major flow patterns in the mountain. Our analyses show that large-scale lateral flow could take place in the UZ under ambient conditions, as a result of capillary barriers formed at the contacts of heterogeneous rock layers. This lateral flow runs generally toward the east (in the southern part) or southeast (in the northern part), which is consistent with the dip of the layer contacts. About 90% of the total lateral flow is found to be conducted by only a few critical rock layers. Faults that penetrate these rock layers act as vertical capillary barriers that stop the lateral flow. The combined effect of horizontal and vertical capillary barriers resulted in reduced percolation flow through repository horizon in general but focused downward flow along those penetrating faults. The model results were found to be consistent with the field water saturation. The findings of this study are consistent with a previously published two-dimensional (2-D) analysis and recent published modeling results using field-observed Cl– data.

Abbreviations: CHn, Calico Hills nonwelded • ECM, effective-continuum model • PTn, Paintbrush nonwelded • TCw, Tiva Canyon welded • TSw, Topopah Spring welded • UZ, unsaturated zone • 2-D, two-dimensional • 3-D, three-dimensional


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
 RESULTS AND DISCUSSIONS
 CONCLUSION
 REFERENCES
 
DIVERSION AND FOCUSING of water flow in the unsaturated zone may occur in many natural subsurface systems. The main factor controlling such diversion and focusing is the spatial distribution of the flow barrier (e.g., the capillary barrier formed at the contact interfaces between the strata of different soils and rocks with contrasting hydrological properties). In general, both vertical and lateral heterogeneities exist in any natural layered formation. For example, a fault system, if present, could alter the spatial distribution of flow barriers by breaking the continuity of those barriers, as well as creating new barriers, through both displacement of strata and alteration of hydrological properties. In a fractured unsaturated zone, fractures of varying intensity and geometry in different strata could further complicate the formation and distribution of these barriers. As a result, lateral water movement and focusing of downward flow may develop in the unsaturated zone and lead to a complicated flow pattern, very different from what would occur in an unsaturated, layered, porous medium system.

Problems involving capillary barriers along sloping, layered unsaturated soils were first addressed by Zaslavsky and Sinai (1981a)( 1981b). In the last few decades, several quantitative analyses of lateral water flow in layered porous media have been presented (Miyazaki, 1988; Ross, 1990, 1991; Steenhuis et al., 1991; Fayer et al., 1992; Oldenburg and Pruess, 1993; Yeh et al., 1994; Stormont, 1995; Morel-Seytoux et al., 1996; Wilson, 1996; Warrick et al., 1997; Pan et al., 1997; Webb, 1997; Ho and Webb, 1998; Morel-Seytoux and Nimmo, 1999). Most investigations have focused on developing analytical or numerical approaches for analyzing flow diversion as the effect of either a single capillary barrier at the layer contact interface or multiple but parallel substrata interfaces. What is lacking in the literature is a systematic analysis of 3-D flow diversion and focusing in an unsaturated zone. In addition, understanding of the interactions between barriers that are nonparallel or intersecting remains very limited.

During the early site characterization of the Yucca Mountain UZ as a potential repository of high-level radioactive waste (Montazer and Wilson, 1984), the capillary barrier and lateral flow concept was proposed. However, in situ direct measurement of the lateral flow is technically difficult, if not impossible, mainly because of the dry conditions of Yucca Mountain. On the other hand, modeling studies, with proper conceptualization and model parameters incorporated, can (i) provide a complete picture of the complicated flow system for integrating field measurements, (ii) indicate the measurable variables to be used as indirect evidence to support or deny a particular mechanism, and (iii) guide field experiments and future modeling efforts to improve our understanding of the system. As a result, several numerical models (Rulon et al., 1986; Wittwer et al., 1995; Moyer et al., 1996; Wu et al., 1998, 1999, 2000a) have been used to explore the capillary-barrier phenomena, but they came to quite different conclusions in terms of the amount of lateral flow that would occur in fractured tuffs. With improved understanding of the site, scientists recently undertook a modeling study (Wu et al., 2002b) of Yucca Mountain using several 2-D numerical models. Their study identified some key hydrogeologic conditions required to form capillary barriers and indicated possible large-scale lateral flow in the Yucca Mountain UZ.

Our objectives are (i) to present an extension of the analytical solutions of Warrick et al. (1997) to fractured rock, using an effective-continuum model (ECM), (ii) to present a systematic analysis of the complicated 3-D flow diversion and focusing in the UZ at Yucca Mountain using site-specific data, and (iii) to compare the roles of the analytical solutions and the numerical analysis in the systematic approach. Our focus will be on flow at the ambient, steady-state conditions.


   HYDROGEOLOGICAL CONDITIONS AND CONCEPTUALIZATION
TOP
 ABSTRACT
 INTRODUCTION
 HYDROGEOLOGICAL CONDITIONS AND...
APPROACHES
 RESULTS AND DISCUSSIONS
 CONCLUSION
 REFERENCES
 
Located in the arid western United States, the thick UZ at Yucca Mountain, Nevada, is currently under consideration by the USDOE as a potential repository for storage of high-level radioactive waste. As shown in Fig. 1 , the domain of the UZ model encompasses approximately 40 km2 of the Yucca Mountain area (Hinds and Pan, 2000; Wu et al., 2000b). Vertically, the UZ is between 500 and 700 m thick; it overlies a relatively flat water table in the vicinity of the potential repository area. The potential repository would be located in the highly fractured Topopah Spring Tuff, approximately 200 m or more above the water table.



View larger version (36K):
[in this window]
[in a new window]
  		Fig. 1. Plane view of the unsaturated zone model domain, showing the model boundary, the potential repository outline, major fault locations from Geological Framework Model (GFM) 3.1, the paths of two underground tunnels—the Exploratory Studies Facility (ESF) and Enhanced Characterization of the Repository Block (ECRB), selected boreholes, and the location of cross sections used in the capillary barrier modeling studies described in this paper.

 
The net infiltration rate at the upper boundary of the UZ is observed to vary spatially and range from 0 to 15 mm yr–1, with an average value (over the domain) of about 5 mm yr–1 (Hevesi and Flint, 2000). Figure 2 shows the spatial distribution of the net infiltration. Higher infiltration rates are located primarily in the northern part of the model domain and along the mountain ridge east of the Solitario Canyon fault from south to north. The net infiltration rate is below 4 mm yr–1 for the majority of the model domain, and the low regional water recharge provides one of the necessary conditions for development of a thick unsaturated zone, as well as formation of the capillary barriers in the area.



View larger version (67K):
[in this window]
[in a new window]
  		Fig. 2. Map of net infiltration at the surface (present day, mean infiltration).

 
The UZ at Yucca Mountain consists of alternating layers of welded and nonwelded ash-flow and air-fall tuffs. Based roughly on the degree of welding, the geological formations have been divided into five major hydrogeologic units (Montazer and Wilson, 1984): the Tiva Canyon welded (TCw), the Paintbrush nonwelded (PTn), the Topopah Spring welded (TSw), the Calico Hills nonwelded (CHn), and the Crater Flat undifferentiated units. Each major unit can be further divided into subunits in the UZ flow model (Hinds and Pan, 2000), mainly on the basis of the rock matrix properties provided by Flint (1998) and other updated geological information in the current Geologic Framework Model of Yucca Mountain (Clayton, 2000). These layers of rocks generally dip to the east at about 10° or less (Fig. 3) . They also vary significantly in thickness over the model domain, and some layers are completely missing (pinched out) in certain locations. The general trend is for the rock layer to get progressively thinner as it moves to south, especially in the upper northern region (Fig. 4) . As a result, the layer contacts are generally dipping to the east or southeast (northern part), but are not necessarily parallel to the surface or to one another (Fig. 3 and 4). The dip degree can also vary significantly as shown in both cross sections. In addition to the highly heterogeneous nature of the fractured tuffs at the site, there are numerous strike-slip and normal faults with varying amounts of offset, ranging from ten to hundreds of meters (Montazer and Wilson, 1984; Day et al., 1998). These major faults generally penetrate the UZ.



View larger version (31K):
[in this window]
[in a new window]
  		Fig. 3. Geological profile along a cross section through Borehole UZ-14, taken at a Northing coordinate of 235087 m (figure does not extend to the water table).

 


View larger version (26K):
[in this window]
[in a new window]
  		Fig. 4. Geological profile along north–south cross section, taken at an Easting coordinate of 170600 m (figure does not extend to the water table).

 
Many types of field data have been collected from the Yucca Mountain site in the past two decades. These data have been used to formulate the conceptual model and to calibrate the numerical models that describe the mountain's hydrologic system. In these numerical models, the dual-permeability modeling approach is normally used to handle fracture–matrix flow in both fractured tuffs and fault zones because this conceptual model can successfully match many types of field data at Yucca Mountain (Ahlers and Liu, 2000). Calibrated hydraulic properties for both matrix and fractures (Ahlers and Liu, 2000; Bandurraga and Bodvarsson, 1999) are used in this study to build models that mimic the real system. In general, welded tuffs have tighter matrix and more fractures than nonwelded tuffs. Based on the measured hydraulic properties of the borehole samples (Flint, 1998), the variability within the same layer is much less than that between layers. Therefore, the layered structure is considered as the main conceptual model to account for the spatial heterogeneity in the UZ formation at Yucca Mountain. The rock in the same geological layer is assumed to be uniform, except for the layers that underwent secondary geological alterations (e.g., zeolitic/vitric zones in Calico Hill unit). Within zeolitic or vitric zones, proper rock properties are used accordingly (i.e., the same layer may have different rocks). Note that the thickness and the dip angle of each layer are also spatially variable over the domain. The detailed hydraulic parameters of fracture and matrix used in this study are provided in Table 1.


View this table:
[in this window]
[in a new window]
  		Table 1. Hydraulic parameters of fracture and matrix used in this study.{dagger}

 
In this paper, faults are treated in the numerical model as continuous, vertical, highly fractured rock zones (Wu et al., 2002b). They can cause discontinuity of geological layers and may serve as structural barriers to lateral flow under unsaturated conditions. Consequently, vertical capillary barriers may form along these faults and intersect possible lateral capillary barriers. On the other hand, high permeability in a fault zone could facilitate rapid vertical drainage of the water accumulated by capillary diversion (Wu et al., 2002b).


   APPROACHES
TOP
 ABSTRACT
 INTRODUCTION
 HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
RESULTS AND DISCUSSIONS
 CONCLUSION
 REFERENCES
 
Analytical Solutions
Analytical solutions, when available, are more useful for illustrating fundamental characteristics of a capillary barrier system, even though they rely on simplified idealization of complicated subsurface flow processes. We will use analytical solutions to identify rock layer contacts that may act as capillary barriers in the Yucca Mountain UZ under ambient infiltration and to evaluate the potential capillary barrier diversion length in response to various infiltration rates. Such information can also be used to guide how to design a proper grid for numerical modeling (e.g., finer resolutions in the neighborhood of these layer contacts). Because there are significant spatial variations in layer thickness and hydraulic properties between the northern part and the southern part of the Yucca Mountain formation, we select two representative layer combinations, based on data from Boreholes UZ-14 (representing the northern part) and SD-12 (representing the southern part). The profile at UZ-14 extends from the elevation of 1341.1 m (top of PTn21) down to the elevation of 667.5 m (bottom of bf3), while the profile at SD-12 extends from the elevation of 1322.1 m (top of TCw12) down to an elevation of 660.5 m (bottom of bf3). Because general dipping angles of geological layers are 10° or less in Yucca Mountain area (Hinds and Pan, 2000), a dip angle of 10° is used for both profiles to represent the upper bound of the potential capillary barrier effects.

To facilitate the analytical solutions, we use the ECM concept here to illustrate some important behavioral aspects of the flow system. The ECM assumes a local thermodynamic or capillary equilibrium condition between fracture and matrix continua. Under such a flow condition, in terms of the ECM formulation, the Richards' equation can be written in the same form as for flow in a single-continuum porous medium (Wu, 2000). Specifically, for steady-state percolation in an unsaturated zone, consisting of dipping and parallel layers, it can be written as (pseudo one-dimensional flow) (Philip, 1991):

[1]
where h, n, and ß are capillary pressure head, the normal coordinate, and the angle of the layer interface, respectively. Note that n is measured upwards and perpendicular to the interface. K(h) is effective hydraulic conductivity, as a function of capillary pressure:

[2]
where Kf,s and Km,s are the saturated hydraulic conductivity of fracture and matrix continua, respectively, while krf and krm are the relative permeability to water for fracture and matrix continua, respectively. Note that under the capillary equilibrium assumption, pressure heads in fractures and matrix are locally the same.

In the ECM formulation, the problem, as described by Eq. [1], for flow through unsaturated fractured rocks becomes equivalent to that through porous soils (Warrick et al., 1997). Therefore, the analytical solutions derived by Warrick et al. (1997) may be directly extended to the case of fractured media under the ECM approximation. The pressure profile for an unsaturated zone of p parallel layers is then described by

[3]
where Ki(h) is the effective hydraulic conductivity of the ith layer, q is the infiltration rate, ß is the dip angle, and n is the normal distance from the lowest layer contact. The terms ni and h*i are the normal distance and the pressure head at the ith layer interface, respectively, while n1 = 0 and K1 = q. The layers are counted from bottom up (i.e., the first layer is the lowest layer). Total horizontal flow Qh of a profile with p layers of rock can be calculated as (Warrick et al., 1997)

[4]
where the dimensions of Qh are meters squared per year. Note that these results do not depend on any particular form of hydraulic conductivity K(h) relationships. In the following discussion, the van Genuchten functions (van Genuchten, 1980) are used to describe both fracture and matrix continua. The effective hydraulic conductivity is calculated according to Eq. [2], based on the fracture and matrix properties calibrated with the field data (Ahlers and Liu, 2000; Bandurraga and Bodvarsson, 1999). It was found that an ECM model could be a good approximation of a dual-permeability model in simulating the capillary barrier effects in many unsaturated situations (Wu et al., 2000a).

Numerical Analysis
Numerical modeling can be used to analyze more complicated flow diversion and focusing systems than analytical approaches because a more appropriate dual-permeability model is used without the limitation of parallel, infinite rock layer assumption. Wu et al. (2002b) presented such an analysis of the detailed structures and mechanisms of the flow diversion and focusing systems with 2-D cross-sectional simulations. They found that significant capillary barrier effects might exist and result in large-scale lateral flow within the PTn unit. They also found that capillary barrier formation in unsaturated fractured rock is determined mainly by a combination of matrix–matrix and fracture–fracture flow fields.

In this study, we will focus on the overall 3-D percolation patterns in Yucca Mountain as the result of flow diversion and focusing with a full 3-D, site-scale model (more than a million grid cells). Finer vertical resolutions ({Delta}z) were used for the critical rock layers revealed in the analytical solutions. The maximum {Delta}z varies (3.5–69 m) between layers (the map view of the 3-D grid is shown in Fig. 5) . The resulting dual-permeability grid consisted of 1077522 grid cells and 4047209 connections. Hydraulic properties used was site-specific and calibrated with a site-scale numerical model (about 100000 grid cells) using field data (Ahlers and Liu, 2000; Bandurraga and Bodvarsson, 1999). Because of the large size of the grid, the numerical analyses presented in this study were performed using the parallel-computing version of the TOUGH2 unsaturated flow module (Wu et al., 2002a). The dual-permeability model was used to represent the fractured porous media. A traditional approximation of fracture–matrix flow (Warren and Root, 1963) was used, but modified using an active-fracture model (Liu et al., 1998) to incorporate fingering flow effects through fractures. The resulting discretized finite-difference equations are highly nonlinear and are solved using the Newton–Raphson iterative scheme.



View larger version (93K):
[in this window]
[in a new window]
  		Fig. 5. Plane view of the three-dimensional site-scale grid (more than 1 million grid cells).

 

   RESULTS AND DISCUSSIONS
TOP
 ABSTRACT
 INTRODUCTION
 HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
 RESULTS AND DISCUSSIONS
CONCLUSION
 REFERENCES
 
Analytical Solutions
Results obtained with the analytical solutions (Eq. [3] and [4]) showed that the capillary barrier–diverting flow of the entire profile consisted of flow from only a few critical rock layers and varied with the infiltration rate. At a lower infiltration rate, only nonwelded tuffs (with fewer fractures than welded tuffs) played an important role in laterally diverting infiltration water. Under a higher infiltration rate, however, the welded tuffs (with more fractures) started to play a more significant diversion role. The spatial distribution of flow diversion between rock layers depended significantly on the particular layer combinations and infiltration rates.

In the profile of UZ-14, the two most important rock layers (based on the magnitude of flow diversion in each layer), PTn23 and PTn21, contributed more than 99% of the total horizontal flow of 8.20 m2 yr–1, corresponding to an infiltration rate of 5 mm yr–1 (Table 2). Among them, PTn23, with a thickness of 3.9 m, conducted 88% of the total horizontal flow, and the PTn23–PTn24 contact was found to be the largest capillary barrier in the profile, which is consistent with findings by Wu et al. (2002b). The capillary barrier diversion length of the entire profile (i.e., the ratio of the total horizontal flow rate to the infiltration rate) was about 1600 m. The total horizontal flow of the entire profile was calculated as the algebraic summation over all layers, while upslope horizontal flow (i.e., negative Qh) took place in certain rock layers but with much smaller magnitudes than the down-dip flow.


View this table:
[in this window]
[in a new window]
  		Table 2. The five most critical rock layers in the UZ-14 profile.

 
Figure 6a shows a profile of the effective hydraulic conductivity corresponding to an infiltration rate of 5 mm yr–1. As shown in the figure, there is a base value for the hydraulic conductivity (i.e., the value of the infiltration rate), accompanied by many jumps in the conductivity when crossing contacts between the different layers. According to Eq. [4], because the hydraulic conductivity is a monotonic function of water-pressure head within each layer, a positive jump in the hydraulic conductivity of the profile (i.e., jumping above the base value) indicates a positive horizontal flow (i.e., down-dip flow). For a very wet case (an infiltration rate of 1000 mm yr–1), total horizontal flow of the entire profile increased to 66.0 m2 yr–1, whereas the capillary diversion length decreases to 66.0 m.



View larger version (21K):
[in this window]
[in a new window]
  		Fig. 6. Calculated effective hydraulic conductivity of the UZ-14 profile corresponding to an infiltration rate of (a) 5 mm yr–1 and (b) 1000 mm yr–1, respectively. The three most critical rock layer contacts are labeled.

 
The critical rock layers changed significantly as well. TSw31, with a thickness of 2 m and abundant fractures, becomes the most critical layer and conducts about 63% of the total horizontal flow of the entire profile (Table 2). The most critical rock layer (PTn23) for an infiltration rate of 5 mm yr–1 becomes secondary in the case of an infiltration rate of 1000 mm yr–1 and contributes less than one-half of the lateral flow from TSw31. This implies that the fractures played more important roles in flow diversion under wet conditions than dry conditions. Figure 6b shows the corresponding distribution of the effective hydraulic conductivity. Obviously, the wetter condition significantly increased the magnitude of the jumps at the TSw31–TSw32 contact and decreased those at the PTn23–PTn24 and PTn21–PTn22 layer contacts (Fig. 6a and 6b). As a result, the relative importance of these layers for the overall capillary barrier effect also changed, responding to different infiltration rates. Although the total contribution of the five most critical layers in the wet case decreased slightly (to 99.73%), horizontal flow was more evenly distributed among these critical layers in the wet case (q = 1000 mm yr–1) than in the dry case (q = 5 mm yr–1).

The capillary barrier effect in the southern part (represented by the SD-12 profile) was found to be weaker than in the northern part (represented by the UZ-14 profile) for the same dip angle and infiltration rate of 5 mm yr–1. In dry conditions, the total horizontal flow of the entire profile was 6.72 m2 yr–1, and the corresponding capillary diversion length was 1344.8 m (Table 3). However, the trend was the opposite under wet conditions. The total horizontal flow was 105.73 m2 yr–1, with a corresponding capillary diversion length of 105.73 m. In other words, the capillary diversion distances in the SD-12 profile were 82% (dry condition) and 160% (wet condition) of those in the UZ-14 profile, respectively. This is because of the difference in composition of the most critical rock layers within the two profiles. The most critical layer, PTn23, identified in the UZ-14 profile, pinched out in the SD-12 profile, whereas the zeolitic rocks of the Calico Hills formation in the UZ-14 profile turned into vitric rocks in the SD-12 profile. As a result, in dry conditions, PTn21 became the most important layer and conducted about 46% of the total horizontal flow. Notably, CH1, with a thickness of 22.5 m, and CH5, with a thickness of 14.3 m, together conducted about 54% of the total horizontal flow, whereas these layers are insignificant in conducting lateral flow in the UZ-14 profile (Table 3). This is because the zeolite content in CH1 and CH5 is much higher at UZ-14 than at SD-12, which results in very low permeability of CH1 and CH5 (including other Calico Hill units) at UZ-14. As a result, only about 50% or less of the down-dipping diversion flow takes place at elevations above the potential repository horizon (at the elevation of 1092 m) in the SD-12 profile (Table 2 and 3), in contrast to the situation in the UZ-14 profile. These south–north differences in capillary barrier composition and mechanism were not identified in the recent 2-D numerical analyses of Wu et al. (2002b) because the cross sections studied in their research were located in the northern part of the domain only. Obviously, full 3-D modeling is necessary to understand the complete picture of the flow diversion and focusing in the highly heterogeneous Yucca Mountain UZ.


View this table:
[in this window]
[in a new window]
  		Table 3. The five most critical rock layers in the SD-12 profile.

 
Numerical Analysis
Figures 7b and 7c show maps of the simulated vertical downward flux at the repository level and the water table, respectively, corresponding to the recharge (the net infiltration) at the surface (Fig. 7a). As shown in Fig. 7a, the net infiltration rate is near zero for most of the region (the red background), while the higher infiltration rates (up to three times the mean infiltration rate) are concentrated along the ridges and in the northern part. At the repository level (below PTn), this flux distribution pattern changed significantly (Fig. 7b). The background percolation rate increased slightly because of normal lateral diffusion of the moisture, especially in the northern part. Remarkably, the high-percolation areas shifted significantly laterally and tended to focus into the fault zones because of the capillary barriers and the faults in the PTn. Flow diversion and focusing continued below the repository level, a feature that can be seen by comparing Fig. 7b and 7c.



View larger version (70K):
[in this window]
[in a new window]
  		Fig. 7. Simulated downward water flux at two different elevations corresponding to the current, mean infiltration rate: (a) net infiltration at surface, (b) at repository level, and (c) at water table level. (d) Plane view of the three-dimensional grid used.

 
The mechanism causing flow diversion below the repository level in the northern part is different from that in the southern part. According to the analytical solutions (Tables 2 and 3), the vitric CHn units (e.g., CH1–CH2 or CH5–CH6) below the repository level play an important diversion role only in the southern part. The permeability of the zeolite CHn units in the northern part was too low to conduct significant diversion flow. On the other hand, the 3-D numerical model incorporated the localized low permeability in the lower TSw units and zeolitic CHn units (see Table 1) to reflect the perched water body found in many places in the northern part. Therefore, the diversion flow below the repository level in the northern part is saturated lateral flow caused by the regular permeability barrier, which is not included in the analytical solution. Consistent with the dip directions of the critical layer contacts, the infiltration flow is diverted to the east in the southern part and to the southeast in the northern part.

The roles played by vertical faults can also be seen from these figures; that is, almost all laterally diverted flow is stopped by faults. Fast vertical flow paths are now established along those faults. For those faults that are almost parallel to the diversionary flow direction (e.g., Drill Hole Wash fault, Pagany Wash fault, and Sever Wash fault), the flow-focusing zones are centered in the faults (Fig. 7b and 7c). For those faults that are almost perpendicular to the flow-diversion direction (e.g., Ghost Dance fault and Sundance fault), the flow-focusing zones are often established along the western walls of the faults (upslope side), which is consistent with the findings of Wu et al. (2002b). This difference implies that certain faults, such as the Ghost Dance fault, act more as vertical barriers than others, such as the Drill Hole Wash fault, while the latter also acts as a fast path for lateral flow (along its trail) because of its intersection angle with the diversionary flow.

Some special points need to be made regarding flow diversion and focusing in the area between the Solitario Canyon fault and the Ghost Dance fault in the southern part. The narrow, north–south trending, high-infiltration zone east of the Solitario Canyon fault on the surface (Fig. 7a) shifted to the Ghost Dance (west) fault with a shorter length (consistent with the length of the fault) at the water table (Fig. 7c), but such shifting is only partially completed at the repository level (Fig. 7b) since the critical layer PTn23 is pinched out in this area.

Figures 7b and 7c show that a medium-to-high percolation area formed to the north of the Ghost Dance (west) fault. In contrast to other areas where the high-percolation zones concentrated along the faults, many small but above-background percolation zones combined to form a region with a complicated percolation pattern. The reason for this complication can be deduced from Fig. 4. The high-permeability region corresponds to a transition zone in which the PTn becomes thinner or pinches out and the dip angle becomes flat (PTn) or is reversed (the layers below PTn). All of these features could cause a wide leaking zone within the capillary barrier. However, flow diversion to the east is still significant, as shown by the existence of high-percolation zones along the Ghost Dance fault and the Ghost Dance (west) fault. Spatial variability in the thickness and dip angle causes leakage of the capillary barriers and reduces flow diversion, but does not remove the lateral diversion flow entirely. However, the faults often completely stop the lateral diversionary flow that tends to cross them. In the relatively dry southeast area (Fig. 7a), a focused-downward flow zone is established along the east boundary (the no-flow boundary in the model) at the water table level (Fig. 7c). This zone, however, does not exist at the repository level (Fig. 7b). Instead, the focused-downward flow zone along the southern portion of the "Imbricate" fault at the repository level (Fig. 7b) virtually disappears at the water table level (Fig. 7c). This change implies that (i) large-scale diversion flow can take place underneath an area with little surface infiltration, resulting in a local high-percolation subregion, and (ii) the focused-downward flow along a fault can be diverted up to a distance of 800 m.

Note that the above modeling results indicating the flow-focusing and flow-diversion phenomena are subject to uncertainties over how well the conceptual model and the involved parameters represent the Yucca Mountain UZ (e.g., the dual-permeability model, and the layering-dominant spatial variability and steady-state flow assumptions). Because direct in situ measurement of the water flux is not feasible under such dry conditions, we present comparisons between the simulated matrix liquid saturation against the field measurements (OCRWM, DOE, 1995) at three typical boreholes (Fig. 8) . The results show that the water saturation profiles simulated by the 3-D site-scale model agree reasonably well with the field data measured by USGS, given the relatively coarse grid resolution, except near the surface where the steady-state assumption may not be valid. This implies that the 3-D numerical model used in this study is a good approximation of the moisture conditions at the Yucca Mountain UZ. In a recent 3-D flow and transport modeling study of the Yucca Mountain UZ, Wu et al. (2003) compared predicted results with field observed data at Yucca Mountain. Although both models matched field-observed hydraulic variables (e.g., water saturation and potential data in boreholes) with similar success, they found that the model with significant lateral flow in the PTn unit matched the field observed porewater Cl– concentration data better than the model without significant lateral flow. This is consistent with findings of this study (i.e., the existence of significant lateral flow in the UZ of the Yucca Mountain).



View larger version (28K):
[in this window]
[in a new window]
  		Fig. 8. Water saturation at three boreholes simulated by TOUGH2 (three-dimensional site-scale model) and measured by USGS, respectively: (a) SD-7, (b) SD-9, and (c) UZ-14.

 
Analytical Solutions vs. Numerical Analysis
A systematic modeling approach for analyzing complex flow diversion and focusing processes may benefit from the use of both analytical and numerical solutions. In the above analyses we used analytical solutions to identify critical layers within a given hydrogeological system in terms of capillary barrier effects and a numerical model to analyze in a more realistic fashion the overall processes of flow diversion and focusing in the real 3-D system. As part of such a systematic modeling approach, 2-D numerical modeling can be an important intermediate step between idealized analytical solutions and large-scale 3-D numerical modeling. As shown in our previous paper (Wu et al., 2002b), 2-D numerical modeling can be used to analyze the detailed mechanism and structures of flow diversion and focusing, especially the effects of the spatial variability of the critical layers (e.g., thickness and dip angle) and interactions between horizontal barriers (e.g., the layer contacts) and vertical barriers (e.g., the faults). While no analytical solutions exist for such systems, large-scale 3-D numerical modeling is often unable to provide such detailed information because of limitations on the grid resolution that can be afforded in practice.

Other issues need to be considered when analyzing flow in a fractured rock system. For example, analytical solutions are available only if the fractured rock can be approximated by mean of ECM media (i.e., assuming local water pressure equilibrium between the fracture and the matrix), while numerical modeling is more flexible. However, if the ECM is a proper approximation of the fractured rock, analytical solutions can be used to test the accuracy of the numerical model and to identify where finer resolutions are needed in the numerical grid. Figure 9 shows a comparison of pressure head profiles in the PTn23 layer (overlaying the PTn24 layer) as obtained using analytical and numerical solutions. Results are for a two-layer system with straight layer contact dipping at an angle of 5.6°. The pressure head is plotted vs. vertical distance from the layer contact at the middle of a 2-D profile to avoid lateral boundary effects of the numerical models (note that the analytical solution is actually a pseudo 2-D solution with infinite long layer contacts). The steady-state profiles are for a uniform filtration rate of 5 mm yr–1 at the top boundary. The comparisons show three points. First, the ECM model is as good as the dual permeability (dual K) model for the PTn units during steady-state flow at a relative low infiltration rate (5 mm yr–1). Second, finer resolutions should be given to regions near (above) the layer contact in a numerical grid to achieve more accurate solutions. Third, the numerical models are verified to be accurate by comparisons with the analytical solutions.



View larger version (26K):
[in this window]
[in a new window]
  		Fig. 9. Comparison of pressure head profiles in the PTn23 layer (overlaying the PTn24 layer), calculated using the analytical and numerical solutions.

 
Table 4 summarizes the comparisons between the analytical solutions and the 3-D numerical solutions in terms of the relative contributions of the critical layers to the total lateral flow through two representative profiles: the UZ-14 profile and the SD-12 profile. The lateral flow data were extracted from the 3-D numerical model along the UZ-14 profile (about 1890 m long) and the SD-12 profile (about 1150 m long), respectively. The UZ-14 is within the footprint of the potential repository, while the SD-12 profile is along the east boundary of the repository footprint. Similar to the analytical solutions, the 3-D numerical solutions indicated that more than 90% of the total lateral flow occurred within a few critical layers, albeit somewhat less concentrated (Table 4). However, there are also significant differences since the 3-D numerical model incorporated more realistic spatial variability in the layer thickness, the dip angle, the rock properties, and the infiltration rate—features which cannot be incorporated in the analytical models. Below we reiterate two important aspects of the use of analytical vs. numerical models.


View this table:
[in this window]
[in a new window]
  		Table 4. Comparison of simulated relative contributions to total lateral flow as obtained with analytical and numerical solutions.

 
Localized Permeability Barrier Effects
The 3-D numerical solutions indicate that the TSw37 layer conducted about 30% of total lateral flow through the UZ-14 profile while the analytical solutions show an insignificant portion in the same layer, even though the difference in the thickness is relatively small. The main reason is that the 3-D model incorporated localized rock parameters consistent with the perched water table found in that area. Perched water was found to be associated with low-permeability zeolites in the CHn or the densely welded basal vitrophyre of the TSw unit (Wu et al., 2000b). Therefore, the permeability of the layers immediately underlying the TSw37 in the Borehole UZ-14 area of the 3-D model is very low (see pcM38–pcF38 and pcM39–pcF39 in Table 1). As a result, saturated lateral flow occurred in the TSw37 layer because of permeability barriers below that layer. This mechanism differs from the capillary barrier effect. The analytical solutions in this case could not incorporate these local permeability barrier features for two reasons. One is that these low permeability barriers do not exist everywhere at the site, whereas analytical solutions can only use features that are uniformly distributed within each layer. The other is that a low permeability barrier could cause the entire profile above it to become fully saturated, thus making the analytical solutions not relevant to the UZ at Yucca Mountain. On the other hand, low-permeability barriers cause only a limited saturation (i.e., a locally perched water) in the 3-D numerical model. The same mechanism causes the lateral flow contribution of CH5 to be doubled in the SD-12 profile according to the 3-D numerical model. However, the possible local perched water body surrounding Borehole SD-12 is considered to be very small and mainly caused by the low permeability within the zeolite CH6 layer (see pcF6z–pcM6z in Table 1) in a limited area surrounding the SD-12 borehole (Wu et al., 2000b). In other words, no permeability barrier exists in the TSw units.

Thickness Variation Effects
The 3-D numerical solutions also indicate that the relative importance of the PTn21 layer within the PTn units is more significant (PTn21 = 23.55; PTn23 = 37.86) than what is simulated with the analytical solutions (PTn21 = 11.28; PTn23 = 88.32) in the UZ-14 profile (Table 4). This is mainly because an average thickness of 7.6 m for PTn21 (4.3 m for PTn23) was used in the 3-D numerical model, while a thickness of 1.9 m (4.3 m for PTn23) used for the same layer in the analytical solutions. In reality, the effects of spatial variation in the thickness of the critical layer on flow diversion could be very complicated. Local dipping angles and changes in the directions of the critical layer can lead to local "bottle necks" for lateral flow, thus affecting the direction, the magnitude, and even the existence of lateral diversion flow. Note that the analytical solutions cannot account for these complexities or heterogeneities of the faults. Even a 2-D numerical model has limited capability in this regard. Obviously, a proper 3-D numerical model is necessary to reveal flow diversion and focusing patterns when addressing real world problems. At the same time, analytical solution can be effective in identifying potential critical layers or layer contacts since they do not have any grid resolution problems, while computing requirements are nominal compared with a large scale 3-D numerical model.


   CONCLUSION
TOP
 ABSTRACT
 INTRODUCTION
 HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
 RESULTS AND DISCUSSIONS
 CONCLUSION
REFERENCES
 
A systematic analytical and 3-D numerical modeling study is presented for lateral diversion and flow focusing resulting from capillary barriers in the UZ of Yucca Mountain. An existing analytical solution for analyzing capillary barriers in porous media (Warrick et al., 1997) was extended to fractured porous rock. The new analytical solution is used to identify layers that are critical to lateral flow initiation and to provide guidance for generating a proper 3-D numerical grid. A large-scale 3-D numerical model (with more than a million grid blocks) was developed with site-specific data to analyze very complicated 3-D flow patterns in the Yucca Mountain UZ. Although the net infiltration rate at the surface depends on the topography, percolation flow through the thick unsaturated zone is controlled primarily by a few critical rock layers and faults. According to the analytical solutions (which ignore the effects of local perched water table), capillary diversion under present-day ambient conditions occurs primarily within the nonwelded units (i.e., matrix flow dominant units). Among the critical rock layers, PTn21, PTn23, and vitric CH1 are identified to conduct the most down-dip diversionary flow, while PTn22, PTn24, and vitric CH2 are layers that act as capillary barriers to downward flow. The relative importance of the critical rock layers also changes as the infiltration rate changes. The layer combination TSw31–TSw32, for example, becomes the major contributor to capillary diversion flow in the UZ-14 profile under relatively wet or high-infiltration conditions.

The 3-D numerical modeling results show that much of the downward flow is diverted toward the east (southern part) or southeast (northern part), following the dipping direction of the critical layer contacts. Most faults or fault zones act as vertical barriers to the lateral diversion flow and fast paths for the focused downward flow. In the south, focused-downward flow along faults can be diverted to a distance of as much as 800 m. The spatial distribution of the downward flux varies greatly with depth, owing to flow diversion and focusing processes caused by complicated capillary barrier effects (or permeability barrier effects in some local regions) at Yucca Mountain. As indicated by the modeling results, the percolation flux through the repository horizon may be significantly reduced since large-scale lateral flow occurs in the overlying units, which divert much of the percolation flux into narrow fault zones.

Diversionary flow is usually confined within higher effective hydraulic-conductivity layers by the presence of the lower hydraulic-conductivity layers directly above and below such layers, at the capillary pressures corresponding to the given local percolation conditions. Any spatial variability in layer thickness, dipping angle, and even hydraulic properties of the same layer may cause leakage or total failure of a capillary barrier. According to the 3-D modeling results, such a leakage area exists between the Solitario Canyon fault and the Ghost Dance fault at Yucca Mountain. However, the existence of leakage does not completely remove the capillary barrier diversion in these areas, as evidenced by high downward flow zones still forming along the Ghost Dance fault, the Ghost Dance (west) fault, and the Sundance fault.

Flow in this study is assumed to be at steady state, and the fractured rock is described by means of a dual-permeability model (in numerical simulations) or by the ECM model (in analytical solutions). Formations within the same layer are assumed to be homogeneous except for the CHn units. The effects of episodic flow, heterogeneity within a layer, and presence of the discrete fracture networks (as well as complex fracture–matrix interactions) on flow diversion and focusing in the Yucca Mountain UZ remain unclear and should be addressed in future research. However, good matches between water saturation values in the rock matrix as predicted with 3-D modeling and field measurements indicate that the approaches used in this study are capable of capturing the major flow features in the UZ at Yucca Mountain. The findings of this study are consistent with previous 2-D analyses (Wu et al., 2002b), as well as with the field-observed Cl– data (Wu et al., 2003).


   ACKNOWLEDGMENTS
 
The authors would like to thank Curtis Oldenburg and Dan Hawkes for their internal reviews and thoughtful 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
 HYDROGEOLOGICAL CONDITIONS AND...
 APPROACHES
 RESULTS AND DISCUSSIONS
 CONCLUSION
 REFERENCES
 





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 Similar articles in ISI Web of Science
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 ISI Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Pan, L.
Right arrow Articles by Zhang, K.
Related Collections
Right arrow Flow
Right arrow Fractured Rock

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome<