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 HighWire
Right arrow Citing Articles via ISI Web of Science (7)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Birkholzer, J. T.
Right arrow Articles by Tsang, Y. W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Birkholzer, J. T.
Right arrow Articles by Tsang, Y. W.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Birkholzer, J. T.
Right arrow Articles by Tsang, Y. W.
Related Collections
Right arrow Heat Transport
Right arrow Coupled Flow/Transport Models
Right arrow Fractured Rock
Published in Vadose Zone Journal 3:819-836 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

SPECIAL SECTION: RESEARCH ADVANCES IN VADOSE ZONE HYDROLOGY THROUGH SIMULATIONS WITH THE TOUGH CODES

Modeling Seepage into Heated Waste Emplacement Tunnels in Unsaturated Fractured Rock

Jens T. Birkholzer*, Sumit Mukhopadhyay and Yvonne W. Tsang

Ernest Orlando Lawrence Berkeley National Laboratory, Earth Sciences Division, 1 Cyclotron Road, MS 90-1116, Berkeley, CA 94720
* Corresponding author (jtbirkholzer{at}lbl.gov)

Received 6 October 2003.



    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 MODELING APPROACH
 MODEL RESULTS
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Predicting the amount of water that may seep into waste emplacement tunnels (drifts) is important for assessing the performance of the proposed geologic repository for high-level radioactive waste at Yucca Mountain, Nevada. The repository will be located in thick, partially saturated fractured tuff that—for the first several hundred years after emplacement—will be heated to above-boiling temperatures as a result of heat generation from the decay of radioactive waste. Heating of rock water to above-boiling conditions induces water saturation changes and perturbs water fluxes that affect the potential for water seepage into drifts. We describe numerical analyses of the coupled thermal-hydrological (TH) processes in the vicinity of waste emplacement drifts, evaluate the potential of seepage during the heating phase of the repository, and discuss the implications for the performance of the site. In addition to the capillary barrier at the rock-drift interface—independent of the thermal conditions—a second barrier exists to downward percolation at above-boiling conditions. This barrier is caused by vaporization of water in the fractured rock overlying the repository. A TOUGH2 dual-permeability simulation model was developed to analyze the combined effect of these two barriers; it accounts for all relevant TH processes in response to heating, while incorporating the capillary barrier condition at the drift wall. Model results are presented for a variety of simulation cases that cover the expected variability and uncertainty of relevant rock properties and boundary conditions.

Abbreviations: TH, thermal-hydrological


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 MODELING APPROACH
 MODEL RESULTS
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
SEEPAGE, in the context of this paper, refers to the flow of liquid water into underground openings. In geologic repositories for high-level radioactive waste, seepage defines the amount of water that can contact waste packages in the emplacement tunnels (drifts) and may potentially promote the corrosion of waste canisters. In the case of waste package failure, seepage also defines the mobilization rate of radionuclides. The more water enters the emplacement drifts, the more radionuclides can be picked up and transported into the surrounding rock. Therefore, predicting the amount of seepage into drifts is essential in assessing the long-term performance of the proposed geologic repository for high-level radioactive waste at Yucca Mountain, Nevada.

The proposed repository is to be located in the unsaturated fractured tuff of the Topopah Spring stratigraphic unit at Yucca Mountain, about 300 m below the ground surface. The unsaturated rock layers overlying and hosting the repository form a natural barrier that is expected to reduce the amount of infiltration water entering emplacement drifts. Two processes contribute to retarding and reducing seepage. The first results from the difference in the capillary strength between the formation and the drift opening, which tends to retain water in the formation pore space. As a result, water is likely to be diverted around the underground opening (e.g., Philip et al., 1989; Birkholzer et al., 1999). In the second process, the heat generated by the radioactive waste gives rise to above-boiling temperatures in the drift vicinity for several hundred years after waste emplacement. Consequently, water percolating down toward the repository is subject to vigorous vaporization in the fractured rock, which would retard and reduce the possibility of water arrival at the drifts (Ramspott, 1991). On the other hand, condensed water will form a zone of elevated water saturation above the boiling rock. Water from this zone may be mobilized to flow rapidly down toward the drift, which could promote seepage, particularly at later times when rock temperatures have decreased below boiling. Analyzing seepage during and after the period of above-boiling rock temperatures in the drift vicinity involves evaluation of the capillary barrier behavior at the rock–drift interface and the TH conditions in the fractured rock above the drift crown, including effectiveness of the vaporization barrier.

The amount of seepage into drifts at Yucca Mountain under ambient conditions was evaluated in a combined experimental and modeling analysis (e.g., Finsterle and Trautz, 2001; Trautz and Wang, 2002; Finsterle et al., 2003). Under ambient conditions, when thermal effects are not relevant, the potential for seepage depends on the capillary barrier behavior in the unsaturated fractured rock. A stochastic continuum model was developed for ambient seepage that was capable of predicting the measured seepage rates from liquid-releases tests conducted at Yucca Mountain. Key elements in this model are (i) the inclusion of heterogeneity in small-scale permeability, (ii) the use of calibrated estimates for the effective capillary strength of the fractures, and (iii) the implementation of a specific seepage boundary condition at the rock–drift interface (Finsterle et al., 2003; Ghezzehei et al., 2004). Based on the same conceptual framework, predictive modeling studies were conducted to evaluate the probability and the magnitude of seepage under ambient conditions at the Yucca Mountain site (Li and Tsang, 2003; Tsang and Li, 2003). These ambient seepage predictions are valuable to assess the long-term performance of the repository after the thermal perturbation of the near-drift flow conditions induced by the decay heat has ceased.

Under thermal conditions, seepage is affected by both the capillary barrier behavior at the drift wall and the TH behavior in the near-drift fractured rock. The TH conditions in the unsaturated fractured rock at Yucca Mountain following waste emplacement have been investigated in several modeling studies (e.g., Haukwa et al., 1999; Buscheck et al., 2002; Haukwa et al., 2003). The methodology used in these predictive studies (e.g., use of the continuum approach, dual-permeability concept) was validated against the TH response of large in situ heater tests conducted at the proposed repository site (Tsang and Birkholzer, 1999; Birkholzer and Tsang, 2000; Mukhopadhyay and Tsang, 2002; Mukhopadhyay and Tsang, 2003). While these models provide valuable insight into the expected TH conditions in the vicinity of emplacement drifts, the possibility of seepage occurring under thermally perturbed conditions, referred to hereafter as thermal seepage, was not specifically addressed.

