Abstract

The Radioactive Waste Management Directorate of the UK Nuclear Decommissioning Authority has prepared a generic disposal system safety case (DSSC) that covers a range of possible host rock environments. In many of the waste packages considered in the DSSC, the formation of gases by chemical and microbial processes is likely to occur. In order to demonstrate safety, it is necessary to understand the rates at which the gases are generated and their subsequent migration from a disposal facility after closure. This paper is concerned with modelling gas migration through a fractured higher strength host rock. A first set of simulations compares alternative approaches to modelling gas migration through a fractured rock. The approaches differ in their representation of the interaction between the fractures and the rock matrix. As expected, the gross features of many of the simulations are very similar, with a single continuum approach in which the porosity is set equal to either the fracture porosity or the matrix porosity providing the bounding cases. Gas migration is slower for those simulations where the gas can access more of the rock matrix. A final simulation, with a heterogeneous permeability field, is compared with the other simulations, again showing a very similar evolution.

Introduction

The Radioactive Waste Management Directorate (RWMD) of the UK Nuclear Decommissioning Authority (NDA) has prepared a generic disposal system safety case (DSSC) (Nuclear Decommissioning Authority, 2010a) that considers the safety of the construction and operation of a geological disposal facility (GDF) for three relevant host rock environments (higher strength rocks, lower strength sedimentary rocks and evaporites). One of the factors affecting safety is gas generation; chemical and microbial processes will lead to the generation of gases in many of the waste packages that will be emplaced within the GDF.

In a higher strength rock, as is considered in this paper, the gases could migrate through discontinuities in the rock and make their way back to the biosphere; a concern is that radioactive gases (e.g. 14C-labelled methane) could migrate to the biosphere if a sufficient volume of carrier gas (i.e. hydrogen) were generated to form a gas pathway. In contrast, both lower strength sedimentary rocks and evaporites have high ‘gas entry pressures’ and would obstruct the gases; in these environments, the gases could accumulate within the GDF and therefore the pressure would build up, which could affect the integrity of the surrounding host rock and/or force radioactive species in solution away from the GDF. Thus, it may be necessary to show that the migration of the gases is well understood in order to demonstrate the safety of the GDF.

This paper is concerned with modelling gas migration through a higher strength host rock. For such a host rock, any fluid movement will be predominantly through discontinuities (e.g. fractures or joints) in the rock (Long et al., 1996). Several conceptual models have been developed for describing flow and transport through a fractured rock (Lichtner, 2000; Berkowitz, 2002).

In order to simulate two-phase (i.e. gas and groundwater) flow through a fractured rock, it may be necessary to model the fractures, the rock matrix between the fractures, and their interaction. The approaches that are available for modelling this problem differ in the way in which they treat the porosity and permeability of the fractures and the rock matrix. The most realistic would be an explicit discrete fracture model. However, the computational intensity of this approach, as well as the difficulty in acquiring detailed knowledge of the fractures at a given site, mean that a continuum model is a more practicable approach for simulating two-phase flow.

The single continuum method represents the fractures and the rock matrix by a single, effective continuum. The dual continuum method includes dual porosity, dual permeability, and the more general ‘multiple interacting continua’ (MINC) models. The main features of these models can be summarized as follows:

  1. In the dual porosity model, the domain is composed of low-permeability matrix blocks embedded in a connected network of fractures. Global flow and transport occur only through the network of fractures, which is represented by an effective continuum. The matrix blocks are treated as sources/sinks in the effective continuum.

  2. In the dual permeability model, global flow and transport occur through the matrix blocks as well as through the network of fractures.

  3. The MINC model is an extension of the dual porosity model. It resolves the gradients that drive the flow between the fractures and matrix blocks by discretizing the matrix blocks into a nested sequence of volume elements. In response to an imposed disturbance in the fracture network (e.g. a pressure drop or a change of phase composition), fluid can migrate in the matrix blocks either outward towards the fractures or inward away from the fractures.

