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 ISI Web of Science (1) Citing Articles via Google Scholar Google Scholar Articles by Birkholzer, J. T. Articles by Zhang, Y. Search for Related Content PubMed Articles by Birkholzer, J. T. Articles by Zhang, Y. GeoRef GeoRef Citation Agricola Articles by Birkholzer, J. T. Articles by Zhang, Y. Related Collections Heat Transport Dual Porosity/Permeability Models Fractured Rock

ORIGINAL RESEARCH

The Impact of Fracture–Matrix Interaction on Thermal–Hydrological Conditions in Heated Fractured Rock

Ernest Orlando Lawrence Berkeley National Lab., 1 Cyclotron Rd., MS 90-1116, Berkeley, CA 94720
^* Corresponding author (jtbirkholzer{at}lbl.gov)

Received 1 June 2005.

	ABSTRACT

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Dual-continuum models have been widely used in modeling flow and transport in fractured porous rocks. Among many other applications, dual-continuum approaches were used in predictive models of the thermal–hydrological conditions near emplacement tunnels (drifts) at Yucca Mountain, NV, the proposed site for a radioactive waste repository in the USA. In unsaturated formations such as those at Yucca Mountain, the magnitude of mass and heat exchange between the two continua—fracture network and rock matrix—depends on the small-scale flow characteristics in the fractures, because channelized finger-type flow strongly reduces the interface area between the matrix surfaces and the flowing liquid. This effect may have important implications, for example, during the time period that the fractured rock near the repository drifts would be heated above the boiling point of water. Depending on the magnitude of heat transfer from the matrix, water percolating down the fractures will either boil off in the rock region above drifts or may penetrate all the way to the drift walls and possibly seep into the open cavities. In this study, we conducted a sensitivity analysis using different approaches to treat fracture–matrix interaction in a three-dimensional dual-continuum setting. Our simulation example was a laboratory heater experiment described in the literature that provides evidence of rapid water flow in fractures, leading to seepage into the heater hole despite above-boiling conditions in the adjacent fractured rock. We showed that the experimental finding can only be reproduced when the interface area for heat transfer between the matrix and fracture continua is reduced to account for flow channeling. Our analysis also suggests that the conditions in the laboratory experiment are unique and not representative of the expected thermal–hydrological conditions near emplacement drifts at Yucca Mountain.

Abbreviations: CNWRA, Center for Nuclear Waste Regulatory Analyses • TH, thermal–hydrological

	INTRODUCTION

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
ANALYSIS of the proposed underground repository for high-level radioactive waste at Yucca Mountain, NV, requires knowledge of water flow in the unsaturated fractured rock at elevated temperatures because of the heat emanating from the decaying waste. The temperatures near the emplacement tunnels (drifts) are expected to be above boiling for several hundred years after the waste is stored. If water penetrates through the above-boiling rock regions above the repository and seeps into the emplacement tunnels, it can promote the corrosion of waste packages and enhance the release of radionuclides. However, because the liquid water, which is mostly flowing in the well-connected rock fractures, will be subject to effective vaporization from the adjacent rock matrix, the above-boiling rock zone forming above the repository may significantly reduce the possibility of water contacting waste packages (Ramspott, 1991; Wilder, 1993).

The thermal–hydrological (TH) conditions in the unsaturated environment at Yucca Mountain following waste emplacement have been evaluated in various modeling studies (e.g., Pruess and Tsang, 1993, 1994; Haukwa et al., 1999; Buscheck et al., 2002). Most of these studies use some form of dual-continuum approach that separately describes the TH conditions in the fractures and the rock matrix because of their strongly different properties. It was generally suggested that the above-boiling rock regions expected near waste emplacement drifts will form an effective vaporization barrier, preventing percolating water from arriving at the drift walls. Thus, as demonstrated in Bechtel SAIC Company (2004a) and Birkholzer et al. (2004), thermal seepage of water into drifts will not be of concern during the above-boiling phase of the repository. (In the context of this study, thermal seepage is defined as seepage of water into drifts under thermally perturbed hydrological conditions.)

Recent experimental results, however, may pose a challenge to the conclusions about thermal seepage at Yucca Mountain, suggesting that the exclusion of liquid water from above-boiling rocks may not be complete. A laboratory heater test conducted by the Center for Nuclear Waste Regulatory Analyses at the Southwest Research Institute in San Antonio (hereafter referred to as the CNWRA experiment) indicated that water rapidly flowing in vertical fractures was able to penetrate the above-boiling rock region and seep into a horizontal cylindrical opening (Green and Prikryl, 1998, 1999; Green et al., 2003). In other words, thermal seepage was shown to be possible despite above-boiling temperatures. In contrast to the in situ heater tests at Yucca Mountain, where natural percolation is relatively small, the CNWRA experiment used artificial water release from the top of the experimental apparatus, testing the thermal seepage potential for rather extreme flow conditions.

The question arises whether the simulation models used for predicting TH conditions at Yucca Mountain would be capable of explaining the observation of thermal seepage in the CNWRA experiment. One concern about these models is the use of continuum approaches, which generally involves some sort of averaging and homogenization of heterogeneous formation properties, and tends to underestimate the probability of small-scale processes, such as the formation of channelized, finger-type flow paths. As a result, the amount of heat available to vaporize the water flowing in fractures may be overestimated because the full geometric contact area is used for energy transfer from the adjacent rock matrix rather than the area where the flow actually occurs. In general, rapid channelized flow events may be more likely to overcome the vaporization barrier. This possibility was demonstrated, for example, using analytical (Phillips, 1996) and semianalytical solutions (Birkholzer, 2003; Birkholzer and Ho, 2003; Birkholzer et al., 2003) that investigated the fate of discrete finger-type flow events penetrating a vertical fracture surrounded by an above-boiling rock matrix. Since these solutions use various geometric and conceptual simplifications, they are mostly limited to conceptual studies and cannot be used for the complex modeling applications necessary for predicting Yucca Mountain's future TH conditions.