Therefore, in this paper, we focus on the combined effect of vaporization and capillary barriers on thermal seepage. We develop a drift-scale heterogeneous process model, referred to hereafter as the thermal seepage model, that integrates the key elements of the aforementioned ambient seepage models into the conceptual framework for describing the TH conditions. This provides a unique and consistent methodology for evaluating the evolution of seepage over the entire time period that is relevant for the performance of the proposed repository, from short-term "hot" conditions to long-term "back-to-ambient" conditions, including the transition between the two regimes. Using the thermal seepage model, transient seepage rates can be directly calculated for (i) boiling conditions, when vaporization effects are dominant; (ii) for subboiling conditions, when thermally induced reflux may be relevant; and (iii) for ambient conditions, when seepage mainly depends on the capillary barrier behavior. The model is applied to several simulation cases that cover the expected range of TH conditions at Yucca Mountain. The main objectives of our model analysis can be summarized as follows:

  1. Evaluate the thermally-perturbed flow conditions in the drift vicinity and discuss impact on seepage.
  2. Analyze the duration and the effectiveness of the vaporization barrier.
  3. Analyze the effectiveness of the capillary barrier for thermally perturbed conditions.
  4. Evaluate the potential for increased seepage as a result of enhanced reflux from the condensation zone toward the drifts.
  5. Predict the evolution of seepage during and after the thermal period.


    MODELING APPROACH
 TOP
 ABSTRACT
 INTRODUCTION
 MODELING APPROACH
 MODEL RESULTS
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Modeling of the TH processes in the fractured rock was conducted using TOUGH2, a general-purpose simulator for multidimensional coupled fluid and heat flow of multiphase, multicomponent fluid mixtures in porous and fractured media (Pruess, 1987; Pruess et al., 1999). TOUGH2 accounts for the movement of gaseous and liquid phases (under pressure, viscous, capillary, and gravity forces), transport of latent and sensible heat, phase transition between liquid and vapor, and vapor pressure lowering.

Model Concept
Relevant Thermal-Hydrological Processes
The heat generated by the decaying radioactive waste gives rise to significant TH changes in the fractured rock, which last several thousand years after waste emplacement (Birkholzer et al., 2003). Figure 1 schematically illustrates these processes occurring in the unsaturated fractured rock surrounding the emplacement drifts. Initially, before the start of heating, the rock matrix in the fractured tuff is about 80 to 90% saturated with water, most of which is stagnant because of the small matrix permeability and the strong capillary forces in the matrix pores. The fractures, on the other hand, form a well-connected, highly permeable medium, but are essentially dry under ambient conditions. When heating starts after emplacement of the waste, the stagnant water in the matrix pores becomes mobile through boiling. As the pore water in the rock matrix vaporizes, the vapor moves away from the vicinity of the drift through the permeable fracture network, driven primarily by the pressure increase caused by boiling. In cooler regions away from the drift, the vapor condenses in the fractures, where it can drain either toward the heat source from above or away from the drift into the zone below the heat source. With continuous heating, a hot dryout zone will develop close to the heat source. This dryout zone forms a barrier for the condensate flowing back toward the drifts. The effectiveness of this barrier depends on the intensity of heating as well as the characteristics and the amount of reflux. Water is more likely to penetrate into the hot dryout zone if the condensate reflux occurs in the form of channelized preferential-flow, as caused by local permeability contrasts in the fracture network.



View larger version (63K):
[in this window]
[in a new window]
 
Fig. 1. Schematic of expected near-drift thermal-hydrological processes in response to repository heating (not to scale).

 
In the thermal seepage model developed in this paper, the conceptual framework for simulating the thermally induced coupled processes was adopted from previous drift-scale TH models for the Yucca Mountain (Birkholzer et al., 2003; Haukwa et al., 2003). In these models, the fractured rock is described using a dual-permeability concept, assuming two separate, but interacting continua that overlap each other in space, one describing flow and transport in the fracture network, the other describing flow and transport in the matrix (Doughty, 1999). A continuum representation of the fractures is appropriate because the fracture density is high and a well-connected fracture network forms at the scale of interest (Tsang and Birkholzer, 1999). Global flow and transport occur within both the fracture continuum and the matrix continuum, while local fracture–matrix interaction occurs between the two continua as a result of local pressure and temperature differences. Appropriate treatment of fracture–matrix interaction is important for the prediction of thermally induced water movement. The magnitude of fracture–matrix interaction defines, for example, whether condensate in the fractures will mainly be imbibed into the matrix, where it would become essentially stagnant or largely remain in the fractures, and where it could contribute to fast flow processes and enhanced seepage.

Dual-permeability concepts often assume that flow takes place through all the connected fractures and is uniformly distributed within individual fractures. In this case, the entire fracture–matrix interface area would be available for coupling of flow between the matrix and fractures, implying relatively large fracture–matrix interactions. In natural unsaturated flow systems, however, fracture flow may not be uniformly distributed because (a) flow channels may form within a fracture, and (b) only a subset of all fractures may be actively contributing to the flow processes. To account for this reduced coupling between the fracture and the matrix continua, an "active fracture model" is used that modifies the fracture–matrix interface areas for flow and transport based on the "activity" of the unsaturated connected fractures (Liu et al., 1998). The active fracture model proposes to use a fracture–matrix reduction factor proportional to a power function of liquid saturation, with the power function coefficient calibrated from measured data.

Including the small-scale heterogeneity of the fractures—the main flow conduits in the welded fractured tuff at Yucca Mountain—is essential for understanding the impact of channeled fast flow on the potential for thermal seepage. Therefore, the permeabilities of the fracture continuum are geostatistically varied to incorporate the significant heterogeneity of this parameter observed in field measurements. In this paper, the spatial variability of the fracture continuum is estimated from small-scale air-permeability testing (see Rock Properties below). These tests have demonstrated that the well-connected fracture network has permeability variations covering three to four orders of magnitude.

Capillary Barrier Behavior
Independent of thermal effects, a significant fraction of the water that percolates down the unsaturated zone at Yucca Mountain is prevented from seeping into the emplacement drifts as a result of the capillary barrier at the interface between the rock and the opening (Finsterle et al., 2003; Li and Tsang, 2003). For an emplacement drift of given shape, the seepage threshold (i.e., the critical percolation flux below which no seepage occurs) and the amount of seepage once this threshold is exceeded depend on the local flow conditions in the drift vicinity. These conditions are mainly determined by the local percolation flux reaching the drift and by the properties of the fractured rock in the vicinity of the opening. The key properties influencing seepage are (i) the capillary strength and (ii) the permeability of the fracture network near the drift wall. For flow diversion to occur, the fracture system must have sufficient connectivity and permeability to provide the necessary effective conductivity in the tangential direction around the drift. Small-scale heterogeneity increases the probability of locally breaching the capillary barrier because it promotes the existence of flow channeling and local ponding conditions close to the drift wall (Birkholzer et al., 1999). As a result, the capillary barrier fails, and water is able to drain into the opening.

