HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Birkholzer, J. T. Articles by Bodvarsson, G. S. Search for Related Content PubMed Articles by Birkholzer, J. T. Articles by Bodvarsson, G. S. GeoRef GeoRef Citation Agricola Articles by Birkholzer, J. T. Articles by Bodvarsson, G. S. Related Collections Heat Transport Multiphase Fluid Flow Fractured Rock

ORIGINAL RESEARCH

Evaluating the Moisture Conditions in the Fractured Rock at Yucca Mountain

The Impact of Natural Convection Processes in Heated Emplacement Drifts

^a Ernest Orlando Lawrence Berkeley National Laboratory, Earth Sciences Division, 1 Cyclotron Road, MS 90-1116, Berkeley CA 94720
^b Sandia National Laboratories, Albuquerque, NM 87185
^c University of California, Berkeley, CA
^* Corresponding author (jtbirkholzer{at}lbl.gov)

Received 14 December 2005.

	ABSTRACT

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION MODELING APPROACH MODEL APPLICATION AND RESULTS SUMMARY AND CONCLUSIONS REFERENCES

 
The energy output of the high-level radioactive waste to be emplaced in the proposed geologic repository at Yucca Mountain, NV, will strongly affect the thermal–hydrological (TH) conditions in the near-drift fractured rock. Heating of rock water to above-boiling conditions will induce large water saturation changes and flux perturbations close to the waste emplacement tunnels (drifts) that will last several thousand years. Understanding these perturbations is important for the performance of the repository, because they could increase, for example, the amount of formation water seeping into the open drifts and contacting waste packages. Recent computational fluid dynamics analysis has demonstrated that the drifts will act as important conduits for gas flows driven by natural convection. As a result, vapor generated from boiling of formation water near elevated-temperature sections of the drifts may effectively be transported to cooler end sections (where no waste is emplaced), where it would condense and subsequently drain into underlying rock units. Thus, natural convection processes have great potential for reducing the near-drift moisture content in heated drift sections, which has positive ramifications for repository performance. To study these processes, we have developed a new simulation method that couples existing tools for simulating TH conditions in the fractured formation with modules that approximate natural convection and evaporation conditions in heated emplacement drifts. The new method is applied to evaluate the impact of in-drift natural convection on the future TH conditions at Yucca Mountain in a three-dimensional model domain comprising a representative emplacement drift and the surrounding fractured rock.

Abbreviations: CFD, computational fluid dynamics • TH, thermal–hydrological

	INTRODUCTION

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION MODELING APPROACH MODEL APPLICATION AND RESULTS SUMMARY AND CONCLUSIONS REFERENCES

 
IN THE LAST FEW DECADES, extensive scientific investigations have been conducted at Yucca Mountain, NV, to explore whether the site is suitable for geologic disposal of high-level radioactive waste (Bodvarsson et al., 1999). If constructed, the repository will be located in partially saturated fractured rock several hundred meters above the groundwater table, at a depth of about 300 m below the ground surface. In the current design, cylindrical waste canisters will be distributed along horizontal tunnels of 5.5-m diameter and several hundred meters in length. The radioactive decay will produce a heat flux that will significantly change the thermal and hydrological environment at Yucca Mountain, affecting both the host rock and the conditions within the tunnels. It is expected that this thermal perturbation will prevail for thousands of years after waste emplacement (e.g., Buscheck et al., 2002; Haukwa et al., 2003). For some time, temperatures within and near emplacement drifts will be above boiling, which will give rise to complex multiphase flow and heat-transfer processes. Pore-water vaporization and subsequent condensation will lead to a large saturation and flux redistribution in the fractured rock. Vapor will also enter the emplacement drifts, thereby increasing the relative humidity near the waste packages. Understanding and predicting these changes—in the natural system as well as in the drifts—is essential for evaluating the future performance of the repository in terms of canister corrosion and radionuclide containment.

The thermally driven flow processes to be expected in the fractured rock at Yucca Mountain have been investigated in various modeling studies (e.g., Pruess et al., 1990a, 1990b; Haukwa et al., 1999, 2003; Buscheck et al., 2002; Birkholzer et al., 2004). The methodology used in these predictive studies 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, 2003). Model conceptualizations have mainly focused on large-scale average behavior for the entire repository (e.g., Haukwa et al., 1999) and on local domains featuring individual emplacement drifts with two-dimensional representations (e.g., Haukwa et al., 2003; Birkholzer et al., 2004), or have combined these two approaches in a multiscale treatment (Buscheck et al., 2002). All these models assume that gas flow processes along emplacement drifts can be neglected, because individual drifts either were not represented at all or were treated as closed systems without axial flow and transport components. As a result, the models predict that the majority of the vapor produced from boiling and evaporation of formation water remains in the fractured rock. In reality, there may be considerable vapor transport from the fractured formation into the emplacement drifts. The importance of this vapor transfer depends largely on the transport mechanism that moves vapor away from the vapor-producing sections to the cold end of the drifts.