Regarding isothermal flow and transport, various researchers have improved classic dual-continuum approaches to account for channelized flow in fractures and the related reduction in fracture–matrix mass transfer. The improved approaches adjust the interface area between fractures and matrix, depending on the flow characteristics in the fractures. The adjusted interface area represents the fraction of the full geometric interface that actively participates in fluid and solute exchange (e.g., Ho, 1997; Doughty, 1999; Liu et al., 1998, 2003). Different approximate formulations were proposed for fracture–matrix interface reduction, including, for example, multiplication with a constant factor, with a power function of fracture liquid saturation, or with fracture relative permeability, summarized in Doughty (1999), as well as multiplication with a factor derived from the active fracture model developed by Liu et al. (1998).

Significantly less research work regarding interface reduction approaches has been conducted for heat transfer problems, possibly because the concept of channelized fracture flow in heated rocks involves a local thermodynamic disequilibrium between the liquid and the gas phase in the fractures, a conceptual challenge for most multiphase simulators. In a prior study (Birkholzer et al., 2004), we introduced the so-called thermal seepage model, a dual-continuum model specifically designed to evaluate the potential of drift seepage during the period of thermal perturbation at Yucca Mountain. As a conservative choice regarding the potential for thermal seepage, the model uses a small, constant factor reducing the amount of heat transferred between heated matrix blocks and fracture water (Bechtel SAIC Company, 2004a). The predictive simulations conducted with the thermal seepage model suggested that, for the conditions at Yucca Mountain, channelized flow is not likely to result in seepage during the time that the matrix temperature is above boiling, because water is unable to penetrate far into the above-boiling rock region despite the strongly reduced fracture–matrix interface area. Seepage may only occur at subboiling temperatures, depending on the amount of percolation flux arriving at the drifts and the near-field properties of the fractured rock (Birkholzer et al., 2004).

The objective of this study was to improve our understanding of thermal seepage at above-boiling conditions by evaluating the applicability and impact of fracture–matrix interface reductions in heat-transfer problems. We adopted the general modeling concept of the thermal seepage model used for Yucca Mountain predictions (Birkholzer et al., 2004), applied different fracture–matrix area formulations known from ambient flow and transport studies, and tested the sensitivity of the TH conditions to these formulation choices. The sensitivity study was conducted using the CNWRA experiment as a simulation example. Our goals were (i) to demonstrate that the thermal seepage model is capable of explaining the observations from the laboratory experiment where thermal seepage occurred despite above-boiling rock temperatures, (ii) to identify the fracture–matrix area formulation best suited for this purpose, and (iii) to evaluate whether the test conditions and observations in the CNWRA experiment are representative of the TH conditions expected at Yucca Mountain.

	BASIC PROCESSES AND DUAL-CONTINUUM MODELING CONCEPTS

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
Thermal-Hydrological Processes in Fractured Porous Rock at Yucca Mountain
The fractured tuff in the hydrogeological units hosting the proposed repository at Yucca Mountain is characterized by strong intrinsic heterogeneity. The tuff matrix has a rather large porosity between 0.1 and 0.2, but a very small permeability (on the order of 10^–18 m²). On the other hand, the formation is intensely fractured, with fracture spacings of a few decimeters or less, and the continuum permeability of the fractures is many orders of magnitude larger than that of the rock matrix. Since the fractures occupy only a very small volume fraction of the formation, the porous matrix stores the vast majority of the mass and energy in the system and accounts for most of the heat conduction. Initially, before the start of heating, the matrix pores are about 80 to 90% saturated with water, whereas the fractures are essentially dry (Birkholzer and Tsang, 2000).

The heat generated by the decaying radioactive waste at Yucca Mountain will induce strong TH processes in the fractured rock, which are schematically illustrated in Fig. 1 (Birkholzer et al., 2004). After emplacement of the waste, heat is transferred from the waste packages to the drift walls and subsequently migrates radially outward, mostly by conduction in the matrix. As temperatures in the fractured rock exceed the boiling point of water, the stagnant water in the matrix pores becomes mobile by vaporization. Most of the vapor moves away from the vicinity of the drift through the permeable fracture network. With continuous heating, an above-boiling dry-out region will develop near the drifts, where the matrix pores and the fractures are essentially desaturated, while condensation of vapor occurs in the fractures farther away from the heat source. Condensation above the heated areas may result in gravity-driven flow of water back toward the heated drift. Whether refluxing condensate boils off in the dry-out zone or is able to reach the drift wall depends on the extent of the above-boiling region, the intensity of heating, and the small-scale flow characteristics within the fractures. As mentioned above, water is more likely to penetrate into the above-boiling dry-out zone if the condensate reflux occurs in the form of channelized, finger-type flow. Water that reaches the drift wall may or may not seep into the open drift, depending on the effectiveness of capillary forces at the rock–drift interface that tend to keep water in the formation (e.g., Finsterle et al., 2003; Ghezzehei et al., 2004).

View larger version (72K): [in this window] [in a new window]   Fig. 1. Schematic of expected thermal–hydrological conditions near emplacement drift at Yucca Mountain (not to scale).

View larger version (87K): [in this window] [in a new window]   Fig. 2. Schematic illustration of channelized flow in a fracture–matrix system.

Channelized flow in fractures may not only affect the interface area to be used in dual-continuum models, but also the details of mass or heat transport in the adjacent matrix. For flow occurring within the entire fracture plane, the mass and heat flow in the adjacent matrix can often be approximated as a one-dimensional process perpendicular to the fracture surface. This can be different for finger-type flow in the fracture, because the local perturbation in the matrix may be radial rather than one-dimensional.

Dual-Continuum Models and Fracture-Matrix Interaction Approaches
Dual-continuum models were developed in the early 1960s to study flow processes in fractured porous rock (e.g., Barenblatt et al., 1960). The basic idea is that the heterogeneous components of a fractured porous formation are separated into two overlapping continua, one representing the fracture system, one representing the porous rock matrix, locally coupled by exchange terms to describe interaction between the two systems. During the past decades, dual-continuum models have been widely applied to various simulation problems in fractured porous rock, mostly involving flow and contaminant transport. Because the matrix permeability is orders of magnitude smaller than in the fractures, the global (matrix-to-matrix) flow and transport processes within the matrix continuum have often been neglected. In this case, the matrix continuum provides a local storage function for the fracture continuum, which in turn accounts for all global flow and transport processes. These simplified dual-continuum models are referred to as dual-porosity approaches.