This paper compares predictions of gas migration through a higher strength host rock obtained using several of these alternative continuum models.

Methodology

West Cumbria is used to develop our understanding of gas migration through a higher strength rock, because the geology and hydrogeology were characterized in studies during the 1990s and these data are readily available.

Model domain

A two-dimensional vertical cross-section is modelled because this is considered to be sufficient to compare the different approaches available for modelling the problem. The geological structure of the model is shown in Fig. 1. The numerical grid is aligned with the x and z directions (indicated in Fig. 1), except in the region above the GDF where the grid is aligned with the Basal Part of the North Head Member (Baker et al., 1997).

Hydraulic properties

A comprehensive list of the hydraulic properties of the various hydrogeological units can be found in reports from the Nirex 97 assessment (Baker et al., 1997). However, site-specific forms for the capillary pressure and relative permeability functions that characterize two-phase flow were not determined. The functions used in this paper are realistic forms appropriate to real rocks, but are not site specific. The functions are based on the Brooks and Corey formulation (Brooks and Corey, 1964). For grid blocks where the groundwater is considered to flow predominantly in the fractures, the gas entry pressure is set to a small value. For grid blocks in which the groundwater is considered to flow mainly through the rock matrix, however, the gas entry pressure is assumed to be correlated with the permeability (Ingram et al., 1997). Although diffusion is unlikely to be a significant process, it is included in the model for completeness. The diffusive flux depends on the tortuosity, for which values are given in Baker et al. (1997).

Gas generation

The gas generation rates used in this paper are taken from the calculations for a recent assessment (Nuclear Decommissioning Authority, 2010b). The calculations assume that the waste containers will not isolate the wastes, which will be distributed homogeneously in the GDF.

Hydrogen is the gas produced in bulk. Figure 2 shows that immediately after closure of the disposal facility, which is assumed to happen at 2150, there is a short-term peak in the hydrogen generation rate from the corrosion of reactive metals (i.e. aluminium, Magnox and uranium) followed by slower generation from the corrosion of steel. For the geological setting being considered, this hydrogen generation rate is sufficient for a gas phase to form.

Methane and 14C-labelled methane are formed in trace amounts. Methane is generated from the microbial degradation of organic material, but this methanogenesis process is inhibited by the presence of nitrate and sulfate until these chemical species have been consumed almost 700 years after construction of the GDF. 14C-labelled methane is generated from corrosion of irradiated metals and graphite.

TOUGH2

The computer program TOUGH21 was used to model the migration of gas through the geosphere. All of the calculations were carried out using TOUGH2 with the EOS7R fluid property module. The EOS7R module was customized to track water, brine, hydrogen, methane and 14C-labelled methane.

Results

Five separate simulations were performed, using, respectively: (1) the single effective continuum approach, with the porosity of the fractured rocks set equal to the fracture porosity; (2) the single continuum approach, with the porosity set equal to the matrix porosity; (3) a dual porosity model; (4) a MINC model (note that this model also accounts for the process of rock-matrix diffusion); and (5) the single effective continuum approach with a heterogeneous permeability field.

Despite the different approaches, all of the simulations show a broadly similar evolution. The essential difference between the simulations is that where the gas can access more of the rock matrix, the gas migrates more slowly and so formation of the dissolved gas plume is retarded.

Single effective continuum model

In the single effective continuum model, the fractured host rock has a porosity of 1.5 × 10−5, a permeability of 1.24 × 10−17 m2 and a gas entry pressure of 102 Pa; the rock matrix is not modelled.

For the first 110 years of the simulation, the GDF is assumed to be open and ventilated. The presence of the GDF at atmospheric pressure before closure leads to a region of drawdown around the GDF. As a result nearly fresh and saline water are drawn towards the GDF.