Recent computational fluid dynamics (CFD) simulations of gas flows within emplacement drifts suggest that the in-drift conditions at Yucca Mountain will be heavily affected by heat-induced natural convection (Webb et al., 2003; Dalvit-Dunn et al., 2004; Webb and Reed, 2004; Webb and Itamura, 2004). Large-scale convection cells have been predicted to form inside the drifts in the axial direction, with gas velocities on the order of centimeters per second. Such large-scale axial circulation patterns would provide an effective mechanism for vapor transport from heated to unheated drift sections. Based on detailed CFD simulations, Webb and Itamura (2004) approximated the simulated circulation flow patterns as an axial binary diffusion process (counterdiffusion of vapor and air), using effective mass-dispersion coefficients directly estimated from the CFD flow fields. To determine the magnitude of vapor transport, the calculated dispersion coefficients were used as input to a semianalytical model for one-dimensional heat and vapor transport along drifts (Webb and Reed, 2004). The analysis demonstrated the importance of convective mixing between vapor-rich (hot) and vapor-poor (cool) drift sections. While the in-drift processes were captured in detail, however, the presence of the fractured rock was approximated using precalculated boundary conditions, thus neglecting the important feedback between the in-drift and the fractured-rock TH conditions. An alternative methodology with a more sophisticated treatment of the near-drift formation was developed by Danko and Bahrami (2003, 2004), using an iterative procedure to decouple the drift and rock-mass formulations. They applied their model to evaluate the efficiency of ventilation in emplacement tunnels.

In this study, we developed a new coupled simulation method that simultaneously handles (i) the mass and energy transport processes in the fractured rock; (ii) the mass and energy transport processes in the open drifts; and (iii) the heat, fluid, and gas exchanges at the interface between the formation and the cavity. The model concepts used for the partially saturated mass and heat transport in the fractured formation are based on existing models as described, for example, in Haukwa et al. (2003) and Birkholzer et al. (2004). Following Webb and Itamura (2004), in-drift convection was approximated as a binary diffusion process, with effective dispersion coefficients estimated from supporting CFD analyses. Mass and heat transfer between the formation and the drift were approximated by empirical boundary-layer correlations (Kuehn and Goldstein, 1976, 1978; Raithby and Hollands, 1985), based on a methodology developed by Webb and Reed (2004). The governing equations for the integrated drift and rock-mass model were solved in a fully coupled, noniterative manner, ensuring consistency between the TH conditions in the fractured rock and those in the open drift. This approach is similar to a comparative modeling study described in Buscheck (2005).

Our new, coupled simulation method was applied to study the future TH conditions at Yucca Mountain in a three-dimensional model domain comprising a representative emplacement drift and the surrounding fractured rock. Our main objective was to better understand how natural convection in drifts affects the near-drift TH processes, particularly with respect to the saturation changes and flux perturbations arising from the elevated temperatures. Sensitivity studies were conducted using two sets of effective dispersion coefficients, corresponding to different degrees of convective circulation within the open tunnels. Transient simulations were performed for the first 5000 yr after waste emplacement, comprising a short-term "hot" period, in which drift and rock temperatures are above boiling, as well as a long-term "cool" period, where drift and rock temperatures return to ambient values. Model results were compared with simulations in which the effect of natural convection was neglected. While conclusions were drawn with respect to the potential impact of our findings on repository performance, it is important to realize that our modeling approach (see below) includes a number of limiting assumptions and geometrical simplifications that are valid for a comparative evaluation, but should not be interpreted as a realistic quantitative representation of the system behavior at Yucca Mountain.

	MODELING APPROACH

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION MODELING APPROACH MODEL APPLICATION AND RESULTS SUMMARY AND CONCLUSIONS REFERENCES

 
Relevant Yucca Mountain Design Features
The current U.S. Department of Energy design for the geologic repository at Yucca Mountain envisions emplacement of waste canisters in a large number of parallel drifts spaced 81 m apart, as schematically illustrated in Fig. 1 . Typically, these drifts will have an emplacement section about 600 m long followed by unheated end sections (turnout sections), which can be as long as 100 m or more (design norm is 97 m). In contrast to the access drifts and shafts, which will be backfilled with crushed tuff material, the emplacement drifts (including their end sections) will remain open, allowing natural convection and turbulent mixing to occur. The cylindrical waste canisters, about 5 m long and between 1.75 and 2.11 m in diameter, will be distributed horizontally along the drifts, with an axial distance of 0.1 m between them. The canisters will sit on support pallets placed on the drift invert, which will be made of crushed tuff. A drip shield will cover the waste packages and protect them from rock fall and water seepage.

View larger version (71K): [in this window] [in a new window]   Fig. 1. Schematic of emplacement design at Yucca Mountain.

Thermal–Hydrological Processes Expected Near and Within Heated Emplacement Drifts
The following gives a brief overview of the heat-induced flow and transport processes expected at Yucca Mountain.

Fractured Rock Mass
The heat emanating from the waste packages will be effectively transferred to the drift walls and into the near-drift formation (Birkholzer et al., 2004). Thermal radiation will dominate the heat transport within the drift, while heat conduction will be most important in the fractured rock. As the partially saturated formation near the drifts heats up to above-boiling temperatures, the initially mostly stagnant pore water in the matrix will become mobile through vaporization (see Fig. 2a ). This will create a pressure increase, which will drive the vapor into the highly permeable fractures that are pervasive in the repository rock layers. The vapor will then move away from the boiling region through the fracture network, mostly driven by the gas pressure gradient, and a hot desaturated zone (no or only small amount of liquid water present) will develop close to the heated drift sections (e.g., Pruess et al., 1990a, 1990b; Buscheck et al., 2002; Birkholzer et al., 2004). It can be assumed that the gas pressure within the drifts will remain at atmospheric conditions, because any pressure increase would vanish in the axial direction owing to the good communication of the drift with the main access tunnels and the large rock volumes. Thus, a portion of the vapor produced in the boiling rock region will be driven back toward the drifts, thereby increasing the vapor concentration in the heated drift sections. The remaining portion will move away in a radial direction toward cooler rock regions, where the vapor eventually will condense on the fracture walls. Such condensation may give rise to elevated liquid fluxes in the fractured rock, a concern for performance assessment because of the potential for thermal seepage. At later stages, when temperatures have decreased below boiling, the rock mass near the drifts will gradually rewet (Fig. 2b). For drifts with average temperature conditions, rewetting will occur roughly a few hundred years after boiling has ceased (Birkholzer et al., 2004) a time period in which natural convection processes within the drifts are still important.