The local coupling between the fracture and matrix continua has often been described by the quasi-steady transfer term initially introduced by Warren and Root (1963), assuming a linear relationship between the state variables of the fracture and matrix continua, respectively. This means that the state variables in the matrix are represented by average values, independent of their actual spatial distribution within a matrix block. Since quasi-steady approaches do not capture details of the local-scale flow and transport processes near the fracture–matrix interface, the transfer term often includes a best-fit coefficient, loosely correlated to the geometry of the matrix blocks (e.g., Gerke and van Genuchten, 1993). More sophisticated coupling concepts have been developed using higher order exchange terms (e.g., Zimmerman et al., 1990; Dykhuizen, 1990; Zimmerman et al., 1993) or solving a one-dimensional equation for local flow and transport in representative matrix blocks (e.q., Bibby, 1981; Huyakorn et al., 1983). While these improved coupling concepts allow a better representation of transient responses in the matrix (particularly the early-time behavior with steep gradients), they generally assume that there is no global matrix-to-matrix flow and transport to allow for simplified representations of the local state variable; i.e., they can only be applied in combination with dual-porosity formulations.

At Yucca Mountain, the global (matrix-to-matrix) flow and heat transport processes in the matrix cannot be ignored. For example, while the matrix-to-matrix contributions are moderately relevant for variably saturated flow (mostly at low overall water content), they are extremely important for heat transfer, since heat conduction occurs predominantly in the matrix. Thus, applications of the dual-continuum method to flow and heat transport at Yucca Mountain have usually involved a general dual-continuum setting with quasi-steady transfer terms for local interaction between fractures and matrix. While the quasi-steady concept has limitations in capturing local-scale details in the rock matrix, these limitations are arguably less important for dual-continuum simulations than the main focus point of this study, which is the impact of fracture–matrix interface reduction as a result of strongly focused, finger-type flow in the fractures.

For two-phase flow in gaseous and liquid phases, the gas and liquid transfer between the fracture and matrix continua can be readily described by quasi-steady transfer terms that have been extended from the original Warren–Root model to account for variably saturated flows. For example, the flow rate of liquid water from the fracture continuum to the matrix continuum can be written as follows (e.g., Ho, 1997; Gerke and van Genuchten, 1993; Liu et al., 1998):

[1]

where S_L is liquid saturation, K_L(S_L) is the unsaturated hydraulic conductivity for liquid water between a fracture element and the corresponding matrix element, A_FM is the geometric area between all fractures and matrix blocks lumped into a computational element, P_cF(S_L) and P_cM(S_L) are the capillary pressures of fracture and matrix elements, respectively, d is the average dimension of the matrix blocks, and ß is the best-fit coefficient. Finally, R_L denotes the interface reduction factor for liquid flow, which accounts for the impact of fingering or nonactive fractures on fracture–matrix interaction. Note that expressions similar to Eq. [1] can be written for the gas flow between the two continua, using an interface reduction factor for gas flow, R_G. However, the interface reduction for gas flow is often neglected in unsaturated media problems because its effect is small (Doughty, 1999).

Different formulations for R_L have been proposed, including a constant factor, a power function of liquid saturation, or the liquid-phase relative permeability (Ho, 1997; Doughty, 1999). In the first formulation, the constant multiplier is typically treated as a fitting parameter to be determined from comparison with measured data. The other two formulations include some functional dependence on the saturation in the fractures to accommodate changes in the hydraulic state, which is physically more meaningful than the constant-factor approach. A rigorous method for deriving R_L was developed by Liu et al. (1998), based on the hypothesis that the interface reduction factor (as well as the characteristic functions for unsaturated flow) should depend on the fraction of so-called active fractures. Active fractures are defined as those fractures or portions of fractures that actively conduct water at a given hydraulic state. Liu et al. (1998) suggested expressing the active fraction as a power function of water saturation, with the power function coefficient a positive constant depending on the properties of the corresponding fracture network.

In principle, the dual-continuum formulation for heat flow should follow a methodology similar to that developed for liquid and gas flows. For example, in correspondence to Eq. [1], the heat transferred from the matrix to the cooler liquid water flowing in fractures should be expressed as

[2]

where is thermal conductivity and R_H,L is the interface reduction factor for heat exchange between the matrix and the water in the fractures. The term R_H,L may be different from R_L in certain cases due to fracture coatings, which could have a different impact on liquid transfer vs. heat transfer across the interface. Also, under above-boiling conditions, vapor films may evolve between the flowing water and the fracture wall, potentially changing the heat-transfer conditions. One may assume, however, that the different interface-reduction formulations developed and applied to liquid and gas flows should also provide reasonable choices for heat exchange. The best-fit coefficient ß is assumed to be identical for liquid transfer and heat transfer.

While heat is transferred from the rock matrix to the liquid phase in the fractures according to Eq. [2], there is a simultaneous heat transfer between the rock matrix and the gas phase in the fractures, which can be written as

[3]

Here, R_H,G is the interface reduction factor for heat exchange between the matrix and the gas in the fractures. Remember that for finger-type flow of cool water into an initially dry fractured rock system, the water temperature in the fracture continuum may be much cooler than the temperature in the adjacent matrix, while the gas temperature remains fairly well equilibrated with the matrix temperature. Thus, rigorously applying Eq. [2] and [3] in a numerical simulation would involve allowing for a disequilibrium between the phase temperatures in each computational grid block, a condition that cannot be captured by the local thermal equilibrium assumption of most multiphase, multicomponent simulators.

This conceptual difficulty explains why the necessity of reducing the fracture–matrix interface area in response to the hydraulic state of fracture flow is largely ignored when considering nonisothermal conditions. For example, Ho (1997) and Doughty (1999) conducted systematic evaluations of various fracture–matrix interface reduction approaches for moisture and gas transport, but both opted to assume a full geometric interface area for heat transfer between the rock matrix and the fractures. The latter is equivalent to assuming that the water and gas phases are well mixed and uniformly distributed across the geometric fracture area, and that the two phases are always in local thermal equilibrium. Doughty (1999) acknowledges that this simplified treatment was chosen not because it was considered physically rigorous, but "due to the complexity of the problem." Green et al. (2003) and Birkholzer et al. (2004) used constant-factor reductions for heat-transfer problems, but did not indicate whether the reduction was deemed necessary because of flow channeling in the fractures. Other factors—such as fracture coatings, vapor films, or initial overestimation of the geometric interface—would also require reduction in heat transfer between fractures and matrix, but would not cause a local disequilibrium between gas and liquid phases.