The thermal seepage model incorporates the key conceptual model elements that have been developed in ambient seepage studies to accurately describe the capillary barrier behavior at the rock–drift interface (Finsterle et al., 2003; Ghezzehei et al., 2004). Most importantly, the small-scale permeability variation of the fracture continuum near the drift wall needs to be described (see Rock Properties below). Another key element of the model is the use of calibrated values for the capillary-strength parameter close to the drift wall, as estimated in comparison with measured seepage rates obtained in ambient-temperature liquid-release tests conducted at Yucca Mountain (Finsterle et al., 2003; Ghezzehei et al., 2004). The experiments were performed by sealing a short interval in a borehole located above the crown of a test tunnel and then releasing water at a constant rate into the isolated test interval. Any water that migrated from the borehole to the ceiling and dripped into the opening was captured and weighed, providing seepage as a function of time and release rate. Typically, the measured seepage rates were much smaller than the release rates, indicating that significant diversion of flow around the niches occurred. Seepage data from multiple liquid-release tests were jointly inverted to estimate the capillary-strength parameter that best described the capillary barrier and flow diversion behavior at the test location. Data from various test locations provided a range of calibrated values that represents the expected variability among the repository rocks.

A third key element is the specific seepage boundary condition implemented for the fracture continuum at the rock–drift interface. In the model, drifts are represented as open cavities with a zero capillary-strength parameter. For a vertical connection between the fracture continuum at the drift crown and the drift, the seepage flow rate qs is given by

[1]
where k is effective permeability of the unsaturated fractures, {rho} is density, µ is dynamic viscosity, Pc denotes the capillary pressure in the fractured medium, and g is gravitational acceleration. Note that the effective permeability is related to saturation and capillary pressure according to van Genuchten's characteristic functions (van Genuchten, 1980). The distance {Delta}z denotes the distance between the last formation node and the first drift node, the latter placed inside the drift immediately at the rock–drift interface. It follows that downward seepage (qs > 0) occurs only when the following condition is satisfied:

[2]
where {rho}g{Delta}z defines the threshold capillary pressure at the last node adjacent to the opening. According to Eq. [2], the fracture continuum close to the drift wall does not need to be fully saturated for seepage to commence.

Conceptually, setting a non-zero nodal distance {Delta}z at the rock–drift interface accounts for the possible presence of discrete fracture segments intersected by the drift opening. If these segments do not extend far enough laterally or do not have a lateral connection to other fractures, the water carried in these segments cannot bypass the opening. As a result, the probability for seepage would increase beyond the value expected for a porous formation, depending on whether the gravity-driven flow in the discrete fracture segment can overcome the capillary barrier. Finsterle et al. (2003) proposed using a value of {Delta}z = 0.05 m as representative for the fracture geometry observed in drifts at Yucca Mountain. For consistency, the same vertical distance is applied in this study. Using this value, the threshold capillary pressure according to Eq. [2] is about 490 Pa at ambient temperatures. Note that seepage from the rock matrix into the drift is unlikely because of its strong capillarity; thus, seepage from the matrix is not considered in the thermal seepage model.

Model Testing
The modeling framework of the thermal seepage model was tested (i) using data from heater experiments to assess model validity with respect to the coupled TH processes and (ii) using data from ambient liquid-release tests to assess model validity with respect to the capillary barrier behavior. Although details are not presented in this paper, it was shown that the model could accurately represent the measured TH response from the heater tests as well as the measured seepage rates from the liquid-release tests. While this supports the conceptual model for thermal seepage, some uncertainty remains because the thermal part and the seepage part of the model had to be tested separately. This is because the small natural percolation at the heater test sites was not favorable for seepage. As a result, no seepage was observed in these tests during and after the thermal period. The niche liquid-release tests, on the other hand, were conducted at injection rates large enough to result in seepage, but were only performed at ambient temperatures. Ideally, validation of the thermal seepage model would require a heater test operated at artificially enhanced percolation fluxes, to observe the seepage potential for extreme flow conditions and elevated temperatures. Without such testing, validation of the thermal seepage mode must be indirect, as described above.

Model Domain
The TH behavior of the fractured rock is simulated in a two-dimensional vertical cross section through the unsaturated zone at Yucca Mountain (Birkholzer et al., 2003). With the focus on near-drift conditions, it is important to represent the drift vicinity with refined discretization. Therefore, the grid spacing is smallest at the drift wall ({approx}20 cm) and increases with distance from the drift. The drift radius is 2.75 m. Since the thermally disturbed zone extends far into the overlying and underlying geological units, and since proper assignment of boundary conditions must be ensured, the vertical model domain comprises the entire unsaturated zone, using the ground surface as the upper model boundary and the groundwater table as the lower model boundary. To increase the computational efficiency of the simulations, we assumed symmetry to reduce the model domain in the lateral direction, with the symmetry plane parallel to the drift axis. The current repository design of parallel drifts spaced at 81 m can be represented as a series of symmetrical, identical half-drift models with vertical no-flow boundaries between them. Thus, the numerical mesh can be reduced to a half-drift model with a width of 40.5 m, extending from the drift center to the midpoint between drifts (see schematic cross section with model boundaries and discretization in Fig. 2) . Different cross sections (submodels) have been studied, depending on the location of the repository drifts and the geologic unit in which they reside. The half-drift model presented in this paper is representative of emplacement drifts located in the center part of the repository.



View larger version (47K):
[in this window]
[in a new window]
 
Fig. 2. Schematic cross section of Yucca Mountain showing the geometry of the model domain (not to scale). Close-up view shows discretization of drift vicinity in the thermal seepage model. The local coordinate system has its origin (x = 0, z = 0) at the center of the drift.

 
Considering that the scope of this study requires several simulation cases to account for the variability in rock properties and boundary conditions important for thermal seepage, a full three-dimensional simulation of drift-scale coupled processes was not feasible because of computational limitations. The main deviations between a three-dimensional and a two-dimensional model occur at the end of each emplacement drift and at the edges of the repository. Additional three-dimensional effects may stem from heat output variation between individual waste packages along the drift, which can lead to cooler spots where seepage may occur earlier. In part, such effects are accounted for in this study by considering various sensitivity cases for the thermal load, that is, by analyzing "cooler" and "hotter" two-dimensional models. Also note that, with respect to flow channeling and the effectiveness of the capillary barrier for seepage into drifts, a two-dimensional representation is more critical in most cases of heterogeneous fracture permeability fields because the potential diversion of flow in the third dimension is neglected.

Boundary Conditions
Both the top and bottom boundaries of the model domain are far enough away from the repository units so as not to be affected by the thermal input. (Repository units are those geologic units at Yucca Mountain in which the waste emplacement drifts will reside.) Accordingly, these boundaries are assigned Dirichlet-type conditions with fixed temperature, gas pressure, and liquid saturation values. The ground surface of the mountain at the top of the model domain is represented as an open atmosphere. The groundwater table at the bottom of the model domain is represented as a flat, stable surface saturated with water. As already mentioned above, the drift opening is represented with a zero capillary-strength parameter, and seepage between the formation and the drift is calculated using Eq. [1].