View larger version (92K): [in this window] [in a new window]   Fig. 2. Schematic of expected TH processes along emplacement drift. The figure depicts part of a drift with emplacement and end sections (modified from Webb and Reed, 2004).

Conceptual Model for Thermal–Hydrological Processes in Fractured Rock
In this study, the conceptual framework for simulating the TH conditions in the fractured rock was adopted from existing TH models for Yucca Mountain (e.g., Buscheck et al., 2002; Haukwa et al., 2003; Birkholzer et al., 2004). The relevant TH processes in the heated rock mass are the convective and diffusive movement of gaseous and liquid phases of water and air (under pressure, viscous, capillary, and gravity forces), transport of latent and sensible heat, phase transition between liquid and vapor, and vapor-pressure lowering (Pruess et al., 1999). Fluid flow was described with a multiphase extension of Darcy's law. The description of thermodynamic conditions is based on a local equilibrium model of the three phases (liquid, gas, and solid rock). The ideal gas law is assumed for the water–air mixture. Henry's law is assumed for the solubility of air in water. Vapor-pressure lowering effects are included, whereas geomechanical and geochemical processes are neglected.

The fractured rock is described using a dual-permeability concept, assuming two separate but interacting continua that superpose with each other in space. One continuum describes flow and transport in the fracture network; the other describes flow and transport in the matrix (Doughty, 1999). A continuum representation of the fractures is appropriate within the repository formation at Yucca Mountain because the fracture density is high and a well-connected fracture network forms at the scale of interest (Tsang and Birkholzer, 1999). Dual-permeability concepts assume that there are two elementary volumes at each spatial location representing the fractures and the matrix, respectively, each having a pressure, saturation, temperature, or vapor concentration value. Global flow and transport occur within both the fracture continuum and the matrix continuum, while local fracture–matrix interaction (e.g., energy transfer, imbibition, matrix diffusion) occurs between the two continua as a result of local differences in TH conditions. Thus, disequilibrium conditions between fractures and matrix and the corresponding mass and heat transfer can be efficiently modeled without explicitly accounting for all individual fractures and matrix blocks.

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 (i) flow channels may form within a fracture, and (ii) only a subset of all fractures may be actively contributing to the flow processes. To account for the possibility of flow channeling and the related reduction in 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).

Conceptual Model for In-Drift Thermal–Hydrological Processes and Coupling with Fractured Rock
In principle, the mass and heat transport processes occurring in an open drift could be modeled with a CFD simulator that would solve the mass, momentum, and energy conservation equations, including their turbulent contributions (e.g., Webb et al., 2003; Dalvit-Dunn et al., 2004; Webb and Reed, 2004; Webb and Itamura, 2004); however, solving the turbulent velocity fields expected in heated drift sections would require fine spatial and temporal resolution. Not only would this result in highly time-consuming simulation runs, but would also necessitate complex coupling approaches because of the large discretization disparities between the drift and the fractured rock. In this study, we were not interested in a detailed representation of in-drift gas flow processes. Rather, we were interested in the global large-scale vapor transport in the drift and the impact that in-drift natural convection has on the TH conditions in the near-drift fractured rock. This requires a reasonable approximation of axial vapor transport along the drifts, but does not require detailed resolution of convective flow patterns. Thus, for the purpose of this study, we followed the methodology described in Webb and Itamura (2004) and assumed that the axial transport of vapor and air can be characterized by a binary diffusion process of the air–vapor mixture, using effective mass dispersion coefficients calculated from complementary CFD flow field simulations.

To estimate effective mass dispersion coefficients for convective mixing, Webb and Itamura (2004) conducted detailed CFD simulations for a 70-m-long section of a representative emplacement drift at Yucca Mountain, using the FLUENT CFD simulator (Fluent, Inc., 2002). Two simulation cases were considered, using different temperature boundary conditions to account for the three-dimensional nature of the repository temperatures: One case assumed a uniform boundary temperature along the drift section, representing a drift in the repository center. The other case assumed a temperature gradient of a few degrees centigrade between the two ends of the drift section, representing a drift at the edge of the repository. The effective mass dispersion coefficients calculated for these cases—at 300, 1000, and 3000 yr after emplacement—are summarized in Table 1. Higher dispersion coefficients were obtained for the temperature-gradient case, reflecting the larger circulation patterns observed in the CFD simulations. In both cases, the resulting coefficients are considerably larger than the standard molecular diffusion coefficient of a vapor–air mixture of 2 x 10^–5 m² s^–1. As pointed out by Webb and Reed (2004), the two cases provide first-order estimates for analyzing the impact of convective mixing in emplacement drifts. Various simplifying model assumptions were made that suggest that the derived coefficients are probably low. Note, for example, that the CFD simulations did not include buoyancy forces induced by the density difference between vapor and air, since a single-component gas phase was assumed.