In this study, we developed an approximate solution for interface reduction in heat transfer problems that allows the use of local-equilibrium simulators. The solution is based on the assessment that an accurate representation of the gas-phase temperature in the fractures is not necessary. Therefore, it is possible to use Eq. [2] for the heat transfer between the matrix and the liquid phase, while ignoring the heat transfer between the matrix and the gas phase in Eq. [3]. Both phases in the fractures are then represented by the liquid-phase temperature, thereby avoiding the necessity of simulating disequilibrium phase temperatures. While this approximation accurately accounts for the limited heat exchange between the rock matrix and the liquid water in the fractures, it misrepresents the gas-phase temperature in the fracture. This misrepresentation has a minor impact on the TH conditions in the fractured rock because of (i) the small relative error in energy storage, and (ii) the small relative error in interphase heat transfer (see Appendix). The proposed method allows using any of the interface reduction formulations developed for liquid transfer processes.

	THE CNWRA HEATER EXPERIMENT

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
The CNWRA experiment (Green and Prikryl, 1998, 1999; Green et al., 2003) was conducted to assess the TH conditions around a heat source in a fractured porous formation with gravity-driven water flux from above. The purpose of the experiment was to better understand the future situation near emplacement drifts at Yucca Mountain during the heating phase of the repository. A 1.2 by 0.6 by 1.2 m test cell was assembled using solid rectangular concrete blocks (0.05 by 0.6 by 0.05 m; see Fig. 3a ) with measured TH property values. The spaces between the blocks represented two sets of horizontal and vertical fractures. The assembly of concrete blocks formed a "fracture–matrix" system with a regularly structured, permeable fracture network next to low-permeability, porous "matrix blocks."

View larger version (55K): [in this window] [in a new window]   Fig. 3. Schematic of (a) Center for Nuclear Waste Regulatory Analyses heater experiment and (b) simulation domain with finite volume grid. The simulation domain—reduced to one-quarter of the laboratory test cell—is rotated around the z-axis to allow for better visualization. The local coordinate system of the simulation domain is located at the top of the test cell in the center of the liquid-release cylinder. The y-axis is aligned with the drift axis.

Three different thermal experiments were conducted that essentially differed by the duration of heating and water release. The second experiment is evaluated here, which started with 5 d of heating without water release, followed by 125 d of heating and infiltration, and completed by 10 d of infiltration during linear ramp-down of the heat source (Green and Prikryl, 1999; Green et al., 2003). We refer to the first 130 d of the experiment as the heating phase and the final 10 d of the experiment as the cooling phase.

During the experiment, temperature was recorded with 100 thermocouples. Within the fractured medium, a maximum temperature of 202°C was reported at a distance of 0.0125 m into a fracture immediately below the heater (Green and Prikryl, 1999). Above the heater, a maximum temperature of 175°C was measured at the drift crown. Thermocouples in the central fracture indicated the location of the boiling point isotherm (i.e., 100°C) at 10 cm above the drift crown at early stages (after 10 d of heating and 5 d of infiltration). With continuing infiltration, the boiling isotherm migrated downward to roughly 5 cm above the drift crown.

The presence of precipitates deposited in the borehole during the experiment indicated that water had advanced toward the heater and dripped into the drift despite the above-boiling temperatures in the adjacent matrix. As Green and Prikryl (1999) pointed out, the precipitates (which included aphthitalite and barite as well as calcite) accumulated during the heating phase of the experiment, indicating that both vaporization and capillary barriers had failed. Green and Prikryl (1999) concluded from their experiment that the presence of an above-boiling dry-out zone in the fractured rock did not preclude the dripping of water into the heated drift.

	SIMULATION OF THE CNWRA EXPERIMENT

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
In this study, simulation of the CNWRA experiment was based on a general modeling concept similar to the thermal seepage model, which was specifically designed to examine whether or not drift seepage could occur at Yucca Mountain during the above-boiling phase of the repository, while considering the possibility of strongly channelized flow in the fractures. Details on the model approach are given in Birkholzer et al. (2004) and Bechtel SAIC Company (2004a). Using a similar general approach for simulation of the CNWRA experiment allowed us to evaluate the validity of the modeling framework for a test case in which the barrier capabilities of heated fractured rock are challenged by significant water drainage from above. In other words, we could test whether the simulation approach used for predicting TH conditions at Yucca Mountain is capable of explaining the observation of seepage in the CNWRA experiment.

As mentioned above, one of the key features of the thermal seepage model is the dual-continuum approach. Thus, in our simulation study, the CNWRA test cell was represented as a dual-continuum system, with a matrix continuum representing the concrete blocks and a fracture continuum representing the gaps between them. At Yucca Mountain, adequacy of the continuum concept is based on the considerable density of the natural fracture network, which is well connected on a scale much smaller than a typical emplacement drift, i.e., the geometry of interest for seepage. In the CNWRA experiment, the fracture geometry was rather artificial, with continuous fractures and equal fracture spacing. In general, these characteristics are well suited for a continuum representation. On the other hand, because of the small distance between water release and borehole as well as the small borehole diameter compared with the fracture spacing, the experiment results may be dominated by a few vertical and horizontal fractures, therefore suggesting that a discrete representation of these dominant features may be equally appropriate.

The main focus of our study was to understand the importance of fracture–matrix heat transfer concepts in evaluating heat transfer problems in fractured porous formations. We furthermore sought to demonstrate that an established modeling methodology, based on a dual-continuum representation of the fractured porous rock, is capable of reproducing the main findings of the laboratory experiment. It was not the aim of this study to simulate the CNWRA experiment in all of its qualitative and quantitative detail. We therefore employed a dual-continuum model rather than a discrete fracture model, aware of the limitations of this modeling approach. Since the global flow and heat transport processes are important in both continua, a dual-porosity formulation (where matrix-to-matrix flow and heat transport would be neglected) was not feasible. It follows that the interaction between fracture and matrix continua needed to be described by the quasi-steady exchange terms introduced in Eq. [1] through [3]. Similar continuum representations of the CNWRA experiment were employed in the complementary modeling studies of Green and Prikryl (1999) and Green et al. (2003).