After closure of the GDF, when the groundwater flowing into the GDF is no longer drained, the pressure of the gas in the GDF starts to build up and quickly returns to near hydrostatic values. The distribution of salinity also starts to return to its equilibrium distribution. However, possibly because of the effect of the migrating gas, the distribution of salinity does not return all the way to its equilibrium distribution by the end of the simulation (at 1000 years after construction of the GDF).

Once the gas pressure in the GDF has increased sufficiently (i.e. to a fraction of hydrostatic pressure), a free gas phase is able to migrate out of the GDF. The gas migrates upwards through the host rock and overlying rocks until it reaches the Basal Part of the North Head Member (which has a low permeability, 4.89 × 10−17 m2, and a high gas entry pressure, 1.42 × 106 Pa). A gas pocket forms, and the gas pressure rises until it overcomes the gas entry pressure of the overlying rocks. The gas phase is then able to migrate into the faulted undifferentiated St Bees Sandstone and upwards, but does not break through at the surface. Instead, all of the gas dissolves in the groundwater flowing in the near-surface aquifer (Fig. 3).

Initially, the gas phase rises due to buoyancy and the gas pathway is close to vertical. At later times, the gas migrates towards the Fleming Hall Formation. This is due to dissolved gas being transported in the groundwater (which flows from right to left in the figures).

Dissolved gas from the GDF is transported in the groundwater through the host rock and overlying rocks (Fig. 4). As the gas rises, the pressure drops and, as a result, some of the gas can come out of solution and the region of non-zero gas saturation below the Basal North Head Member is extended (Fig. 3).

As the gas phase does not break through at the ground surface (but instead dissolves in the groundwater flowing in the near-surface aquifer), the radiological risk associated with the migration of 14C-labelled methane to the biosphere is determined to be very small.

Porosity set equal to the matrix porosity

A simulation in which the porosities of fractured rocks were set equal to the matrix porosities (e.g. the porosity of the fractured host rock was increased to 1.1 × 10−2), instead of to the fracture porosities, was used to illustrate the maximum possible effect of interactions between fluids in the fractures and in the rock matrix. In other words, the interactions between fluids in the fractures and in the rock matrix were assumed to be sufficiently quick to establish a local equilibrium.

Although the gross features of this model are the same as for the previous case (porosity set equal to the fracture porosity), the details are different. The increased porosity of the host rock has a number of related effects:

  1. The gas phase migrates through the host rock at much smaller gas saturation.

  2. TOUGH2 assumes that for a gas phase to form in a grid block, the pressure associated with the dissolved gas must exceed the liquid pressure. That is, the concentration of dissolved gas in the grid block must reach a critical value, dependent on the Henry's law constant and the liquid pressure, before a gas phase will form. (In order
    Fig. 4.

    The mass fraction of dissolved gas at 1000 years after construction of the GDF for the single (effective) continuum model, with porosity set equal to the fracture porosity.

    Fig. 4.

    The mass fraction of dissolved gas at 1000 years after construction of the GDF for the single (effective) continuum model, with porosity set equal to the fracture porosity.

    Fig. 5.

    The mass fraction of dissolved gas at 1000 years after construction of the GDF for the single continuum model, with porosity set equal to the matrix porosity.

    Fig. 5.

    The mass fraction of dissolved gas at 1000 years after construction of the GDF for the single continuum model, with porosity set equal to the matrix porosity.

    to simplify the formulation of simulators for two-phase flow, it is usual to assume, as here, that the system is in local thermodynamic equilibrium). As a result, increasing the porosity delays gas migration. Compared to the previous case, less gas migrates to the Basal Part of the North Head Member at any given time.
  3. As a consequence of the previous effect, less gas migrates into the overlying faulted undifferentiated St Bees Sandstone. The formation of the dissolved gas plume is retarded, and the plume flows at greater depth.

Figure 5 shows the predicted dissolved gas plume at 1000 years after construction of the GDF. This may be compared with the dissolved gas plume from the previous case (porosity set equal to the fracture porosity), which is shown in Fig. 4.