1. All drift elements representing the open enclosure (see Fig. 3) were assigned the effective mass dispersion coefficients developed by Webb and Itamura (2004). For the transient TH simulations, the coefficients determined at 300 yr were used for the postclosure period up to 600 yr after emplacement, the coefficients determined at 1000 yr were used for the period from 600 to 2000 yr, and the coefficients determined at 3000 yr were used for the remainder of the simulation period up to 5000 yr (see Table 1). The coefficients were assumed to be uniform throughout the entire length of the drift. For simplification, we furthermore assumed that the radial dispersion coefficients were identical to the axial ones. Note also that we calculated thermal conductivity in the drift based upon a Lewis number of one; i.e., we described the simultaneous natural-convection heat and mass transfer with thermal diffusivity values identical to the chosen mass diffusivity values.
   
2. Drift elements representing the open enclosure were assigned thermophysical properties (density, heat capacity) of the vapor–air mixture at prevailing temperatures and pressures.
   
3. Drift elements representing the open enclosure were assigned Darcy-type permeability values large enough to ensure that the pressure in the emplacement drifts remained close to atmospheric throughout the boiling period. In other words, the drift pressure response to the boiling-induced pressure increase in the fractured formation should be negligibly small. After some sensitivity testing, this was achieved by using a permeability value on the order of 10^–8 m², about four orders of magnitude larger than the permeability of the fractured rock mass. (The convective gas flow velocities simulated using these permeability values were several orders of magnitude smaller than the natural convection flows modeled in the CFD simulations by Webb and Itamura [2004]. This is important, because the effect of natural convection was already accounted for as a binary diffusion process.)
   
4. The waste package, the drip shield, and the small air space between them were treated in the simulation runs as a lumped entity with equivalent thermal properties (based on averaging the respective thermal properties of the waste package, drip shield, and air). Thus, the flow and transport processes occurring under the drip shield and in the gap between subsequent waste packages were neglected in this study. Also, the decay heat of the individual waste packages was averaged across the emplacement section of the drifts. For simplicity, we use the term waste package for the lumped waste package, air gap, and drip shield unit.
   
5. Thermal radiation was accounted for using a diffuse gray-body approximation, with the hotter waste package surfaces emitting and the cooler drift wall and invert surfaces receiving. The Stefan–Boltzmann law with the appropriate view factors calculated with Hottel's method (Siegel and Howell, 1972) allowed for a detailed representation of the surface-to-surface heat transfer. Wall-to-wall radiation was neglected.
   
6. The heat and mass transfer between the gas flow in the drifts and the confining drift surfaces was calculated from empirical boundary-layer correlations given in the literature. It is well known that such transfer is affected by the presence of a surficial fluid boundary layer, where the motion of particles is retarded compared with the free stream velocity (e.g., Incropera and DeWitt, 2002). Extensive experimental and theoretical work has been conducted in the past to study and describe these processes, and various empirical correlations have been developed for heat and mass transfer through the boundary layers forming in simplified geometries. Following Webb and Reed (2004), the correlations used in this study were based on work by Kuehn and Goldstein (1976, 1978) for free convection in the annular space between long horizontal eccentric cylinders, and by Raithby and Hollands (1985) for free convection between long horizontal upward-facing plates. Using these correlations, the new version of TOUGH2 calculates the heat transfer between the open drift and the confining bodies as a function of empirical heat-transfer coefficients and the local temperature gradient. Diffusive mass transfer was calculated from similar correlations, using an empirical mass-transfer coefficient and the local vapor concentration gradient. This was done for heat transfer between the waste package and the open drift, for heat and mass transfer between the open drift and rock mass (for both fracture and matrix continua), and for heat and mass transfer between the open drift and the invert. Further details on the methodology and the correlations used are given in Webb and Reed (2004).
   

View larger version (41K): [in this window] [in a new window]   Fig. 3. Schematic showing the geometry of the three-dimensional model domain (not to scale). Close-up view shows discretization of drift and drift vicinity.

Model Geometry
Three-dimensional simulation runs were performed for a simplified geometrical representation of an emplacement drift located in one of the southern panels of the repository (Fig. 3 ). In the vertical direction (z direction), the model comprised the entire unsaturated zone, having the ground surface as the upper model boundary and the groundwater table as the lower model boundary. Two symmetry assumptions were used to increase computational efficiency of the simulation runs. In the axial drift direction (y direction), symmetry allowed reducing the model to half of the drift length (plus sufficient volumes of fractured rock beyond the end of the drift to provide adequate boundary conditions). Thus, the simulated drift comprised half of the typical emplacement section length (300 m) followed by an unheated section away from the symmetry axis, here assumed to be 90 m long. The total length of the model domain in the y direction was 520 m, which included the 390-m drift plus a 130-m-long volume of fractured rock beyond the end of the drift.

Following the grid design developed by Spycher et al. (2003), symmetry assumptions can also be used to reduce the model domain in the x direction, perpendicular to the drift axis. The current repository design of parallel drifts can be represented as a series of symmetrical, identical half-drift domains with vertical no-flow boundaries between them. Thus, the numerical mesh could be reduced to a lateral width of 40.5 m, extending from the drift center to the midpoint between drifts. The model geometry was further simplified by neglecting the curvature of the drift turnout and by assuming strictly horizontal top and bottom model boundaries. These simplifications with respect to the large-scale geometry of the three-dimensional model have negligible impact on the near-drift and in-drift flow processes targeted in this study.