Several dual-continuum simulation cases were conducted in our study, applying different formulations for the fracture–matrix interface area reduction to better understand the related sensitivities. For comparison, we also performed one simulation run where the full geometric interface area was available for fracture–matrix heat transfer. The simulations were conducted using the multiphase, multicomponent simulator TOUGH2 (Pruess et al., 1999), with some modifications made to incorporate interface reduction formulations for heat exchange between the fractures and matrix. Since TOUGH2 assumes local thermodynamic equilibrium between the liquid and gas phase in each computational cell, we used the approximate treatment of gas-phase temperature proposed in the section "Dual-Continuum Models and Fracture–Matrix Interaction Approaches."

Model Setup, Boundary Conditions, and Thermal–Hydrological Properties
The CNWRA experiment was simulated in a three-dimensional model domain. Symmetry along two vertical planes justified a model domain reduced to one-quarter of the laboratory test cell. The one-quarter model and the numerical discretization used for the dual-continuum grid are shown in Fig. 3b. The grid consisted of 7800 grid elements.

Most of the test conditions and properties used in our simulation model were directly based on the literature describing the experiment and related analyses (Green and Prikryl, 1998, 1999; Green et al., 2003). Initially, the test cell was assumed to be dry at a gas pressure of 0.1 MPa and a temperature of 20°C. The transient simulation ran during the 130 d of heating and 10 d of cooling in the experiment. Conductive heat loss was considered by explicitly modeling the insulation layer around the test cell. To allow for gas release, the boundary elements were assigned a small gas permeability. Liquid flow across the boundaries was not permitted, consistent with Green et al. (2003).

The TH properties assigned to the fracture and matrix continua are given in Table 1. Uniform properties were used throughout the test cell except for the drift. All relevant material properties for the concrete matrix were taken from measurements (Green et al., 2003), with the exception of matrix permeability, which was reduced by one order of magnitude compared with the measured value. Green et al. (2003) suggested that such a reduction in matrix permeability could better reproduce observed matrix saturation patterns. Since measured data were not available, fracture properties were mostly estimated, based on (i) the known fracture network geometry, (ii) estimates of the average geometric fracture aperture (0.4 mm as reported in Green and Prikryl [1998]), and (iii) data from studies under similar settings. Overall, the properties listed in Table 1 are similar to those used in Green et al. (2003).

View this table: [in this window] [in a new window]   Table 1. Hydrogeological and thermal input values for simulation model.

Dual-continuum models require information on the respective volume fractions of fractures and matrix blocks in the domain, the geometric interface area between fractures and matrix blocks, and the average dimension of a typical matrix block. For the regular geometry in the CNWRA experiment, these properties were directly derived from geometric considerations. As mentioned above, a best-fit coefficient is needed in a quasi-steady transfer model because the local-scale flow and heat transport processes in the matrix (e.g., the steep gradients at the fracture–matrix interface) are not captured in detail. Ideally, these coefficients should be based on model calibration; however, since the literature on the CNWRA experiment does not provide sufficient information for a comparative analysis between the model results and measured data, we selected a best-fit coefficient based on literature recommendations for a two-dimensional fracture network with slab-like matrix blocks (e.g., Warren and Root, 1963; Ho, 1997).

Model Results
Temperature Conditions before Water Release
All simulation cases using different liquid and heat transfer formulations had identical temperature fields during the first 5 d of heating (no water release), since the same TH properties were used. Also, since the test cell was initially dry, no flow channeling in fractures and no related reduction in fracture–matrix interaction occurred. The simulated fracture and matrix temperatures were in local equilibrium. Matrix temperature contours on the three-dimensional quarter-model are shown in Fig. 4 for the situation at 5 d of heating, just before the onset of water release. The test cell had heated up considerably; the simulated temperature just above the crown of the drift in the center of the cell (x = 0 m, y = 0 m) was 170°C, which compares well with the 175°C measured in the experiment (Green and Prikryl, 1998). The boiling-point isotherm was located between 8 and 11 cm above the drift wall, depending on the axial location.

View larger version (25K): [in this window] [in a new window]   Fig. 4. Temperature contours after 5 d of heating, just before the onset of water release. Black contour line shows 100°C isotherm.

The temporal evolution of simulated fracture and matrix temperatures, as well as fracture liquid saturation, are presented in Fig. 5a and 5b, for a location in the fractured medium just above the drift crown in the center of the test cell. The distance of the boiling isotherm from the drift crown as a function of time is given in Fig. 5c. Note that the drift crown is about 0.52 m below the top of the test cell, where water infiltration began after 5 d of heating. It is obvious from the temperature plot that the water release of 1000 mL d^–1 had a significant impact on the thermal conditions in the test cell. All simulation cases exhibited strong temperature drops at the start of the infiltration period, the characteristics of which depended on the considered scenario.

View larger version (25K): [in this window] [in a new window]   Fig. 5. (a) Temperature and (b) saturation evolution just above the crown of the drift in the center of the test cell (x = 0 m, y = 0 m); (c) gives the evolution of the distance of the boiling isotherm above the drift crown, measured in the center of the test cell (x = 0 m, y = 0 m). The time axis scale is identical in all three graphs; M = matrix, F = fracture.

In Case A, the matrix temperature at the drift crown continued to rise for a few days before the impact of infiltration became evident. With time, temperature decreased steadily to 108°C at the end of the heating phase. More sudden responses in the matrix temperature occurred in Cases 1a and 1b, with significant temperature drops soon after the onset of water infiltration. As will be shown below, these accelerated responses were caused by the fast movement of the liquid front down to the drift vicinity, a result of reduced imbibition. The sudden responses were followed by a relatively slow temperature decrease as time progressed; at the end of the heating phase, the matrix temperatures in Case 1a (113°C) and Case 1b (119°C) was actually higher than in Case A.