Most important for the TH behavior of the proposed repository and the potential for thermal seepage are two other boundary conditions of the model. The first is the thermal-operating mode as defined by the heat output characteristics of the radioactive waste. The thermal load of the waste packages is placed into the appropriate model elements in the open drift. The second is the surface infiltration, which is applied using a source term for the first rock gridblocks just below the top boundary. The amount and distribution of surface infiltration ultimately defines the ambient background percolation flux in the unsaturated zone of Yucca Mountain. We have conducted several simulations with alternative heat loads and surface infiltration rates to cover the expected range of variability in these boundary conditions. The heat load scenarios account for the spatial variability in the temperature conditions over the repository stemming from differences between individual waste packages, emplacement time differences among repository sections, and three-dimensional edge effects. The infiltration scenarios cover the spatial variability of surface infiltration at Yucca Mountain, and also consider the impact of long-term future climate changes.

From the suite of thermal operating modes analyzed, three thermal scenarios were selected for presentation in this paper. The "reference mode" applies an initial heat load of 1.45 kW per meter drift length, a thermal load representative of the average thermal conditions for the current repository design (Bechtel SAIC Company, 2003). The initial heat load decreases with time as a result of diminishing radioactive decay, as depicted in Fig. 3 . The figure also reflects the current design, which considers forced drift ventilation for the first 50 yr after emplacement (preclosure period). As a result, about 86% of the heat is removed during this period. The other thermal operating modes are studied as sensitivity cases. The first sensitivity case, referred to as the high-temperature mode, addresses uncertainty about the efficiency of the ventilation system. It uses the same initial heat load and decay curve as the reference mode, but assumes that the heat-removal efficiency of the ventilation system is only 70% for the first 50 yr after emplacement. The second sensitivity case, the low-temperature mode, addresses the heat output variability between individual waste packages. It considers a thermal load of 1.2 kW m–1, with the heat removal efficiency of the ventilation system similar to the reference mode. For the heat decay in this case, the time-varying values of the reference mode are linearly scaled by a factor of 1.2/1.45. Together, the three sensitivity cases cover temperature conditions in the drift vicinity that can be as high as 140°C (for the high-temperature mode) or that will exceed boiling temperature by about 12°C (for the low-temperature mode).



View larger version (11K):
[in this window]
[in a new window]
 
Fig. 3. Thermal load of the reference mode as a function of time. Time zero represents the time of emplacement of the radioactive waste.

 
Consistent with the future climate analyses for Yucca Mountain (Houseworth, 2001), the thermal seepage model considers three long-term climate periods with constant net infiltration: the present-day climate (up to 600 yr from emplacement of waste), the monsoon climate (600–2000 yr), and the glacial transition climate (more than 2000 yr). During these periods, the expected average percolation fluxes at the repository horizon are fairly small, on the order of 6, 16, and 25 mm yr–1, respectively (Birkholzer et al., 2003). These average values are used as the base case scenario in our study. Additional infiltration scenarios have been studied with the thermal seepage model to cover the potential variability of percolation fluxes within the repository footprint. The additional scenarios shown in this paper represent high-percolation regions; they are defined by multiplying the base case percolation fluxes by factors of 5, 10, and 20, respectively. Scenarios with percolation rates smaller than the expected average were also analyzed, but are not included in this paper because seepage was only observed for much larger percolation fluxes (see Prediction of Seepage at Ambient Conditions below). Also note that the thermal behavior of the low-percolation cases was almost identical to the average case. This demonstrates that percolation fluxes close to or less than the base case (i.e., a few millimeters per year for the present-day climate) are too small to affect the TH conditions in the near-drift environment. In low-percolation cases, the majority of the thermally induced vapor and condensate fluxes stems from pore water resident in the matrix rather than percolating water.

Rock Properties
Hydrological and thermal properties for the matrix and fracture continua of the various geological units in the model domain are obtained from Liu et al. (2003). Most of these layer-averaged properties are based on calibration against surface-based borehole measurements such as saturation data, water-potential data, pneumatic pressure data, ambient temperature data, and geochemical data, with the exception of the seepage-specific effective fracture capillarity of the repository unit analyzed in this paper. In agreement with the conceptual model for thermal seepage modeling, the effective fracture capillarity was chosen based on the range of calibrated parameter values from liquid-release testing (Finsterle et al., 2003; Ghezzehei et al., 2004). The fracture capillarity used in this paper is 589 Pa, representing an average value of the range reported in Finsterle et al. (2003). A summary of the rock properties for the repository unit is given in Table 1.


View this table:
[in this window]
[in a new window]
 
Table 1. Matrix and fracture properties used for repository unit.

 
In addition to the layer-averaged properties, spatial variability in fracture permeability is accounted for by assigning stochastically distributed permeability values to each fracture model element in the vicinity of the drift. The parameters needed to define the small-scale variability of the fracture continuum were based on geostatistical analyses of permeability values determined by air injection testing in boreholes drilled from test tunnels in the repository units (Wang et al., 1999). The tests were performed by isolating 0.3-m (1-ft) sections of the boreholes, using an inflatable packer system, and then injecting compressed air at a constant rate into the isolated interval. The permeability of each packed-off interval was derived from the measured pressure response once steady-state conditions were reached. Considering that the matrix is much less permeable than the fractures and that the fractures are essentially dry, these results can be considered estimates of the fracture permeability on the scale of the interval length.

The particular parameters used for the stochastic permeability fields in this study unit were derived from 63 data points measured in one of the side niches of the test tunnel at Yucca Mountain (Niche 4788). These permeability values are approximately lognormally distributed with a standard deviation of 0.84 in log10 (Finsterle et al., 2003). Evaluation of empirical semivariograms show weak spatial correlation, suggesting that a random correlation structure is adequate. Note that, while the variation of permeability is based on the small-scale air-injection data, the mean value of permeability is taken from the abovementioned calibration effort using surface-based boreholes. The layer-averaged calibrated value for the repository unit studied in this paper (i.e., 0.33 x 10–12 m2) is expected to be more representative over the entire repository area than the local niche measurements, since several surface-based borehole locations are available within the repository footprint. Spatially heterogeneous fracture permeability fields were generated using the specified mean and standard deviation values and mapped to the numerical mesh displayed in Fig. 2. Note that the fracture permeability was only varied within a radius of 20 m (about four drift diameters) measured from the drift centerline. Twenty meters is sufficient for the purpose of thermal seepage modeling, but small enough not to increase excessively the computational load for the simulations. Three different realizations of stochastic fracture permeability fields were generated and are analyzed in this paper.


    MODEL RESULTS
 TOP
 ABSTRACT
 INTRODUCTION
 MODELING APPROACH
 MODEL RESULTS
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Prediction of Seepage at Ambient Conditions
Before conducting the TH model analyses, predictive simulations at ambient conditions were conducted to obtain the probability of ambient seepage with the given rock properties and percolation flux boundary conditions. These simulations used the same conceptual model, numerical grid and properties as the thermal simulations, including the heterogeneous fracture permeability fields and the seepage-specific fracture capillary-strength parameter. Steady-state simulations were conducted separately for each of the three climate stages (present-day, monsoon, and glacial transition) and each of the four flux multiplication factors (1, 5, 10, and 20), in the absence of the thermal load. The magnitude of seepage derived from these simulations depends only on the capillary barrier behavior at the rock–drift interface. The ambient seepage rates were then used as reference values for comparison with the simulations performed in the presence of heat, thus providing evaluation of the effectiveness of the vaporization barrier.