The numerical grid used in the simulations is shown in Fig. 3 for a vertical cross-section orthogonal to the drift axis. With the focus on in-drift and near-drift conditions, it is important to represent the drift and its vicinity with a refined discretization. In this cross-section, the drift domain comprises one finite volume for the waste package (for the lumped waste package, air gap, and drip shield unit), one finite volume for the invert, and 18 finite volumes for the annulus. Remember that in the rock-mass domain, there are two superposing continua—the fractures and the matrix blocks, respectively—with identical discretizations. In the third dimension (along the drift axis), the numerical grid comprised 28 vertical slices of varying thickness, ranging from 5 to >100 m. The entire three-dimensional grid had about 10000 finite volumes with about 35000 connections between them.

Model Boundary Conditions
Both the top and bottom boundaries of the model domain are far enough away from the repository units such that they are not significantly affected by the thermal input. (Repository units are those geologic units at Yucca Mountain in which the emplacement drifts will reside.) Accordingly, these boundaries were assigned constant Dirichlet-type conditions for temperature, with boundary values based on the geothermal gradient. Gas pressure at the top was fixed based on average barometric conditions at Yucca Mountain. Surface infiltration was modeled by imposing uniform percolation flux values that vary with time. Consistent with future climate and infiltration analyses for the Yucca Mountain (Bechtel SAIC Company, 2003b), the model considered three long-term climate periods with constant percolation: the present-day climate (up to 600 yr from emplacement of waste), the monsoon climate (600–2000 yr), and the glacial transition climate (>2000 yr). During these periods, the imposed average percolation fluxes used in this study were 6, 16, and 25 mm yr^–1, respectively, representing average infiltration conditions. The groundwater table at the bottom of the model domain was modeled as a flat, stable surface saturated with water and impermeable to gas flow. All vertical boundaries were no-flow for gas, liquid, and heat.

A Dirichlet-type boundary condition was also defined just outside the turnout section, where the emplacement drift is connected to the main access drift via a bulkhead door. Assuming good communication to the atmosphere (i.e., assuming that access drifts are backfilled with a high-permeability crushed tuff material), the main access drift was represented using a fixed gas pressure value based on the ambient gas pressure at depth. Also fixed were temperature (using the initial geothermal temperature at depth) and vapor concentration (equilibrium vapor pressure at a given temperature). These boundary conditions were constant with time; the effect of short-term barometric pumping was ignored. The bulkhead door separating the main access drift from the emplacement drift was modeled as a low-permeability, low-conductivity zone between the respective domains. (To test the impact of less effective pneumatic communication between the emplacement drift and the access drift, we have conducted additional simulation runs in which the presence of the access drift was neglected. The simulation results were hardly affected by this change in boundary condition. The fractured rock near the emplacement drifts is permeable enough to provide sufficient conduits for gas pressure release even without the presence of a nearby access drift.)

The decay heat of the radioactive waste was imposed as a uniform line load to the lumped waste package unit, thus neglecting the variability of thermal output between individual waste packages. As shown in Fig. 4 , the heat load was time dependent, exponentially decaying from an initial value of 1.45 kW m^–1 drift length. The figure also reflects the considerable heat removal by forced drift ventilation during the first 50 yr after emplacement (the preclosure period). Similar to previous studies (e.g., Buscheck et al., 2002; Birkholzer et al., 2004), forced drift ventilation was not explicitly modeled in our simulation. Rather, the heat losses by forced drift ventilation were incorporated by reducing the effective heat output in the simulation to about 14% of its original value. Design estimates for the ventilation system suggest that, on average, 86.3% of the decay heat can be removed during the preclosure period (Walsh et al., 2004).

View larger version (11K): [in this window] [in a new window]   Fig. 4. Thermal load as a function of time (from Bechtel SAIC Company, 2003c). Time zero represents the time of waste emplacement.

View this table: [in this window] [in a new window]   Table 2. Hydrogeological and thermal input values of fractured rock mass.

	MODEL APPLICATION AND RESULTS

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION MODELING APPROACH MODEL APPLICATION AND RESULTS SUMMARY AND CONCLUSIONS REFERENCES

 
Our model simulations studying coupled TH processes in the drifts and the fractured rock mass covered the first 5000 yr after emplacement of radioactive waste. Being interested in the effects of natural convection, we focused on the postclosure period, which follows the forced-ventilation preclosure period of 50 yr, and ran the simulations until the impact of natural convection became less important because of decreasing temperatures. Two main simulation cases were considered using two sets of effective dispersion coefficients: Case 1 represents strong convective mixing in emplacement drifts; Case 2 represents moderate convective mixing (see effective dispersion coefficients in Table 1). For comparison, we also studied a simulation case where axial transport along the drift was neglected (Case 3). This was achieved by removing all open drift elements from the numerical grid, thus making the drift cavity an impermeable boundary to the fractured rock mass. The only remaining in-drift elements in Case 3 were those representing the lumped waste package units, which were connected to the drift wall to allow for radiative heat transfer.

The initial conditions for the postclosure simulation runs were derived from simplified model simulations for the preclosure period, where forced ventilation removed a significant amount of the decay heat. As pointed out above, the heat losses from forced ventilation were accounted for by reducing the effective heat generation of the radioactive waste according to Fig. 4. The possible impact of forced ventilation on the moisture content in the rock mass was neglected in these simplified initialization runs.