The simulated fracture temperatures closely followed the matrix temperatures in Cases A and 1a, with some minor deviations at the end of the heating phase caused by the impact of nearby heat pipes. This is in strong contrast to the fracture temperature evolution observed in Case 1b, i.e., the simulation with reduced interface area for liquid and heat. The fracture temperature at the drift crown dropped to the boiling point immediately after water release and remained at 100°C during the entire heating phase. This drop created a strong disequilibrium between the matrix and the fracture temperatures, with matrix temperatures up to 60°C higher at early heating stages. The immediate drop in the fracture temperature of Case 1b corresponds with the early arrival of liquid water at the drift crown (see Fig. 5b). Thus, the simulation with reduced interface area for heat transfer allowed a situation in which fast-flowing water could penetrate the above-boiling fractured rock and reach the drift crown, even while the temperature in the adjacent matrix blocks was as high as 160°C. At early times, the amount of water arriving at the drift crown was quite small, with saturation just above residual. Later, the fracture saturation at the drift crown built up, indicating that the impact of vaporization was reduced with decreasing matrix temperatures. The maximum saturation of 0.1 was reached between 70 and 130 d of heating. Seepage of water into the drift is possible under such conditions, and ultimately depends on the effectiveness of capillary forces at the rock–drift interface. In contrast, Cases A and 1a predict a fully effective vaporization barrier, with zero fracture saturation at the drift crown during the entire heating phase. Thus, the simulation cases assuming a full geometric interface area for heat transfer cannot explain why dripping of water was observed during the heating phase of the experiment.

Note that Case 1b—with reduced thermal interaction between fracture and matrix—features the highest matrix temperatures during most of the heating phase, except for the early stages. Less cooling of the matrix blocks occurred, since the heat exchange with the cool water in the adjacent fractures was diminished. This effect can also be seen in the temporal evolution of the boiling point isotherm in Fig. 5c. All simulation cases exhibited a fast early response, as the boiling zone receded to less than half of its initial extent. At later stages, however, the boiling isotherm remained farther away from the drift crown in Case 1b (4–5 cm at the end of the heating phase) compared with Cases A and 1a (2–3 cm at the end of the heating phase).

All simulation cases showed a strong temperature decrease during the 10-d cooling phase, down to less than 40°C at the end of the experiment (140 d). As temperature decreased below boiling, the differences between the simulation cases virtually disappeared. Water saturation at the drift crown built up immediately to maximum values of about 0.25. Thus, during the cooling phase, seepage of water into the drift is likely in all three cases; however, for analysis of the above-boiling period, we are only interested in the 125 d of simultaneous heating and water release.

We now compare the three-dimensional flow patterns of water infiltrating down the fracture network toward the heated drift. Fracture saturation contours at 10, 50, and 130 d after the onset of heating (i.e., after 5, 45, and 125 d of water release) are depicted in Fig. 6 for the upper part of the three-dimensional model domain. For comparison, the plots also include the location of the boiling isotherm determined from simulated matrix temperatures. All cases exhibited fast, focused, gravity-driven downward flow in response to water release at the top of the test cell. For example, after 5 d of infiltration, the liquid front had traveled most of the vertical distance between the release point and the drift. The colder water flowing in the fractures heated up on its way down; simultaneously, matrix temperature decreased in response to water arrival. As water encountered above-boiling temperatures, the liquid front slowed down. Vaporization and subsequent condensation led to lateral flow diversion and drainage of water toward the bottom boundary. (Note that Green et al. [2003] reported flow diversion around the drift in a zone 15 cm wide, as concluded from mineral precipitation patterns observed after disassembling the test cell. This agrees favorably with the saturation distributions shown in Fig. 6.) Saturations below the drift remained low because the flow region was sheltered from the downward infiltration.

View larger version (62K): [in this window] [in a new window]   Fig. 6. Fracture saturation contours after 10, 50, and 130 d of heating for simulation runs (a) Case A, (b) Case 1a, and (c) Case 1b. Black contour line shows 100°C isotherm.

Outside of the boiling region, Cases 1a and 1b have virtually identical saturation patterns. The reduction in heat transfer mostly affected the effectiveness of vaporization, because boiling requires considerably more energy than merely heating of water. There are, however, obvious differences compared with Case A regarding the flow patterns observed outside of the boiling zone. Since Case A assumed a full interface area for fracture–matrix liquid transfer, more water imbibed into the adjacent matrix compared with the other cases, which retarded the early-time movement of the liquid front (see Fig. 6a, 6b, and 6c at 10 d of heating). As a result, Case A features the least extended liquid front and the largest boiling region at 10 d of heating, when water release had just started. On the other hand, the same simulation case exhibits the smallest boiling region of all cases at the end of heating, because the interaction between the cool fracture water and the matrix was more effective. The temperature even dropped below the boiling point near the unheated end of the drift, and water accumulated at the drift wall at 130 d of heating (which opens the possibility of seepage). However, Case A is not a likely scenario because all temperature measurements in the CNWRA experiment suggest that the entire drift crown remained above boiling during the heating phase of the experiment (Green and Prikryl, 1998, 1999).

Model Results for Other Sensitivity Cases
We evaluated additional sensitivity cases using a range of interface reduction formulations as listed in Table 2. These include constant-factor reductions using values of 0.1, 0.01, and 0.002, as well as cases where the interface reduction was calculated from the active fracture model. The values used, 0.3 and 0.569, are based on calibration to Yucca Mountain field measurements (Liu et al., 2003; Bechtel SAIC Company, 2004b). The resulting interface reduction factors of all cases are plotted as a function of saturation in Fig. 7 . Some cases are associated with a significant reduction in the geometric interface, most notably the constant-factor case with 0.002 and the active fracture case with = 0.569. Similar to the model results for saturation-dependent interface reduction, above, we conducted two types of simulations for each case, one with interface reduction for liquid only, the second with interface reduction for liquid and heat.

View this table: [in this window] [in a new window]   Table 2. Sensitivity cases with different interface reductions and maximum fracture saturation values in grid elements adjacent to drift wall (upper half), simulated during the heating phase of the Center for Nuclear Waste Regulatory Analyses experiment.

View larger version (25K): [in this window] [in a new window]   Fig. 7. Fracture–matrix interface reduction factor given as a function of saturation. For cases using the active fracture model, the saturation of active fractures, SLe,a, is shown.