Table 2 gives the ambient seepage results for the respective percolation flux scenarios in the form of seepage percentage. The seepage percentage is defined as the amount of water that seeps into the drift, divided by the amount of water that percolates with constant and uniform infiltration rate through a cross-sectional area corresponding to the footprint of the drift. The smaller this percentage, the larger the diversion capacity of the capillary barrier forming at the drift wall. Only results from one representative realization are shown in Table 2 because the other realizations show similar trends while demonstrating moderate variability. The simulation results are consistent with previous predictive analyses (e.g., Li and Tsang, 2003), where seepage was found to be extremely unlikely at average percolation conditions, and where the potential for seepage was shown to increase with the percolation rate. According to Table 2, no seepage will occur at any time for the case with a flux multiplication factor of one (with average percolation fluxes of 6, 16, and 25 mm yr–1). For multiplication factors of 5 and 10, seepage is not predicted for the present-day climate, but is expected to occur during the monsoon and glacial transition climate stages. Seepage will always occur for the extreme flux multiplication factor of 20, with percolation fluxes of 120, 320, and 500 mm yr–1 during the three climate stages, respectively. The maximum seepage percentage is 23.1%, indicating that even at high infiltration rates, the capillary barrier is capable of diverting the majority of the percolation flux.


View this table:
[in this window]
[in a new window]
 
Table 2. Predicted ambient seepage percentage for three climate stages and different flux multiplication factors.

 
Thermal-Hydrological Conditions after Waste Emplacement
Thermal-hydrological simulation runs were conducted for a period of 4000 yr after emplacement, applying the various thermal load and percolation flux boundary conditions introduced above in Boundary Conditions. Our analysis and discussion of the simulation results emphasizes those drift-scale flow phenomena that are relevant for thermal seepage. All results presented are given for one selected realization of the stochastic permeability field (i.e., the same realization as used in Table 2), if not indicated otherwise. To evaluate the role of vaporization processes in limiting seepage, we shall mainly discuss those cases where seepage conditions exist at ambient state. As shown in Table 2, these are the cases with comparably large percolation fluxes (flux multiplication factors of 5, 10, and 20).

Sample results are first provided in the form of temperature histories for all fracture gridblocks along the perimeter of the drift (Fig. 4 and 5) . The rock temperature at the drift wall is an important indicator of the effectiveness of the vaporization barrier. The figures show that the simulated temperature conditions are strongly driven by the percolation flux and the thermal load applied. For the reference thermal mode, the heat generated from the waste canisters results in maximum rock temperatures at the drift wall between about 120 and 130°C, depending on the amount of percolation flux. With a flux multiplication factor of 1, the period of above-boiling rock temperature is about 1000 yr, and the rock temperature at the end of the simulation period is still as high as 65°C. Increasing the percolation flux leads to cooler temperatures and a shorter boiling period in the fractured rock at the drift wall, thus reducing the time that vaporization can be effective in preventing seepage. This means that an increase in percolation flux would reduce the efficiency of both the capillary and the vaporization barrier.



View larger version (37K):
[in this window]
[in a new window]
 
Fig. 4. Rock temperature at the drift wall for different percolation flux multiplication factors and reference mode thermal load. For each scenario, the temperature histories in all gridblocks along the drift perimeter are depicted in the same color. Heat-pipe signatures can be observed as temperature plateaus at boiling temperature ({approx}96°C at Yucca Mountain elevation).

 


View larger version (28K):
[in this window]
[in a new window]
 
Fig. 5. Rock temperature at the drift wall for different thermal loads and a flux multiplication factor of 10. For each scenario, the temperature histories in all gridblocks along the drift perimeter are depicted in the same color.

 
The magnitude of percolation also defines the prominence of heat-pipe signatures observed in the temperature response. A heat pipe is a nearly isothermal zone maintained at about the boiling temperature, characterized by a continuous process of boiling, vapor transport, condensation, and reflux of water (Pruess et al., 1990). Such reflux processes are likely to increase the potential for seepage. Strong heat-pipe signatures occur for simulation cases with flux multiplication factors of 10 and 20, as indicated by the temperature plateaus at the boiling point of water. The duration and spatial extent of the heat pipes varies locally as a result of heterogeneity, giving rise to considerable differences between the gridblocks along the drift wall with respect to the evolution of temperature close to the boiling point and the duration of the boiling period.

Figure 5 shows the rock temperature evolution for the three different thermal load cases—the reference mode, the high-temperature mode, and the low-temperature mode. All three cases result in boiling conditions in the drift vicinity, with the maximum temperature and the duration of the boiling period depending on the respective heat load. Differences between the reference mode and the high-temperature mode (with less effective heat removal by ventilation during the first 50 yr after emplacement) occur mostly during the early heating phase, with the high-temperature mode having a higher peak temperature value and a slightly longer boiling period. Later, the differences between the two cases are marginal. In contrast, the temperatures in the low-temperature mode (smaller thermal load generated by the waste package) are considerably smaller than those of the other cases for the entire simulation period of 4000 yr. This shows that variability in ventilation efficiency (represented by the high-temperature mode) will be important only for the first several hundred years after emplacement, while variability in the heat output of the waste (represented by the low-temperature mode) will result in temperature differences that may persist for thousands of years.

As a first assessment of the potential for thermal seepage, we analyze the moisture distribution processes in response to the boiling of rock water at different times since waste emplacement. Since the fracture permeability is orders of magnitude higher than the matrix permeability, the thermally induced movement of water and gas occurs primarily in the fractures. As rock water in the matrix boils, vapor is driven into the fractures and away from the boiling zone. Condensation, causing water saturation to build up, may lead to elevated water fluxes in the fractures. The contours of fracture saturation and temperature as well as the fracture flux vectors in Fig. 6, 7, and 8 illustrate this behavior. Note that the figures show results only for the drift vicinity, while the entire simulation grid encompasses several hundred meters of rock above and below. The size of the vectors correlates with the flux magnitude at the interface between gridblocks of the numerical mesh. (The same scale for the flux vectors is used for the three figures.) Example results are presented for the reference thermal mode and a simulation run with a flux multiplication factor of 10.



View larger version (73K):
[in this window]
[in a new window]
 
Fig. 6. Fracture saturation (flooded contours), temperature (contour lines), and liquid flux vectors for flux multiplication factor 10 and reference thermal load at 100 yr after waste emplacement. The close-up view shows liquid saturation in the immediate drift vicinity.

 


View larger version (74K):
[in this window]
[in a new window]
 
Fig. 7. Fracture (flooded contours), temperature (contour lines), and liquid flux vectors for flux multiplication factor 10 and reference thermal load at 500 yr after waste emplacement. The close-up view shows liquid saturation in the immediate drift vicinity.

 


View larger version (76K):
[in this window]
[in a new window]
 