Simulation results for the postclosure period are first provided in the form of temperature profiles along the drift, extracted in a volume element just above the center of the drift. Figure 5 shows different times during the postclosure period—100, 500, 1000, 2000, and 5000 yr after emplacement—for Cases 1 and 2. The differences between the two cases are rather subtle, indicating that the transport of sensible and latent heat with the gas flow in the drift is not very important for the in-drift temperature distribution. Temperatures for the case with strong convective mixing are slightly lower in the heated section and slightly higher in the end section of the drift, a result of the more effective natural convection transport of heat in the axial direction.

View larger version (21K): [in this window] [in a new window]   Fig. 5. Temperature profiles along emplacement drift, extracted in a volume element just above drift center. The distance along the drift is measured with y = 0 in the center of the heated section of the drift, the symmetry boundary of the half-drift model. The heated section of the drift ends at y = 300 m (indicated by vertical line), followed by a 90-m-long unheated end section.

To evaluate the heat-related moisture redistribution processes, we have plotted in Fig. 6 and 7 the TH conditions for a vertical cross-section along the drift. We selected the simulation results at 500 yr after emplacement, representative of the expected situation toward the end of the superheated period at Yucca Mountain. For better visualization, we focused on the rock mass in the drift vicinity (with vertical extension from –40 z 40 m), and used a vertical length scale five times the horizontal.

View larger version (38K): [in this window] [in a new window]   Fig. 6. Simulated thermal–hydrological conditions for Case 1 (strong convective mixing) in vertical cross-section along the drift after 500 yr. (a) Colored contours show fracture saturation; contour lines show temperature in the rock mass and in the drift. Arrows depict relative magnitude and direction of liquid fluxes. (b) Colored contours show matrix saturation; contour lines show vapor concentration (vapor mass fraction) in the matrix and in the drift.

View larger version (39K): [in this window] [in a new window]   Fig. 7. Simulated thermal–hydrological conditions for Case 2 (moderate convective mixing) in vertical cross-section along the drift after 500 yr. (a) Colored contours show fracture saturation; contour lines show temperature in the rock mass and in the drift. Arrows depict relative magnitude and direction of liquid fluxes. (b) Colored contours show matrix saturation; contour lines show vapor concentration (vapor mass fraction) in the matrix and in the drift.

Before evaluating heat-related flow and transport processes, we should briefly discuss the ambient (initial) situation in the partially saturated repository units at Yucca Mountain, where natural infiltration is very small (on the order of a few millimeters per year). As pointed out above, a present-day infiltration of 6 mm yr^–1 was assumed in the model. As a result of the small infiltration into the unsaturated zone, the fractures are essentially drained of water and nonconductive, with the saturation value close to residual. The matrix pores, on the other hand, can hold a significant amount of water owing to strong capillary forces, with saturation values of the order of 0.85. Despite the relatively large saturation, however, the water in the matrix is almost stagnant at ambient conditions because of the very low matrix permeability. Before heating, the temperature is about 25°C at the repository depth, and the gas pressure is close to atmospheric pressure at the elevation of Yucca Mountain (about 89000 Pa). The ambient vapor concentration in the gas phase is very small, corresponding to a relative humidity of 100% at 25°C.

According to our simulation results, these ambient conditions will be strongly altered at 500 yr after emplacement by boiling and subsequent condensation of pore water. Both the matrix and the fracture saturation contours in Fig. 6 show that a large volume of rock has desaturated in the vicinity of the heated section of the drift, with the region of reduced liquid saturation extending several meters into the formation. The desaturated region below the drift is larger than above, where further saturation decrease is somewhat limited by the percolation flux arriving at the boiling zone and the gravity-driven reflux of condensate. Note that desaturation is less common in the matrix continuum than the fractures, for two reasons. First, there is significantly more ambient water in matrix pores that would need to boil off. Second, vaporization leads to a pressure increase in the matrix, since the small matrix permeability limits the vapor release from the pores. As pressure builds up, the boiling point of water increases, allowing the presence of liquid water in the matrix at rock temperatures >96°C. (For comparison, the gas-pressure increase from 89000 Pa, at atmospheric, to 110000 Pa, the predicted maximum increase in the rock matrix, would result in a boiling temperature of about 102°C.)

Where does the vapor produced from boiling migrate and where does condensation occur? We answered these questions by analyzing regions of elevated saturation (i.e., elevated above the ambient saturation values) as well as by evaluating vapor concentration contours, expressed as vapor mass fraction. Figure 6b indicates that the vapor concentration in the gas phase has increased significantly in the heated rock regions, with maximum values above 0.7 a few meters away from the drift. Part of the vapor moves effectively within the highly permeable fractures out into the distant rock mass. As a result, elevated vapor concentrations of 0.3 and more can be seen at distances of 30 m above and below the heated drift. That condensation of vapor occurs in these cooler rock regions can be inferred from the elevated matrix saturations, which reach maximum values well above 0.9. This saturation increase is caused by condensation of vapor on the fracture walls and subsequent imbibition into the matrix pores. The fracture continuum does not exhibit a corresponding saturation buildup because the bulk of the condensate moves effectively away in the permeable fracture network as a result of the small fracture capillarity.