In contrast, all cases that allowed for heat transfer reduction exhibited water buildup at the drift crown, with maximum fracture saturations ranging from 0.09 in Case 1b to 0.31 in Case 6b. Locally, these values would be even higher if small-scale heterogeneities in fracture permeability had been accounted for in the model. (As shown in Birkholzer et al. [2004], small-scale heterogeneity can create local saturation maxima at the drift wall, thereby increasing the seepage potential). We may conclude that most, if not all, sensitivity cases with interface reduction for heat transfer arrived at near-drift conditions that make thermal seepage at above-boiling conditions likely. Cases with strong interface reduction for fracture–matrix heat transfer resulted in the highest fracture saturation values at the drift wall and, thus, the highest the potential for seepage.

It follows that the sensitivity cases featuring interface reduction for fracture–matrix heat transfer are generally better suited to represent the CNWRA test results than those assuming that heat is conducted across the full geometric area. But which of the interface reduction formulations for heat transfer is in best agreement with the experimental data? This is a difficult question to answer, for the following reasons. First, the literature on the CNWRA experiment does not provide sufficient information to allow a comprehensive comparison between model results and measured data. Second, all our sensitivity cases—ranging from a moderate to a rather extreme reduction in heat transfer—show water arrival at the drift crown despite above-boiling matrix temperatures. Differences between these cases would be more obvious in a situation where the boiling isotherm was farther away from the drift than the few-centimeters distance observed in the CNWRA experiment. Presumably, with a larger boiling region, water would not reach the drift crown in the moderate sensitivity cases, while extreme interface formulations would still allow breaching of the vaporization barrier.

Implications for Yucca Mountain Models
Our modeling analysis of the CNWRA experiment has demonstrated that the experimental observation of thermal seepage into a heated drift with above-boiling temperatures in the adjacent rock can be reproduced in the numerical simulations. All sensitivity cases that accounted for channelized flow in fractures and the related reduction in fracture–matrix heat transfer were capable of predicting TH conditions favorable for thermal seepage to occur under such conditions. This result provides confidence that the conceptual model discussed in the section "Dual-Continuum Models and Fracture–Matrix Interaction Approaches," using an approximate treatment of thermal fracture–matrix interface reduction, is adequate.

A similar conceptual approach for reducing heat transfer between fractures and matrix was chosen for the predictive analysis of thermal seepage during the above-boiling phase of the Yucca Mountain repository. The thermal seepage model developed for this task used a constant-reduction factor of about 0.002 for thermal interface reduction (Bechtel SAIC Company, 2004a), similar to Case 4b studied here. (Instead of modifying the actual interface area, the fracture–matrix heat transfer reduction was realized by an equivalent adjustment of the thermal conductivity used for the heat exchange between the fractures and the matrix). This rather extreme interface reduction probably overestimates the impact of flow channeling on fracture–matrix interaction and thus exaggerates the potential of water arrival at the drift crown. Despite this extreme modeling choice, water was not able to penetrate far into the above-boiling regions and no seepage occurred in the Yucca Mountain simulations.

We may conclude that the thermal seepage observations in the CNWRA experiment are not representative of the future TH conditions expected at Yucca Mountain. With respect to the possibility of water arrival at the drift crown, the most important differences between the CNWRA experiment and the Yucca Mountain conditions are related to (i) the size of the heated rock volume, in particular the extent of the boiling region, and (ii) the intensity and the temporal evolution of downward flow above the heated drifts. In a previous study, we analyzed the predicted TH conditions relevant for thermal seepage at Yucca Mountain during the above-boiling period (Birkholzer et al., 2003). On average, the boiling region above a typical emplacement drift at Yucca Mountain is expected to cover >4 m of fractured rock during the first 500 yr after emplacement. Conversely, in the CNWRA experiment, water released at the top gave rise to significant cooling of the test cell, and the location of the boiling isotherm receded to a few centimeters above the drift crown. The smaller the boiling zone, the larger the possibility of channelized finger flow reaching the drift.

To evaluate the intensity and the temporal evolution of downward flow above the heated drifts, we compared the characteristics of the downward flow processes in the experiment to the Yucca Mountain setting. Since the average infiltration in the unsaturated rock is very small (5–10 mm yr^–1), the main source of downward flow toward heated drifts at Yucca Mountain would be from condensation of matrix pore water mobilized in the boiling region. It follows that the intensity and evolution of thermal flux perturbation is directly linked to the intensity and evolution of heat. Accordingly, Birkholzer et al. (2003) reported that the condensation-driven downward flux toward heated drifts at Yucca Mountain would be largest during the first few hundred years after emplacement, corresponding to the period when the heat output of the radioactive waste would be highest. During this time, the boiling zone in the fractured rock would be large and would effectively prevent water arrival at the drift, although the maximum downward fluxes of condensate (80 mm yr^–1 for average conditions) would exceed the natural infiltration by more than an order of magnitude. Later, when less heat emanates from the waste packages and the boiling zone becomes smaller, the thermal flux perturbation would decrease accordingly, so that water arrival at the drift is also unlikely. This temporal correlation between flux intensity and extent of the boiling zone at Yucca Mountain was shown to be very important in preventing thermal seepage at above-boiling conditions (Birkholzer et al., 2003).

Conversely, in the CNWRA experiment, downward flow was initiated by forced water release from the top of the test cell, maintained at a constant rate of 1000 mL d^–1 independent of thermal conditions. Using the footprint of the drift (0.15 by 0.6 m) as a reasonable reference scale, this water release rate converts to a Darcy flux of >4000 mm yr^–1, which is much larger than the condensate-driven downward flux expected for average infiltration conditions at Yucca Mountain (Birkholzer et al., 2003). Even if an emplacement drift was located in a high-percolation area and small-scale flow channeling was considered in the fracture continuum, the maximum downward fluxes expected at Yucca Mountain would still be smaller than those in the CNWRA experiment. Birkholzer et al. (2004), for example, reported maximum fluxes of 350 mm yr^–1 for a TH simulation that accounted for small-scale heterogeneity in fracture permeability and assumed a percolation flux 10 times higher than average.