Fig. 8. Fracture (flooded contours), temperature (contour lines), and liquid flux vectors for flux multiplication factor 10 and reference thermal load at 1000 yr after waste emplacement. The close-up view shows liquid saturation in the immediate drift vicinity.

 
The fracture saturations and temperatures show an above-boiling dryout zone extending about 2 to 3 m above the drift crown at 100 yr after emplacement, which is about the time when the peak temperatures are reached (Fig. 6). Outside of the dryout zone, water saturation builds up as a result of condensation. Although this saturation buildup involves only changes of a few percent on average, it is enough to increase water fluxes considerably from the ambient unperturbed percolation. Downward flux occurs from the condensation zone above the drift back to the dryout region as a result of capillary and gravity forces. This water vaporizes, and a region of vapor–liquid counterflow may develop (heat pipe). At the same time, a significant amount of the condensate flows around the dryout zone and percolates downward, effectively removing water from the drift vicinity. Heterogeneity of permeability gives rise to strong variation in fracture saturation and water flux observed outside of the dry-out zone. Fracture saturation varies from about residual saturation up to maximum values of more than 0.11, while water fluxes range from about zero to more than 300 mm yr–1. However, as the saturation contours indicate, the effect of flow channeling is not strong enough to allow considerable penetration of water into the dryout zone. It appears that the vaporization barrier is fully effective as long as rock temperature is above boiling, despite the permeability variation. There is no water flux within the dryout zone, where fracture saturation is zero.

Figure 7 shows the local flow regime in a transition state between above-boiling and below-boiling conditions at 500 yr. The temperatures at the drift wall are just above 100°C, and the dryout zone has receded to a very small region around the drifts. No significant saturation buildup is visible outside of this region, indicating that vaporization and condensation processes are less important than before. Notice that at a few locations along the drift perimeter, the fractures close to the drift wall have started rewetting, as water has penetrated into the small dryout region and reached the drift wall. Whether seepage occurs for these "wet" conditions at the drift wall depends on the capability of the capillary barrier at the rock–drift interface to prevent water from entering the drift.

Figure 8 shows fracture saturations and flux vectors 1000 yr after waste emplacement. The overall background percolation flux at Yucca Mountain is higher at this time, since the monsoon climate period has already set in. The impact of increasing the infiltration boundary condition from 60 to 160 mm yr–1 can be clearly seen in the undisturbed flux vectors further away from the drift. Although the rock temperatures near the drift remain elevated at slightly more that 80°C, the fracture saturations and flux vectors appear similar to an ambient flow system, where the fractures close to the drift wall have rewetted and flow is diverted around the drift as a result of the capillary barrier at the drift wall. This is evident from the slightly elevated fracture saturation directly at the drift, and by the flux vectors showing water flow around the opening. Seepage can occur at these conditions if the fracture saturation is high enough and the associated capillary pressure is higher (less negative) than the threshold capillary pressure defined in Eq. [2]. Whether these conditions are met at the displayed time will be analyzed below in Evolution of Thermal Seepage.

Figure 9 allows a quantitative comparison of the fracture flow above the drifts and the magnitude of heat-induced reflux processes. To better illustrate the general trends of the flow fields at different time steps, we conducted additional thermal runs using a homogeneous fracture permeability field, so that the results are not masked by the stochastic variation. For these homogeneous runs, we plotted the downward water fluxes along a vertical line above the drift crown, close to the centerline of the drift.



View larger version (30K):
[in this window]
[in a new window]
 
Fig. 9. Fracture liquid flux vectors for flux multiplication factor 10 and reference thermal load at different times after emplacement. Solid lines show downward fluxes along a vertical line above the drift crown for a simulation with homogeneous fracture properties (drift crown is at z = 2.75 m). In comparison, the scattered symbols show results from a heterogeneous simulation at 100 yr, giving the flux magnitude in all grid elements above z = 0 m as a function of radial distance from the center of the drift (the drift wall is at r = 2.75 m).

 
Only the first two time steps, at 100 and at 500 yr after emplacement, are clearly affected by heating. At 100 yr, the maximum downward flux is about 155 mm yr–1 just outside of the dryout zone. This value is about 2.5 times larger than the ambient percolation of 60 mm yr–1, demonstrating that condensation gives rise to strong flux elevation above the drifts. However, the vertical fluxes inside the dryout zone are zero, which indicates a fully effective vaporization barrier. At 500 yr, the maximum downward flow is elevated by a smaller factor of about 1.75. It appears that the relative magnitude of flux perturbation reduces with time, corresponding to a decreasing rock temperature and a shrinking dryout zone. Thus, the maximum flux perturbation occurs (at about 100 yr in this simulation case) when heating is most intense and the vaporization barrier should be most effective. No flux elevation can be seen at 1000 yr after emplacement, when temperatures have dropped below boiling, thermal effects on fracture flow have become marginal, and the vaporization barrier capability has ceased. Here, the downward flow is equal to the imposed infiltration rate of the monsoon climate stage (160 mm yr–1), except for a small area immediately above the drift, where the water is diverted sideways as a result of the capillary barrier. Notice that the maximum heat-induced flux at early times is slightly smaller than the monsoon infiltration boundary condition. In this case of enhanced infiltration (flux multiplication factor 10), the relative impact of heating on the downward fluxes is smaller than the impact of the expected climate change.

For comparison, Fig. 9 also displays water fluxes derived from the heterogeneous simulation for the expected conditions at 100 yr after emplacement. Since the maximum fluxes do not necessarily occur along the vertical line chosen for the homogeneous runs, we show the flux magnitude in all grid elements above z = 0 m as a function of radial distance from the center of the drift. The observed trend follows roughly the solid line for the homogeneous results. However, the symbols indicate significant scatter of flux values as a result of variation in permeability, with maximum fluxes about twice as large as the homogeneous results. Also, water has penetrated slightly further toward the drift crown in some locations, as indicated by non-zero flux symbols that can be found in closer proximity to the drift. Nevertheless, the vaporization barrier remains fully effective in the heterogeneous case because all fluxes within a radius of about 1.25 m from the drift wall are still zero.

Evolution of Thermal Seepage
The thermal seepage model is capable of predicting water seepage rates from the formation into the heated emplacement drift as a function of time, simulating the combined effectiveness of the vaporization barrier and the capillary barrier during the thermally perturbed time period. In general, seepage is possible only when (i) liquid water arrives at the drift wall, depending on the vaporization impact, and (ii) the saturation at the drift wall exceeds a given threshold value, defined by the zero capillary condition in the drift and the 0.05-m vertical connection. In the following section, we first discuss the evolution of fracture saturation at the drift wall as an indicator of when and under what conditions seepage can occur. For the given fracture properties, the capillary barrier at the drift wall fails for fracture saturation values above about 0.5. At this threshold saturation, the associated capillary pressure is equal to the threshold capillary pressure defined in Eq. [2]. We then present the calculated evolution of seepage percentage for the various simulation cases and compare them with the ambient seepage results given in Table 2. Only such simulation cases are presented that are likely to have seepage conditions (i.e., high percolation fluxes).