The saturation buildup in the fractured rock mass near the unheated end section of the drift is more pronounced, indicating that a significant amount of the vapor produced in the rock mass enters the drift, migrates in the axial direction (as a result of circulating gas flows from natural convection), and eventually condenses on the cooler drift wall surfaces. Most of the condensate flows from the drift walls toward the drift bottom through the invert and drains downward in the fracture continuum. As a result, fracture saturation increases strongly below the drift, and significant gravity-driven liquid fluxes occur, driving water away into the underlying formations (Fig. 6a). A small fraction of the condensate imbibes from the drift walls into the matrix pores as a result of capillary forces. The impact of such imbibition can be seen in the elevated matrix saturations both above and below the drift end section (Fig. 6b). Despite the saturation increase, the liquid fluxes in the matrix are much smaller than those in the fracture continuum. In fact, when plotted on the same length scale, the flux vectors in the matrix are too small to be visible in Fig. 6b.

While vapor-rich gas is driven from the heated drift sections toward the unheated drift turnout, vapor-poor gas moves back, thereby reducing the in-drift vapor concentration and relative humidity near the heat-producing waste packages. This effect, modeled in this study as an effective binary diffusion process, is present along the entire length of the emplacement section, causing a strong vapor concentration difference between the rock mass and the drift. For example, Fig. 6b shows that the maximum vapor concentration in the drift is between 0.2 and 0.3, compared with a maximum value of >0.7 in the rock. As a result of this concentration difference, vapor migrates from the heated rock mass into the drift not just by pressure-driven convective flow, but also by diffusion. More detail on the magnitude of the vapor exchange between the fractured rock mass and the drift, as well as on the relative importance of convective vs. diffusive transport, are presented below.

We may now compare the simulation results for Case 1 in Fig. 6 with the simulation results for Case 2 in Fig. 7. Whereas the heat-induced moisture redistribution processes are similar in general, Case 2 differs from Case 1 in that the convective mixing between vapor-rich and vapor-poor regions of the drift is less effective. As a result, more vapor remains in the heated rock mass where it is produced, and less vapor migrates from the heated drift section toward the cooler drift turnout. The former is evident in Fig. 7 in the smaller rock-desaturation region above and below the drift, the higher matrix saturation values in the rock-mass condensation zone, and the higher vapor concentrations in both the fractured rock and the heated drift section. The latter shows in the less pronounced fracture and saturation buildup near the drift end section, indicating that less vapor migrates there and condenses on the drift walls. Also, the liquid flux vectors below the drift are considerably smaller in Case 2 than Case 1 (the vector scales in both figures are identical).

Note that the vapor concentration contours above and below the drift show a distinct asymmetry in both Fig. 6b and 7b, with higher values occurring above the drift. One reason for this asymmetry is the already mentioned gravity effect, with less water boiling below the drift than above. More important are the lower temperatures below the drift, a result of the small thermal conductivity of the crushed-tuff invert. Because conduction through the invert is less effective than direct thermal radiation, the temperature increase of the fractured rock below the invert is slower than at other locations along the drift perimeter. Lower temperatures cause a decrease in the saturated vapor pressure, which in turn reduces the maximum possible vapor concentration in the gas phase.

Figures 6 and 7 demonstrate clearly that in-drift natural convection strongly affects the TH conditions in the fractured rock mass at 500 yr after emplacement. Below, we shall evaluate the flow and mass-transport processes for the entire simulation period of 5000 yr. We start with the in-drift conditions for Cases 1 and 2 (i.e., for strong vs. moderate convective mixing). Vapor concentration profiles along the drift are depicted in Fig. 8 at different times during the postclosure period, extracted in a volume element just above the center of the drift. Three conclusions may be drawn: (i) Cases 1 and 2 show distinct differences not just at 500 yr, as shown above, but at all times except at 5000 yr, when the vapor concentration along the drift is small in both cases; (ii) the highest vapor content in the drift occurred at early heating stages, when boiling of pore water was most significant; and (iii) the vapor concentrations in the heated drift section remained elevated long after the rock mass had returned to below boiling.

View larger version (23K): [in this window] [in a new window]   Fig. 8. Vapor concentration along drift just above drift center for (a) Case 1, with strong convective mixing, and (b) Case 2, with moderate convective mixing

View larger version (22K): [in this window] [in a new window]   Fig. 9. Relative humidity along drift just above drift center for (a) Case 1, with strong convective mixing, and (b) Case 2, with moderate convective mixing

The question arises, how would the in-drift moisture conditions compare with Case 3, in which axial transport along the drift was neglected? Since Case 3 was modeled by removing all open drift elements from the simulation, information on the in-drift vapor distribution is not directly available from our simulation results. It is reasonable to assume, however, that the in-drift vapor concentrations would be in close equilibrium with those in the adjacent fractured rock (because there is no axial transport driving vapor away). We have therefore plotted in Fig. 10 the vapor concentration and relative humidity values extracted in the fractured rock just above the drift crown at 500 yr after emplacement for the three simulation cases. Both Cases 1 and 2 exhibit strongly reduced vapor concentration and relative humidity conditions compared with Case 3, indicating the importance of enhanced vapor transport even when moderate convective mixing is assumed. In Case 3, the relative humidity in the fractured rock (and thus in the drifts) is already at 100% throughout a 70-m section near the unheated end of the drift. This section corresponds with the section where the rock temperature has decreased below the boiling point (compare with Fig. 5). In other words, without axial vapor transport, the rock and in-drift gas phase become vapor saturated soon after boiling in the rock has ceased and the fractured rock has started to resaturate.