We may conclude that the test conditions in the CNWRA experiment were extremely favorable for thermal seepage to occur despite above-boiling rock temperatures, compared with the future situation at Yucca Mountain. In the CNWRA experiment, forced water release led to intense gravity-driven downward flow that could easily penetrate through the small boiling region around the heated drift and lead to water accumulation at the drift crown. At Yucca Mountain, the maximum gravity-driven fluxes are expected to be much smaller than in the experiment, while the boiling region is up to 80 times larger.

	SUMMARY AND CONCLUSIONS

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
We conducted a sensitivity analysis studying the impact of fracture–matrix mass and heat transfer in fractured porous rock heated to above-boiling temperatures. Our motivation was to gain a better understanding of processes important for thermal seepage at the proposed geologic repository for radioactive waste at Yucca Mountain, NV, during the period of above-boiling rock temperatures. The sensitivity study was performed using a laboratory heater test (CNWRA experiment) as a simulation example. It was shown in this heater test that fast downward flow in fractures could give rise to thermal seepage despite above-boiling rock temperatures, thereby challenging earlier model predictions for Yucca Mountain in which thermal seepage was essentially ruled out for such temperature conditions.

Based on the modeling framework developed for the predictive evaluation of thermal seepage at Yucca Mountain, the CNWRA experiment was simulated with a dual-continuum approach, which treats the fractured rock as separate, interacting continua for fractures and matrix blocks, respectively. To account for channelized, finger-type flow in the fractures, the degree of fracture–matrix interaction was reduced by applying different interface reduction approaches available from the literature, specifically multiplication of the full geometric interface with a constant factor, multiplication with fracture liquid saturation, and multiplication with factors derived from the active fracture model. We considered two main categories: (i) interface reduction for liquid transfer, but full geometric interface for heat transfer, and (ii) interface reduction for liquid and heat transfer. For comparison, we also studied a case using the full geometric interface without any reduction. The following two main conclusions can be drawn from this analysis.

1. Our simulation results stress the importance of understanding and adequately simulating the degree of heat transfer between the matrix and the flowing water in fractures. Only those sensitivity cases that considered interface reduction for liquid and heat transfer were capable of reproducing the CNWRA test results. These cases—featuring moderate to rather extreme interface reduction for heat transfer—exhibited early and consistent arrival of water at the drift crown, even though the matrix temperatures showed no clear indication of such events and remained above boiling throughout the heating phase. Seepage is possible under these conditions, which is consistent with the experimental findings. In contrast, if interface reduction was not considered at all or was only applied to liquid exchange processes, the saturation front and the boiling isotherm coincided and water was prevented from arriving at the drift crown by a fully efficient vaporization barrier.
   
2. All simulation cases have indicated that the CNWRA experiment was operated at conditions favorable for thermal seepage. Most importantly, these conditions were (i) the rather small boiling region near the drift, and (ii) the strong gravity-driven downward flow in response to forced water release from the top. These conditions are not representative of the future TH situation expected at Yucca Mountain, where the boiling region is much larger and where the downward flow characteristics are different in magnitude and temporal evolution.
   

Attempts to identify the interface-reduction formulation best suited to describe the reduced heat transfer in the CNWRA experiment failed. This is in part a result of the test geometry because the region of above-boiling rock temperature was too small to bring out clear differences between the simulation cases. Thus, future work should be directed at conducting and interpreting large-scale in situ experiments, ideally with an above-boiling rock region of a few meters extent. Such a test should be operated at different levels of water release to create conditions more or less favorable for thermal seepage.

	APPENDIX

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION BASIC PROCESSES AND DUAL... THE CNWRA HEATER EXPERIMENT SIMULATION OF THE CNWRA... SUMMARY AND CONCLUSIONS APPENDIX REFERENCES

 
The approximate approach for thermal fracture–matrix interaction introduced in the section "Dual-Continuum Models and Fracture–Matrix Interaction Approaches" ignores the local disequilibrium between liquid- and gas-phase temperatures that can evolve as a result of channelized water flow in the fracture continuum. To allow application of equilibrium models, the liquid-phase temperature in the fractures is used to describe the thermal state in both phases. This approximate treatment misrepresents the gas-phase temperature, introducing an error in the numerical simulation. The impact of this error is evaluated below for the TH conditions observed in the CNWRA experiment. The situation is similar to the schematic in Fig. 2, where liquid fingers penetrate into an above-boiling fracture–rock system. The liquid is at a temperature close to the boiling point of water at prevailing pressure (T_F,L 100°C), while the matrix temperature is significantly higher (compare Fig. 5a). The gas phase in the fractures is similar to the matrix temperature (T_F,G T_M), but in the approximate scheme is described by the local temperature of the liquid phase. Thus, the error in gas-phase temperature equals the difference between T_M and the boiling temperature of water.

As the first criterion, we calculate the relative error in energy storage. The error in the gas-phase temperature leads to an energy discrepancy in the fracture. This discrepancy should be small compared with the heat stored in the other phases present in the system, i.e., the solid phase and the water phase.

For a given fracture volume V_F and given gas saturation (1 – S_L), the energy discrepancy E_G in the gas phase can be calculated as

[A1]

where _G is gas density and C_G is gas specific heat. The fracture volume is 4 x 10^–4 m³ for a 1-m² fracture with an aperture of 0.4 mm (the aperture value estimated for the CNWRA experiment). We may assume that the gas phase in the fracture consists entirely of vapor and insert the thermophysical properties of vapor at 100°C and 0.1 MPa (_G 0.6 kg m^–3; C_G 2.02 kJ kg-K^–1). Reasonable values of liquid saturation and matrix temperature for the CNWRA experiment can be extracted from Fig. 5. With a liquid saturation of 0.05 and a matrix temperature of 160°C (representative of the test conditions just after water release), the resulting energy discrepancy is 0.03 kJ. (Similar results would be obtained when the gas phase comprises an air–vapor mix or pure air. For example, using _G 0.9 kg m^–3 and C_G 1.0 kJ kg-K^–1 of pure air at 370 K, the resulting energy discrepancy would be 0.02 kJ.)

First compare this energy discrepancy with the energy stored in the solid phase (i.e., the matrix). We can calculate the matrix volume V_M associated with the fracture volume V_F from the fracture volume fraction (0.016 in Table 1), giving V_M = (1 – 0.016)V_F/0.016 = 0.0246 m³. If we ignore the energy contributions from gas or liquid phases stored in matrix pores, the energy E_M contained in V_M is