Figure 10 shows the evolution of fracture saturation for the reference thermal load and a flux multiplication factor of 10. Results are displayed for all gridblocks located immediately at the drift wall and above the drift springline. For comparison with the heterogeneous case, the fracture saturation at the drift crown simulated from a homogeneous run is also depicted. The saturation curves demonstrate that no water arrives at the drift for about 500 yr after emplacement, corresponding to the approximate duration of the boiling period. As rock temperature decreases and the first stepwise change in infiltration occurs at 600 yr, the saturation values build up strongly, while significant variability in saturation becomes evident for the heterogeneous case. When the fractures in the drift vicinity start rewetting, local permeability contrasts promote channelized flow and create local saturation increase close to the walls. At one location, the saturation buildup is large enough to create seepage conditions after about 1400 yr after emplacement, while still in the monsoon climate stage (i.e., the saturation exceeds the threshold value for seepage to occur). With the stepwise increase of infiltration at 2000 yr, the fracture saturations increase again, and seepage starts to occur at a second location along the drift wall.



View larger version (30K):
[in this window]
[in a new window]
 
Fig. 10. Evolution of fracture saturation for flux multiplication factor 10 and reference thermal load. Dashed lines show saturation values in all gridblocks along the drift perimeter for the simulation with heterogeneous fracture properties. In comparison, the solid line gives results from a homogenous simulation, showing the fracture saturation extracted at the drift crown.

 
All saturation values increase monotonically with time, with the two jumps at 600 and at 2000 yr indicating the climate changes from present-day to monsoon and from monsoon to glacial transition climate. There is no distinct saturation peak during rewetting, which would indicate enhanced downward flow of condensate reaching the drift crown. Obviously, the enhanced downward fluxes observed in the condensation zone during the boiling period cannot reach the drift because of the effective vaporization barrier. Later, when the vaporization barrier has ceased with decreasing rock temperature, the maximum fluxes are already back to "normal," so that there is no "peak" saturation. These conditions are favored by the fact that the thermal load history displayed in Fig. 3 shows a very slow decaying behavior; that is, the transition between above-boiling and subboiling conditions is smooth. It is likely that these conditions would change for a system that is heated with full power and then abruptly shut down to zero.

Notice that the maximum saturation observed in the homogeneous case is <0.2, which is significantly smaller than the seepage threshold saturation. Clearly, fracture heterogeneity is a major factor for breaking the capillary barrier at the drift wall. Also note that the first sign of rewetting in the heterogeneous case is seen about 50 yr earlier than in the homogeneous case. This shows that heterogeneity may affect the effectiveness of the vaporization barrier because the boiling period at certain locations can be shorter than the boiling period predicted for homogeneous conditions.

The evolution of seepage into drifts for this simulation case is given in Fig. 11 , showing the total amount of water flow into the drift without distinction between different seepage locations along the perimeter. The magnitude of seepage is provided in percent, relative to the imposed percolation-flux boundary conditions. The transient curve of thermal seepage is compared with the ambient seepage percentages, which are taken from the respective simulation case in Table 2. As pointed out, the ambient seepage results were evaluated by running the thermal seepage model without thermal load until a steady state was achieved, separately for each climate period.



View larger version (24K):
[in this window]
[in a new window]
 
Fig. 11. Evolution of thermal seepage percentage for reference thermal mode and flux multiplication factor 10.

 
In agreement with Fig. 10, seepage starts to occur at about 1400 yr after waste emplacement, while in the monsoon climate period. After boiling has ceased at about 500 yr, the onset of thermal seepage is delayed by several hundred years, depending on the rather slow rewetting behavior of the fractured rock in the drift vicinity. This is not only influenced by the limited amount of percolation flux to the drift, but also by the ongoing drying of rock faces as a result of relative humidity differences between the rock and the in-drift environment. Relative humidity can still be considerably smaller than 100% inside the drifts, as the decaying radionuclides still produce significant amounts of heat.

With the stepwise increase of infiltration at 2000 yr, the seepage percentage increases considerably. At the end of the simulation period, the thermal seepage percentage is at 17%, slightly smaller than the long-term ambient value of 19.5% for the 250 mm yr–1 percolation flux. The flow system at 4000 yr is (i) still affected by the slightly elevated rock temperatures (at about 31°C) and (ii) still adjusting to the infiltration change at 2000 yr; that is, a steady-state situation has not been reached. Similar to the evolution of fracture saturation, the seepage values increase monotonically from the time of seepage onset, without showing a peak value that would indicate the possibility of heat-induced reflux processes giving rise to temporal increases in seepage. In fact, the thermal seepage values are always smaller than the ambient seepage estimates for the respective percolation flux scenario and climate period.

Note that the ambient seepage from steady-state simulation runs performed with the present-day infiltration rate (i.e., 60 mm yr–1) is zero. In other words, even without heating of the repository, the capillary barrier at the drift wall is predicted to be fully effective during the first 600 yr after waste emplacement. That no water arrives at the drift wall during these first 600 yr as a result of the vaporization barrier provides additional confidence in the barrier capabilities of the fractured rock, since two fully effective barriers are operating simultaneously and independently.

Figures 12 through 14 show the evolution of thermal seepage for additional simulation cases. Figure 12 presents results for the second realization of the random permeability field. While the absolute values of thermal and ambient seepage are smaller than for the first realization, the evolution of thermal seepage relative to the ambient values is similar to Fig. 11. Figure 13 demonstrates the effect of the different thermal loads, using the first realization of the random field. The reference mode and the high-temperature mode give almost identical results with respect to thermal seepage. This is because the temperature difference between the two thermal load cases has almost vanished at the time when thermal seepage first occurs (Fig. 5). Using the low-temperature mode has stronger impact on the transient characteristics of seepage. Here, thermal seepage is predicted to start at about 1050 yr, several hundred years earlier than in the other cases, a result of the smaller temperatures and the less effective vaporization of water. Overall, the difference in seepage percentage is relatively small among all thermal load cases.



View larger version (23K):
[in this window]
[in a new window]
 
Fig. 12. Evolution of thermal seepage percentage for reference thermal mode and flux multiplication factor 10 using a different realization of the stochastic fracture permeability field.

 


View larger version (26K):
[in this window]
[in a new window]
 
Fig. 14. Evolution of thermal seepage percentage for reference thermal mode and flux multiplication factor 20.

 


View larger version (25K):
[in this window]
[in a new window]
 