View larger version (28K): [in this window] [in a new window]   Fig. 10. (a) Vapor mass fraction and (b) relative humidity along drift at 500 yr after emplacement. Dashed lines give vapor mass fraction and relative humidity extracted in fracture element just above drift crown. For comparison with Fig. 8 and 9, solid lines give in-drift conditions extracted just above drift center.

View larger version (21K): [in this window] [in a new window]   Fig. 11. Evolution of vapor flux from formation to drift integrated across heated drift section for (a) Case 1, with strong convective mixing, and (b) Case 2, with moderate convective mixing. Vapor flux was normalized by dividing by the ambient percolation flux integrated across the cross-sectional area of the heated drift section.

In addition to the total vapor flux into the drift along the heated section, we also plotted the total liquid flux leaving the unheated end section of the drift into the fracture continuum. The total liquid fluxes leaving the drift into the fractures are slightly smaller than the total vapor fluxes entering the drift. This indicates that a small fraction of the vapor condensing in the unheated end section is imbibed into the matrix pores, where the liquid remains mostly immobile because of the small matrix permeability (compare with Fig. 6b and 7b). The majority of the vapor, however, drains off into the underlying formations.

Figure 11 furthermore distinguishes between convective and diffusive transport of vapor between the formation and the drift. Convective transport results from the pressure increase in the fractured rock as pore water boils off, while diffusive transport is caused by the vapor concentration gradient. In Case 1, diffusion is the dominant transport mechanism, since the strong mixing along the drift ensures small vapor concentrations and thus a large concentration gradient. In Case 2, vapor diffusion is less important, while the magnitude and evolution of convective transfer is comparable to Case 1. Thus, the differences between Case 1 and Case 2 are mostly a result of the more or less effective diffusive exchange. Note that in both cases, convective transport is only effective for the first 2000 yr after emplacement.

Let us now analyze in more detail how the observed in-drift conditions affect the TH situation in the near-drift fractured rock. It is useful to plot the evolution of matrix liquid saturation along the drift for selected times (Fig. 12 ), considering all three simulation cases. The figure shows clearly how important it is to account for in-drift transport in TH simulations of the rock mass at Yucca Mountain. In Case 3, with no vapor transfer between the formation and the drift, the rock mass near the drift resaturates much earlier after the initial impact of vaporization and arrives at much higher long-term saturation levels along the emplacement section. These observations are particularly apparent in comparison with Case 1 (strong convective mixing). At 500 yr, for example, the matrix pores along the emplacement section are still desaturated in Case 1, while Case 3 already exhibits considerable saturation buildup. Further, in Case 1 the maximum saturation at 5000 yr is <0.7, indicating that evaporation remained effective, while Cases 2 and 3 ends up with maximum saturations of 0.94 and 0.97, respectively. Notice the local minimum in the saturation profiles of Case 2 at 1000, 2000, and 5000 yr, which corresponds to a similar minimum in the relative humidity profiles in Fig. 9b. Convective mixing reduces the moisture content in the rock matrix, but only across a limited drift length close to the end of the heated section.

View larger version (21K): [in this window] [in a new window]   Fig. 12. Matrix saturation along drift just above drift crown for (a) Case 1, with strong convective mixing, (b) Case 2, with moderate convective mixing, and (c) Case 3 (no axial vapor transport).

View larger version (17K): [in this window] [in a new window]   Fig. 13. Evolution of desaturation volume of rock mass along heated drift section. Threshold used for desaturation was a matrix saturation <0.425, half of the average ambient saturation.

It is of interest to determine how these changes in TH conditions affect the heat-induced perturbation of liquid fluxes in the fractured rock. Figure 14 shows the downward liquid fluxes at different times along a horizontal profile chosen close to the center of the heated drift section. The profile is perpendicular to the drift axis (x direction); it starts at the sidewall of the drift and extends laterally several meters into the fractured rock. All fluxes are given relative to the ambient natural percolation flux at the considered time; i.e.; fluxes >1 denote a flux enhancement. At early heating stages (i.e., 100 and 500 yr), vaporization of pore water and subsequent condensation lead to a considerable flux enhancement just outside of the desaturation zone where the condensate drained off. Note the strong differences between the different simulation cases, with Case 1 (strong convective mixing) showing much less flux enhancement than the other cases because more vapor moves out of the fractured rock into the drift. Flux enhancement during boiling is a concern for repository performance because it may increase the potential for preferential flow events penetrating the desaturation region and possibly dripping into the heated drift (Birkholzer and Ho, 2003; Birkholzer et al., 2003). Later, when boiling has ceased and the fractured rock near the drift resaturates, the maximum vertical fluxes occurs close to the drift wall. These fluxes are enhanced because water—arriving at the drift crown from natural percolation—is diverted sideways and around the drift as a result of the capillary barrier at the drift wall. Again, the flux enhancement is smallest in Case 1, since liquid water present near the drifts is subject to evaporation. Compared with the other cases, strong desaturation can be seen at 1000 and at 2000 yr in Case 1; less pronounced is the impact of evaporation at 5000 yr. These results correlate well with the evolution of the vapor transfer between the fractured rock and the drift.

View larger version (21K): [in this window] [in a new window]   Fig. 14. Vertical liquid fluxes in horizontal profile from drift wall into fractured rock in a drift section close to the center of the drift. Fracture and matrix fluxes have been summed up and divided by ambient percolation flux at the considered time.

	SUMMARY AND CONCLUSIONS

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION MODELING APPROACH MODEL APPLICATION AND RESULTS SUMMARY AND CONCLUSIONS REFERENCES