Dual porosity and MINC models

In both the dual porosity and MINC models, the fractured host rock has a porosity of 1.5 × 10−5, a permeability of 1.24 × 10−17 m2 and a gas entry pressure of 102 Pa; the rock matrix has a porosity of 1.1 × 10−2, a permeability of 1.0 × 10−20 m2 and a gas entry pressure of 2 × 107 Pa.

It is expected that the single continuum simulations in which the porosity of the fractured rocks is set equal to either the fracture porosity or the matrix porosity will be bounding cases, as these represent the two limits for the interactions between fluids in the fractures and in the rock matrix. This was found to be true, with the results for the dual porosity and the MINC models falling between the two, although much closer to the case with the porosity of the fractured rocks set equal to the fracture porosity.

The three simulations in which the porosity of the fractured rocks is set equal to the fracture porosity all show very similar migration at early times. However, at later times the results begin to differ due to the different assumptions about the interactions between fluids in the fractures and in the rock matrix. In particular, there is a small but noticeable difference between the predictions of the dual porosity model and the MINC model due to the fact that the former assumed a three-dimensional geometry for the fracture network whereas the latter assumed a one-dimensional geometry. This meant that in the dual porosity model more of the rock matrix is closer to the fractures, and therefore more gas is able to access the rock matrix. For those simulations where the gas can access more of the rock matrix, gas migration is slower.

Heterogeneous permeability field

A simulation with a heterogeneous permeability field2 was also performed. The objective of this simulation was to demonstrate a capability to incorporate models of the spatial variability in permeability for each of the different rock units into TOUGH2 calculations, and to compare the results obtained with the results of the previous simulations. Therefore only a single realization of the permeability field was modelled. (In practice, results would be obtained from simulating gas migration through a number of realizations.) The results show a similar evolution to the previous simulations. The small differences between this simulation and the previous simulations are: (1) gas migration out of the GDF is retarded; (2) the gas pathway forms branches on length scales comparable to the grid block size (i.e. comparable to the scale of the heterogeneity); and (3) more gas goes into solution because gas migration out of the GDF is retarded, and therefore the lower plume of dissolved gas that develops towards the end of the simulation is more extensive (Fig. 6).

Conclusions

Different approaches to simulating gas migration from a GDF through a higher strength host rock have been implemented in the computer program TOUGH2. The approaches include both single continuum and dual continuum methods, and, in addition, a simulation with a spatially varying permeability field has been performed.

For the particular geological environment considered here, the results obtained show that the different approaches to modelling gas migration produce similar results on short timescales, when the gas in the fractures has had little time to interact with the rock matrix. On longer timescales, however, the interaction does result in some differences. Gas migrates more slowly for those simulations where the gas can access more of the rock matrix.

Acknowledgements

This work was funded by the Radioactive Waste Management Directorate (RWMD) of the UK Nuclear Decommissioning Authority. The authors thank Cherry Tweed and Simon Norris of RWMD for supporting this work.

1
TOUGH2 is a program for simulating coupled fluid and heat flows for multi-component, multiphase fluid mixtures in porous and fractured media (Pruess et al., 1999). It was developed by Lawrence Berkeley National Laboratory, and was qualified for use on the Yucca Mountain Project. Software qualification involved applying TOUGH2 to a variety of test problems in order to verify and validate the code.
2
Calculations to determine effective hydrogeological properties for the Nirex 97 assessment (Baker et al., 1997) were used to estimate the spatial variability of the permeability for blocks of about 10 m size. Conceptual models for each rock unit were parameterized using measurements on various length scales (i.e. measurements on core samples with scales of centimetres, tests in boreholes on scales of decametres, and field measurements on scales of kilometres). Thus, both effective (i.e. homogeneous) properties and their associated uncertainties were determined.
P1Freely available online through the publisher-supported open access option.