Fig. 13. Evolution of thermal seepage percentage for different thermal modes and flow focusing factor 10.

 
Figures 14 shows the importance of the percolation flux for both vaporization and capillary barrier capabilities, evidenced by the magnitude and the evolution of thermal seepage. The simulation results are obtained using a flux multiplication factor of 20, which results in infiltration rates of 120, 320, and 500 mm yr–1 imposed at the top model boundary for the three climate stages. These fluxes are rather extreme compared with the average percolation flux conditions at Yucca Mountain. As a result, the rock temperatures are generally lower and the boiling period is much shorter than in the other percolation flux cases (Fig. 4). Also, as given in Table 2, ambient seepage is predicted to occur during all three climate periods. In other words, without vaporization barrier, the seepage percentage expected for the present-day climate would be 12.2%, followed by 21.1% during monsoon climate, and 23.1% for glacial transition conditions. However, with the onset of heating, the rock temperatures rise and an effective vaporization barrier evolves shortly after the 50-yr preclosure period. As a result, the seepage percentage drops to zero, and no seepage is predicted until about 600 yr after emplacement. At this time, thermal seepage increases rapidly, mainly because the percolation flux increases from 120 to 320 mm yr–1, while the rock temperatures have significantly decreased and thermal effects have diminished. For the remaining time period, the thermal seepage percentage is relatively close to, but consistently smaller than, the ambient seepage estimates for the considered percolation flux scenario.

Figures 11 though 14 show considerable differences with respect to the magnitude and the evolution of thermal seepage. Despite these differences, there are a few important observations common to all simulation cases that were studied:

Categorization of Key Parameters for Thermal Seepage
Various additional sensitivity cases were analyzed to verify that the above main conclusions hold for the range of seepage-relevant rock properties expected at Yucca Mountain. Additional analyses included, for example, variation of the thermal conductivity, the mean value of fracture permeability, and the capillary-strength parameter. While the detailed results are not presented in this paper, these analyses also provided insight that allow categorization of the key parameters for thermal seepage as follows: (i) parameters mainly affecting the TH conditions, (ii) parameters mainly affecting the capillary barrier behavior, and (iii) parameters with impact on both the TH conditions and the capillary barrier behavior.

Thermal load and thermal conductivity, for example, belong to the first category. Varying these parameters results in considerable differences in the duration of the boiling period and the predicted maximum temperature in the rock. These conditions are important for the initiation time and the evolution of thermal seepage, but do not change the ambient seepage rate (which defines the asymptotic upper limit for thermal seepage at later stages).

Fracture capillary strength belongs to the second category. This parameter has minor impact on the TH behavior in the fractured rock, but significantly affects the asymptotic upper limit for thermal seepage at later stages. As a result of the different seepage threshold saturation, the initiation time and evolution of thermal seepage can also be affected.

The third category comprises parameters that are important for ambient seepage and also affect the intensity of TH coupling. Large percolation fluxes, for example, are typically related to large ambient seepage rates. At the same time, increased percolation flux gives rise to a reduction of temperature and a shorter duration of the boiling period. Thus, for large percolation fluxes, thermal seepage may start earlier and approach larger asymptotic values at later stages of heating. Changes in fracture permeability can also affect both the vaporization and the capillary barrier, but have counteracting effects on these barriers. Large permeabilities are generally beneficial for the performance of the capillary barrier because they allow for more flow diversion around the drifts. The vaporization barrier, on the other hand, is moderately less effective because large permeabilities can cause strong heat-pipe processes that result in lower rock temperatures and a shorter boiling period. This discussion indicates that the relationship between the relevant parameters and the seepage results can become very complicated when coupled TH processes are considered.


    SUMMARY AND CONCLUSIONS
 TOP
 ABSTRACT
 INTRODUCTION
 MODELING APPROACH
 MODEL RESULTS
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
A TH model based on the TOUGH2 simulator was developed to evaluate the potential for seepage into drifts during the time that TH conditions will be perturbed as a result of decay heat emanating from the radioactive waste to be buried in the proposed repository at Yucca Mountain. The model is designed to study the combined effectiveness of two barriers that can retard and reduce seepage: (i) the vaporization barrier caused by above-boiling conditions in the drift vicinity and (ii) the capillary barrier caused by capillary forces at the rock–drift interface. The conceptual framework for simulating the coupled TH processes is similar to the models that have accurately predicted the thermal behavior measured in large-scale underground heater tests conducted at Yucca Mountain (Birkholzer and Tsang, 2000). The specific modeling framework for the capillary barrier behavior is consistent with the method employed to simulate the seepage rates measured in ambient liquid-release experiments (Finsterle et al., 2003). The thermal seepage model is applied to various probable scenarios to cover the expected range of variability in seepage-relevant parameters. Transient seepage rates into the drift are directly calculated from the model.

Results demonstrate that the thermal perturbation of the flow field—giving rise to increased downward flux from the condensation zone toward the drifts—is strongest during the first few hundred years after emplacement, corresponding to the period when rock temperature is highest and the vaporization barrier is most effective. Even for very high percolation fluxes into the model domain, and strong flow channeling as a result of fracture heterogeneity, water is unable to penetrate far into the superheated rock during the time that rock temperature is above boiling, and model results show no seepage.

When temperatures return to below-boiling conditions and fractures start rewetting at the drift, the capillary barrier at the drift wall continues to reduce (or prevent) water seepage into the drift. Seepage is predicted to occur for simulation cases that feature strongly heterogeneous fracture permeability fields, weak fracture capillary strength in the drift vicinity, and relatively high percolation fluxes. In these cases, water starts to seep several hundred to a few thousand years after rock temperature has returned to below boiling, the delay caused by the time needed to rewet the dry fractured rock in the drift vicinity. Once the water is able to break the capillary barrier, the magnitude of seepage increases monotonically with time, and asymptotically approaches the seepage rates estimated for ambient conditions. There is no seepage "peak" as a result of enhanced reflux from the condensation zone back to the drifts.

The findings of this study, which indicate no occurrence of seepage during the period of above-boiling temperatures, can have important performance implications for the proposed high-level radioactive waste at Yucca Mountain. For example, without water seeping into drifts and contacting the waste packages, early canister failure as a result of localized corrosion is less likely. Furthermore, if there were an early failure of the waste canisters (for reasons other than corrosion), the release of radionuclides from the drifts into the rock would be limited to extremely slow diffusive transport processes, because there would be no flowing water that would pick up the released radionuclides. Thus, the vaporization and capillary barrier capabilities of the unsaturated rock in the vicinity of emplacement drifts provide an effective safety mechanism for the repository.

Three different realizations of random fracture permeability fields were analyzed in this study. Heterogeneity in fracture permeability is a key factor for seepage because local permeability contrasts promote channelized flow that (i) may penetrate into the superheated rock zone that could overcome the vaporization barrier and (ii) may create local saturation increases in the fractures close to the drift walls that could overcome the capillary barrier. Ideally, many more realizations are needed to perform a thorough statistical evaluation of the model results. This was not possible because of the vast computational effort involved in conducting TH simulations for various sensitivity cases. However, while the detailed TH processes are affected by the specifics of the stochastic fields, the general trends of thermal seepage and the differences observed between the various sensitivity cases were found to be independent of the respective realization and thus should hold true for a full statistical evaluation with hundreds of realizations. In spite of the limited number of realizations, our study elucidated critical parameters and processes that are important in predicting thermal seepage.


    ACKNOWLEDGMENTS