Introduction

Geological disposal in a multi-barrier system, comprising engineered and geological barriers, is commonly accepted as the most feasible long-term strategy for high-level radioactive waste (HLW) disposal. Consequently, this system has seen a significant amount of international research by disposal implementers, regulators, and the international research community (Selvadurai and Nguyen 1996; Ikonen 2003; NIREX 2003; Chapman 2004; Rutqvist et al. 2005a, b, 2009; Nguyen et al. 2009; Kim et al. 2011). Currently, there are no fully operational HLW disposal sites in any of the 31 countries that benefit from nuclear power (Kim et al. 2011), and disposal concepts vary throughout the world based on different geological conditions. Disposal facilities in clay, salt, and fractured crystalline rock continue to be investigated. Heat generation and high levels of radioactivity are the key criteria for classification as HLW and are therefore the main considerations in repository design (NDA 2010). In fractured rock facilities, the near-field engineered barrier aims to isolate radionuclides from the biosphere via a multi-barrier approach of waste form, canister, and bentonite buffer. The near field, i.e. the immediate vicinity of the waste disposal drifts, is subject to the most pronounced complex interactions between thermal, hydrogeological, mechanical, chemical, and biological processes. Consequently, there has been a concerted international research effort into these phenomena via experiments and coupled process modelling, e.g. the PEBS, DECOVALEX, THERESA, and FEBEX projects, as well as site-specific work led by the implementer (NIREX 1997a, b; Rutqvist et al. 2001, 2005a, b, 2009; Ikonen 2003; Löfman 2005; Nguyen et al. 2005, 2009; Anttila et al. 2008; Tsang et al. 2009; Mayor et al. 2014). The UK currently does not have a selected site. Historically, most research has focussed around the Sellafield area of Cumbria, NW England, associated with the initial investigations for an intermediate-level waste (ILW) repository in the 1980s and 1990s (NIREX 1997a). That work has provided large, high-quality datasets, which have provided encouragement to investigate the thermal effects of a HLW repository at this saturated fractured rock site.

Understanding of thermo-hydrogeological disturbances caused by radiogenic decay in fractured crystalline rock repository conditions has developed over the past few decades with much work undertaken to quantify thermal effects through both experiments and modelling, e.g. at Yucca Mountain in the USA among other sites. However, much of the modelling is site specific and has not necessarily encompassed the processes and conditions that may apply to a UK site. For example, the majority of the Yucca Mountain work focussed on the complex behaviour of boiling and condensation in the unsaturated zone for example (Birkholzer and Tsang 2000; Haukwa et al. 1999, 2003; Buscheck et al. 2003), which is not a consideration for a saturated repository concept such as Sellafield, or the Finnish and Swedish KBS-3 concept. This research led also to an improved understanding of the coupling between thermo-hydrogeological and mechanical processes especially in fractured rock masses (Rutqvist and Tsang 2003; Tsang et al. 2009). While thermo-mechanical processes can clearly have an effect on fracture permeability and may influence the permeability field for a continuum model (Hudsen et al. 2005), this study is primarily focussed on the thermo-hydrogeological processes and does not consider the mechanics. This is in part because the uncertainty in modelling regional scale repository effects, as opposed to near field effects, is largely a function of the uncertainty of up-scaled hydraulic properties rather than the added uncertainty of including the mechanical process (Hudsen et al. 2005). The influences of mechanical processes on the overall conclusions of the study are therefore likely to be small compared to the uncertainty. It has also been suggested that the thermal impact on fractures at repository depth will be limited because the confining stress is high enough that the fractures are nearly completely closed before thermal loading (Hudsen et al. 2005).

A number of studies on saturated fractured crystalline rock repository conditions have also been conducted to investigate potential sites in Finland, Canada, and France (Chan et al. 1998; Goblet and de Marsily 2000; Löfman 2005). Regional fractures were explicitly incorporated in all three cases, and the regional scale effects of heat generation were modelled. It was shown that radiogenic decay heat can increase the rate of flow in higher-permeability fault zones by up to a factor of ten in some cases and that travel times from repository to surface could be significantly reduced, although these durations are generally considered to still be longer than the lifetime of the thermal perturbation (Chan et al. 1998; Löfman 2005). Due to the low permeability in all cases, no convection cells were predicted in these case studies. However, these investigations do not fully cover the scope of the processes required to model a site similar to Sellafield in the UK, where transmissive aquifers are present above the repository location, and groundwater at the repository depth is mobile through connected networks of fractures (Black and Brightman 1996). In the case of the Finnish model, the transfer of heat was considered to be only via conduction because of the low permeability of the host rock and overburden (Löfman 2005), and in the Canadian and French models, the influence of solute concentration on fluid density and viscosity was not considered (Chan et al. 1998; Goblet and de Marsily 2000). Therefore, this study aims to investigate the thermal disturbance on the regional groundwater flow from a hypothetical HLW repository in the Sellafield area with a coupled thermo-hydrogeological model combining mass transport and the conductive and advective transport of heat.

Geological and hydrogeological conceptual model

Sellafield is located in Cumbria, NW England, on a structural transition zone between the Lake District massif in the east, and the East Irish Sea Basin to the west (Black and Brightman 1996). The geology, as shown in Fig. 1a, b, consists of the low matrix permeability, highly fractured Ordovician Borrowdale Volcanic Group (BVG) unconformably overlain by a sequence of Carboniferous Limestone, Permian Brockram Breccia, and higher-permeability Triassic sandstones of the Sherwood Sandstone Group. The structural control on the area has produced large faults with associated damage zones of elevated permeability (Degnan et al. 2003). The fault damage zones have been shown to act as steeply dipping permeable pathways, which strongly influence groundwater flow in the cover sediments, increasing permeability by up to two orders of magnitude (Gutmanis et al. 1998). Damage zones and fault envelopes have been shown to be on the order of 50 m wide in this section. Flow, in the cover sediments especially, is focussed in this wider zone rather than the discrete fault itself (Gutmanis et al. 1998). Consequently, the faults are treated as a geomechanical facies (McDermott et al. 2006a) that incorporates the fault and associated damage zones and are discretised separately as shown in Fig. 1a, b. The conceptual hydrogeological system used in this modelling was developed from multiple borehole and geophysical studies during NIREX investigations, shown in Fig. 1b. Three separate water regimes of varying densities were identified: the Coastal Plains regime (CP), Hills and Basement regime (HB), and East Irish Sea Basin regime (EISB) (Black and Brightman 1996; Bath et al. 2006). The CP regime is topographically driven freshwater with a relatively high flow rate within the transmissive Triassic sandstone cover rocks (NIREX 1995; Black and Brightman 1996; NIREX 1997a). This regime is decoupled from the underlying, slower flux HB regime, which consists of a moderate concentration of solutes and is driven by topography in the east (c.10 km to the east of the cross-sectional boundary) and density differences with the EISB regime in the west (McKeown et al. 1999). The EISB regime is dense brine with high concentration of solutes derived from dissolution of the St Bees Evaporites in the basin to the west. EISB brine has a slow flux and is considered almost stationary in comparison with the CP regime (Black and Brightman 1996; Bath et al. 2006). In the vicinity of the hypothetical repository, flow of the HB regime is deflected upwards by the density contrast with the EISB regime (Fig. 1a) (Haszeldine and McKeown 1995; NIREX 1995; Black and Brightman 1996; NIREX 1997a; McKeown et al. 1999; Metcalfe et al. 2007).

Fig. 1
figure 1

2D geological sections of the Sellafield region showing the location of the hypothetical geological disposal facility (GDF) and the conceptual model (a) and boundary conditions used (b). Geological information is derived from best available data in NIREX Reports S95-012 and S97-012. This section has been discretised and forms the framework for numerical TH modelling. a The hydrogeological conceptual model shows the three separate flow regimes: the Coastal Plains regime, the Hills and Basement regime, and the East Irish Sea Basin regime. Relative flow rates are indicated by the size and direction of the arrows. b Boundary conditions and source terms used in the natural conditions baseline modelling. Boreholes shown are used in calibration. BH3 and BH9 are the model boundaries and are used to inform the boundary conditions. Basal surface is no-flow surface, and upper, west, and east surfaces permit flow

Full size image

Numerical modelling approach and site calibration

Modelling is undertaken using the versatile open source numerical code OpenGeoSys (Kolditz et al. 2012). Simulations are performed with a coupled thermo-hydrogeological model that also solves for the mass transport of conservative solutes. The transient flow of groundwater is calculated from the mass balance equation based on Darcy’s Law after de Marsily (1986):

$$S_{\text{s}} \frac{\partial p}{\partial t} - \nabla \cdot \left( {\frac{k}{\mu }\left( {\nabla p + \rho g\nabla z} \right)} \right) = Q$$

where S s is the specific storage, p is the pressure (Pa), t is time (s), k is the intrinsic permeability of the rock (m2), μ is the dynamic viscosity of the fluid (kg/m s), ρ is the fluid density (kg/m3), g is the acceleration due to gravity (m/s2), z is the depth to the datum (m), and Q is the volumetric flow rate (m3/s). The natural conditions steady-state model in this study was calculated with a specific storage equal to zero.

Heat can be transported via conduction and/or advection. While other studies have considered heat transport by conduction only, this work takes into account the possibility for advective heat transport. These processes are governed by the energy balance equation for heat transport after Kolditz (2002) and McDermott et al. (2006b):

$$c^{\text{m}} \rho^{\text{m}} \frac{\partial T}{\partial t} + c^{\text{w}} \rho^{\text{w}} v \cdot \nabla \cdot T - D_{\text{h}} \nabla^{2} T = \rho^{\text{m}} Q_{T}$$
(1)

where c is the specific heat capacity (J/K), p is the bulk density (kg/m3), T is the temperature (K), v is the advective velocity (m/s), D h is the heat diffusion dispersion coefficient (W/mK), and Q T is a heat source. The superscripts m and w refer to the medium and fluid components of the system, respectively. The heat diffusion dispersion coefficient is defined as

$$D_{{h_{x} }} = D_{\text{e}} + v_{x} \beta_{x}$$
(2)

where \(D_{{{\text{h}}_{x} }}\) is the heat diffusion dispersion coefficient in the x direction, D e is the effective heat diffusion coefficient (W/mK), v x is the advective velocity in the x direction (m/s), and β x is the heat dispersion coefficient in the x direction (J/mK). The effective heat diffusion coefficient is given by

$$D_{\text{e}} = \frac{{n_{\text{e}} \lambda_{\text{w}} + \left( {1 - n_{\text{e}} } \right)\lambda_{\gamma } }}{{n_{\text{e}} c_{\text{w}} \rho_{\text{w}} + \left( {1 - n_{\text{e}} } \right)c_{\gamma } \rho_{\gamma } }}$$
(3)

where λ is the thermal conductivity (W/mK), n e is the effective porosity, and the subscript γ refers to the properties of the rock.The heat dispersion coefficient (β) is given by

$$\beta = \alpha \frac{{c_{\text{w}} \rho_{\text{w}} }}{{c_{\text{m}} \rho_{\text{m}} }}$$
(4)

where α is the dispersivity (m). Mass transport is controlled by diffusion and/or advection and is described in this work by the mass balance equation for mass transport of conservative solutes:

$$\frac{\partial C}{\partial t} + \nabla \cdot \left( {vC - D\nabla C} \right) = N$$
(5)

where C is the concentration of solutes (mg/l), t is the time (s), D is the hydrodynamic dispersion coefficient (m2/s), v is the advective velocity (m/s), and N is a source term (mg/l). However, in this model N is zero, i.e. there is no modelled source or sink of solute. The hydrodynamic dispersion coefficient is defined as

$$D = D_{\text{m}} + \alpha \left| v \right|$$
(6)

where D m is the molecular diffusion, α is the dispersivity (m), and v is the advective velocity (m/s). The aim of this study is not to assess the complex processes of radionuclide fate and transport but to investigate the impacts of the heat generation on the groundwater velocities in the vicinity of the repository and the effects of heterogeneities within the system. Conservative transport is assumed.

Coupling arises from the dependency of the advective velocity calculated from Eq. 1 on the temperature and solute concentration-dependent fluid parameters of density and viscosity. In turn, the rate of mass transport and heat transport, and therefore distribution of solute and heat throughout the model domain, is largely influenced by the advective velocity. The mass transport and fluid flow are calculated iteratively, and heat transport is calculated sequentially with respect to these processes in a staggered scheme.

Fluid density and dynamic viscosity are a function of the temperature, pressure, and the concentration of solutes. In a system where thermal and concentration profiles vary laterally or vertically, such as near Sellafield, these factors significantly influence the groundwater flow regime. As a result, it is important to consider site-specific functions for the spatial distribution of fluid properties as determined by site investigations. Although density and viscosity are both functions of the pressure, water is often considered incompressible (McCutcheon et al. 1992), and over the temperatures and salinities observed in this section, the pressure has minimal effect. Therefore, densities and viscosities in this study were computed using NIREX-derived site-specific relationships of chloride concentration and temperatures in Eqs. (7, 8):

$$\begin{gathered} \rho \left( {C,T} \right) = (998.2063 + 1.266 \times 10^{ - 3} C - 1.19 \times 10^{ - 9} C^{2} ) + \\ \left( { - 0.22644 - 2.82 \times 10^{ - 6} C + 7.20 \times 10^{ - 22} C^{2} } \right)\left( {T - 20} \right) + \\ ( - 0.00361 + 2.68 \times 10^{ - 8} C - 7.17 \times 10^{ - 14} C^{2} )\left( {T - 20} \right)^{2} \\ \end{gathered}$$
(7)
$$\begin{gathered} \mu \left( {C,T} \right) = (1.002 \times 10^{ - 3} + 2.43 \times 10^{ - 9} C + 1.08 \times 10^{ - 14} C^{2} ) \\ \times \exp \left( { - 0.02358\left( {T - 20} \right) + 0.000107\left( {T - 20} \right)^{2} } \right) \\ \end{gathered}$$
(8)

where C = concentration of Cl− (mg/l) and T = Temperature (°C).

Including both the density and the viscosity functions to enable a coupled thermo-hydrogeological model is essential as these parameters heavily influence the flow velocity computed from the groundwater flow equation. This is especially so for the viscosity as this has been shown to be most sensitive over the temperature range modelled (McDermott et al. 2006b).

The two-dimensional model domain (Fig. 1a, b) was discretised from published data (NIREX 1995, NIREX 1997a, b; McKeown et al. 1999), and a finite element mesh was produced using an automatic mesh generator, Gmsh (Geuzaine and Remacle 2009). The modelled hypothetical repository is situated within an area of the BVG free from large fault zones in order to simulate a containment scenario in accordance with (McKeown et al. 1999). It is understood that the modelled repository may be far smaller than the necessary footprint of a combined HLW and ILW repository, but a larger size would only exacerbate the potential impacts on flow in this location as large fault zones would necessarily be intersected by the repository. Fault names and numbers used in this study are consistent with NIREX investigations. Hydrogeological and thermal data are derived from published data and previous modelling studies as shown in Table 1 (NIREX 1997a, b; McKeown et al. 1999). Heat and mass dispersivities are considered isotropic and are of the magnitude of half the smallest element edge length so that the Courant stability criterion of mass transport is satisfied. The heat capacity of water is considered to be constant because the variation over the modelled range of temperatures and pressures is a maximum of 1.33 % and assumptions such as single permeability values for large units are likely to have a much larger effect on the results. As a base case calibration, the model developed in this study is calibrated to match present-day temperature and solute concentration data from NIREX Deep Boreholes BH2, BH4, and BH10 because fluid density is a major driving force in the hydrogeological system (Black and Brightman 1996). Solute concentration is expressed as a density at a constant temperature such that, for example, a concentration of 1,035 kg/m3 is equivalent to 35,000 mg/l salinity. Pressure boundary conditions for the east and part of the west vertical boundaries of the model were determined through a site-specific spatial distribution density function (Eq. 6) and calibrated over the depth of the model using Cl− concentration data from published field data points in BH3 and BH9 (Bath et al. 2006). Pressure boundaries permit, but do not require, flow into or out of the model domain. The East Irish Sea Basin regime was represented by a no-flow constant concentration boundary in the west base of the model domain. The east base of the domain was modelled with a constant concentration boundary condition replicating the down-hole field data and a Neumann heat flux source term consistent with previous modelling (Fig. 1b) (McKeown et al. 1999). A thermal gradient of 20 °C/km calibrated within the downhole measured range from BH2, BH3, and BH4 was imposed as initial conditions for the natural conditions model, but past climatic variations were not considered. A constant temperature boundary of 10.67 °C was imposed over the modelled ground surface between the faults as this is the current average annual temperature at Sellafield (NIREX 1997b). A source term of half the expected annual rainfall was chosen to simulate recharge, and a constant pressure boundary of atmospheric pressure was applied at the modelled ground surface.

Table 1 Summary of the thermal and hydrogeological data used to describe the hydrogeological units (McKeown et al. 1999; NIREX 1997a, b)
Full size table

As with all regional scale groundwater flow modelling, uncertainties exist associated with representing hydrogeological units with single up-scaled parameter values, e.g. permeability, in an equivalent porous media method. This is also important to note when considering solute transport in fractured media; a local effective permeability could be represented by a single fracture with large aperture through which advective transport would be faster than the same permeability represented by a system of smaller fractures (McKeown et al. 1999). Also, in this study area the 2D nature of modelling does not take into consideration the along strike orientation of maximum permeability in the fault damage zones (Gutmanis et al. 1998).

In this study, the initial conditions, i.e. the natural conditions model, were derived from a steady-state model that aimed to reproduce the observed temperature and salinity profiles. These same properties were used to investigate variants, i.e. the transient thermal pulse models.

Natural flow conditions

In order to make any meaningful predictions of the effect of heat generation from a hypothetical repository on groundwater flow, it is necessary to produce a natural flow model that matches field data for baseline purposes, which can then be perturbed as a series of variants by adding heat, or changing repository location. The baseline simulation developed is able to recreate the three flow regimes and the dominant processes controlling flow and provides a good match to solute concentration and temperature in field data (Figs. 2, 3).

Fig. 2
figure 2

Model section was populated with site-specific rock properties (Table 1). An improved ‘fit’ of baseline model results to measured natural data was obtained by calibration to Boreholes 10, 2, and 4 along this section, which extends from BH3 to BH9. Faults within the GDF host rock are given the rock property of surrounding rock, except for F1, F202, and F2 adjacent to the GDF, which have ×10 enhanced permeability, acknowledging field data. These are conservative assumptions. A very simple calibration was achieved by altering the permeability of the lithologies within the 95 % confidence calibration range in NIREX (1997a) and using previous modelling data McKeown et al. (1999). This matches measured data, but is a non-unique result

Full size image
Fig. 3
figure 3

Natural conditions model was run to determine the natural flow paths and concentration of solutes through the section before simulation of the thermal disturbance from disposal of radioactive waste. This cross section shows the good match between the conceptual model (Fig. 1) derived from field observations and data, and the simulated results

Full size image

The major driving forces of groundwater flow in the region are the topographic high to the east (modelled by the field pressure data at BH9 as boundary conditions) and the density contrast between the Hills and Basement (HB) regime and the East Irish Sea Basin (EISB) regime in the west. Rapid near-surface flow towards the west and the coast within the transmissive Calder Sandstone and St. Bees Sandstone and the influence of the fault damage zones is predicted as a consequence. Fluid density contrasts between the HB regime and EISB regime cause the upwards deflection of the HB Regime water body towards, and through, the hypothetical repository site as is consistent with the conceptual model (Black and Brightman 1996). The low permeability of the Brockram Breccia causes a sharp transition in solute concentration (Fig. 2), suggesting that the freshwater CP regime and the HB regime are decoupled. This transition across the low-permeability Brockram Breccia also causes a rise in freshwater head. However, the model predictions do not match the measured field baseline data exactly. The sharpness of the solute concentration transition is underestimated in all cases especially BH10, and the finer scale variations within BH2 are not recreated. In BH2 where there is pressure data from (McKeown et al. 1999), the model achieves a reasonable fit in the vicinity of the repository but the increase in freshwater head is not as sharp as field data and the finer scale variations at depth are not well simulated (Fig. 4). This is due to the single material group assigned to the whole of the BVG which does not allow for these finer variations in freshwater head and solute concentration to be resolved.

Fig. 4
figure 4

Comparison between freshwater head calculated by the calibrated model and the borehole data for BH2 after McKeown et al. (1999)

Full size image

In terms of the thermal profile, a good match is achieved with available data on section in BH2 and BH4. At repository depth, there is a minor discrepancy of c.1.0 °C in both boreholes and a maximum error of 1.2 °C at c.800 m depth. However, in BH2 the thermal gradient at depth is slightly underestimated. This can be explained as a result of the application of global initial thermal conditions at the start of the model. Equilibrium temperature log data from BH2, BH3, BH4, and BH5 (off section) presented in NIREX (1997b) indicate two distinct measured geothermal gradients: a lower geothermal gradient in the cover sediments of c.18–19 °C/km and a higher gradient in the fractured crystalline basement of c.24–25 °C/km. In this simple model, the whole domain is assigned a value of 20 °C/km which gives the best model fit to the data. Therefore, the geothermal gradient at depth is underestimated, and the predicted results are lower than measured. As the error is less than 1.0 °C in most parts of the profile and the thermal trends are generally well fitted, the thermal model is considered a suitable base for a scoping investigation of the large thermal effects of radioactive waste storage which could increase temperatures locally by up to 60 °C.

During calibration, it emerged that the controlling factors on the best fit model to natural conditions were the permeability of the BVG and the initial conditions of geothermal gradient. The St. Bees Sandstone permeability was also found to have a minor impact on the thermal profile and freshwater flush above the repository. In order to match the solute concentration profile to depth, a permeability value for the BVG higher than the calibrated modal values in NIREX studies was used. However, this permeability value was consistent with other modelling work that was able to better recreate the borehole pressure distributions (McKeown et al. 1999). Increasing the permeability of the BVG still further would enable a better prediction of the sharp solute boundary change, but this would negatively influence the temperature profile. This balance between the BVG permeability and temperature profile and the finer details of the solute profile and freshwater head may lend support to further discretisation of the BVG sub-units, e.g. in NIREX calculations (NIREX 1995, 1997a). However, while it is clear that more detailed investigations are undoubtedly possible, this 2D simplified model is adequate to illustrate the principles of the thermal effects of the heat pulse on the regional system.

A transient simulation using the calibrated natural flow conditions model for initial conditions was run for a further 3000 years in order to have a baseline comparison with the thermal pulse model. The results indicate a slight rise in temperature within the repository area of less than 0.80 °C, suggesting that there is minor heat transport from depth caused by the upward flow of water through the repository. The solute concentration distribution was stable and matches well with the conceptual model and field data. The freshwater head distribution and flow paths are shown in Fig. 5.

Fig. 5
figure 5

Natural conditions model was also run for 3000 years into the future to determine the natural flow paths and freshwater head through the section before simulation of the thermal disturbance from disposal of radioactive waste

Full size image

Thermal pulse simulations

A series of simulations were made, commencing with present-day conditions, and moving forward in time initially for 3000 years and up to 6000 years into the future. Time steps were 2 years. Heat generation, decaying through geological time, is modelled as a Neumann source term of the power output data from the NAGRA spent fuel canister specifications based on British Nuclear Fuels high-level waste canisters (Fig. 6) (McGinnes 2002). The repository material properties were not modelled explicitly, and no engineered barrier system effects on thermo-hydrological transport were considered. The thermal power source term was distributed over the repository area such that the maximum design temperature of 100 °C for the bentonite buffer in the engineered barrier is not exceeded, as is consistent with saturated fractured crystalline rock repository concepts (Anttila et al. 2008). This restriction on the maximum bentonite temperature is currently being investigated for various design concepts, e.g. the temperature buffer test (TBT) in fractured crystalline rock where the heater surface reaches c.130–150 °C and temperatures at the buffer–rock interface exceed 60 °C in less than 4 years (Åkesson et al. 2012), and the HE-E heater test at Mont Terri where canister–bentonite temperatures exceeded 130 °C (Gaus et al. 2014). We therefore choose a representative value for a thermal source with the aim to show any thermal effects of heat generation rather than to comment on specific repository designs. No further reduction in the source term is considered to account for the use of a 2D section as opposed to a 3D simulation. This is because the simulated section represents the centre of the hypothetical repository, which could extend for many hundreds of metres to the north and south of this particular section. Consequently, the transport of heat in these directions will be markedly reduced by the heat generated by the rest of the repository in comparison with the transport of heat laterally and vertically; the 2D section can be thought of as being insulated in the north–south direction. Therefore, it is considered an acceptable simplification, for an initial investigation into the possible thermal disturbances at this location, to conduct a 2D simulation.

Fig. 6
figure 6

Decay heat power generation based on the BNFL waste canister data in NAGRA waste specifications (McGinnes 2002)

Full size image

The thermal pulse simulation predicts that the peak temperature in the repository of 95 °C is reached after 80 years and a large-scale thermal pulse develops over a period of a few hundred years as heat is transported away from the repository (Fig. 7). Disturbance of natural flow commences immediately and progressively increases flow rates and temperatures, along pathways similar to the natural case, through timescales of hundreds of years. A 1.0 °C temperature rise in groundwater at the surface expressions of F1 and F202 caused by heat transport from the repository is predicted after c.1100 and c.1200 years, respectively. The upward direction of flow through the repository area towards permeable fault zones of F1 and F202 promotes this advective transport of heat towards the surface. A maximum temperature rise of 13 °C is predicted at the modelled base of F1 after c.1725 years and 1.7 °C at the surface location of F1 after 2100 years. Closer to the repository in F202, a larger maximum rise of 22 °C is predicted after c.882 years and 1.20 °C at the surface after c.1720 years. Figure 7 shows the evolution of the thermal pulse over time in terms of the change in temperature from simulated natural conditions. These diagrams show that the heterogeneity in permeability caused by the fault damage zones exerts a controlling influence over the flow field and heat transport from the repository. In the initial stages, the advective heat transport is focussed on F202, but over time, the regional flow pattern leads to a larger temperature change in F1 at the surface. The freshwater flush to the west, superimposed on the advective heat ascent within the Calder Sandstone, ‘washes’ the thermal pulse to the west transferring the heat transport from F202 to F1. When the heat generation from the repository subsides, the temperature and concentration seen in the fault zones also subside back towards natural conditions.

Fig. 7
figure 7

Cross sections at 500-year intervals of modelled temperature change from natural flow conditions during a thermal pulse. Heat generated by the decay of HLW/SF in the repository causes a thermal pulse to develop. Heat transport is focussed on fault zones of elevated permeability close to the repository. (MaOD Metres above Ordnance Datum)

Full size image

The hydraulic conductivity of a hydrogeological unit is dependent on the intrinsic permeability, the fluid density, and the dynamic viscosity. Consequently, increasing the temperature of the fluid increases the ease with which fluid can flow through the unit (de Marsily 1986). As a result, the velocity of warm water flow increases within the fault zones and flux rates increase (Fig. 8). F202 close to repository depth is affected most by the thermal pulse and also shows the largest increase in flow rate by a factor of over 1.5: a similar magnitude of change as modelled at the Finland crystalline repository site (Löfman 2005), but smaller than at the Whiteshell Research Area, Canada, where fault permeabilities varied by up to six orders of magnitude over the country rock (Chan et al. 1998). The timing indicates that the increased rate of flow is first focussed on F202 closest to the fault where peak flow rate occurs at c.940 years. Reduction in fluid density and viscosity in F202 means that the effective hydraulic conductivity is increased and flow is focussed towards this fault zone in the early stages of heating. The peak flow rate is reached in the F1 cover St. Bees Sandstone nearly 1000 years before the F1 basement BVG, suggesting that the interaction between individual heterogeneities and these structures with the regional flow field is just as important as the presence of the heterogeneities themselves.

Fig. 8
figure 8

Flow rate evolution with time at locations within faults F1 and F202 in m3/s)

Full size image

One factor that could have a potential impact on the effective hydraulic conductivity, which is not considered in this current work, is the mechanical response to the thermal pulse. The thermal expansion of the rock could potentially lead to a reduction in permeability in the fault damage zones by reducing the fracture apertures. However, it is thought that fracture apertures would already be at a residual minimum at repository depth due to the confining pressure and therefore less effected by the thermal pulse (Hudsen et al. 2005). Such effects will be offset by the thermal expansion of the fluid, which increases volume with warmer fluid. Previous work focussing on the thermally induced mechanical effects in a far-field study also predicts a reduction in permeability in the vicinity of a hypothetical repository (Min et al. 2005). However, a TM-coupled model is used that does not consider the effects of the heat on fluid flow, and therefore the change in effective hydraulic conductivity, or indeed the increase in permeability caused by repository excavation (Min et al. 2005). Detailed modelling of these mechanical effects on the thermo-hydrogeological system is not covered here.

Large-scale convection cells are not predicted by this model, as is consistent with previous regional scale site-specific modelling, despite elevated permeabilities in comparison with Whiteshell (Canada) and Olkiluoto (Finland) (Chan et al. 1998; Löfman 2005). Large differences in permeability between the fault zones and the country rock could theoretically lead to convection developing within the fault zones, e.g. López and Smith (1995), but the ranges of permeability in this current work suggest this is unlikely. López and Smith (1995) provide a ‘permeability space’ that predicts the range of host rock vs fault permeabilities at which conduction, advection, and stable and unstable convections are dominant thermal transport processes. A 3D representation would be required to predict convective transport in the fault zones for this site, but the results of advective transport in the fault zones of this simple scoping study coincide well with the expected zones of advection in (López and Smith 1995) and confirm the necessity of including advective heat transport in the model.

The reduction in density and increased velocity of upwards flow causes more saline water to be drawn up from depth (Fig. 9). There is a predicted increase of 0.97 kg/m3 in the solute concentration at the base of F1 but just 0.16 kg/m3 at the surface. However, while this is consistent with previous modelling (Löfman 2005), the potentially significant effects of repository excavation on the salinity profile (Löfman 2005) have not been considered. High solute concentration in the vicinity of the repository could potentially cause adverse effects on bentonite buffer swelling capacity (Hudsen et al. 2005). The reduction in fluid density causes a reduction in freshwater head, which in turn leads to a minor influx of fluid with lower concentration of solutes from the east. This, coupled to the regional flow field, prevents the solute concentration in the repository from increasing by more than 1.5 kg/m3, suggesting that, for this location, heat generation is unlikely to increase solute concentration enough to cause adverse effects. However, it is possible that this effect may be overestimated in this modelling because an increase in pressure in the vicinity of the repository is possible due to the thermal expansion of water in a confined aquifer. This would reduce the effect on freshwater head and the influx of fluid with lower concentration of solutes. Further research developing the model to include a coupled mechanical process and investigate the extent to which the permeable fault damage zones would allow this excess pressure to dissipate would be required to quantify this effect.

Fig. 9
figure 9

Cross sections at 500-year intervals of modelled concentration change from natural flow conditions during a thermal pulse. Heat generated by the decay of HLW/SF in the repository causes a reduction in fluid density and a subsequent rise in more saline water from depth. The reduction in fluid density also causes a reduction in freshwater head and an influx of less saline water from the complex fault system of F2. The regional westerly flow regime flushes the higher-concentration fluid to the west so solute concentration in the repository does not increase more than 1.5 kg/m3. (MaOD Metres above Ordnance Datum)

Full size image

Radionuclide fate and transport modelling is not undertaken, and therefore, the complex processes of migration and lifetime of radionuclides in the geosphere, essential for a full site performance assessment, are not addressed here. However, in fractured crystalline rock disposal concepts, such as Sellafield, groundwater transport is the dominant mechanism for radionuclide migration, and the increase in rate of flow within F1 and F202 results in a substantial reduction in travel time to surface for fluid passing through the geological barrier. As chloride is considered a conservative chemical species in transport modelling, advective transport of chloride was used in this study as a proxy for conservative contaminant transport. There are no modelled sources or sinks of Cl−. This scoping model predicts that a thermal pulse induced by emplaced heat-generating waste promotes significantly enhanced transport towards the surface in the early stages of a repository lifetime. Failure of the multi-barrier safeguards during this stage could significantly influence contaminant travel time. For example, in the thermal pulse scenario, after less than 1,260 years, particles which track water released from the repository western edge reach the St. Bees Sandstone aquifer via faults F1 and F202. This represents over a third of a reduction in travel time compared to natural conditions (over 1960 years) and, as a result of the sustained elevated temperatures from the thermal pulse, a reduction in travel time to the surface can be expected. As is evident from Fig. 7, the thermal pulse at depth persists over a significant period of time although the intense peak falls away relatively quickly. In order to investigate the influence of the slow dissipation of the thermal pulse, particles were released into the model at intervals of 0, 1000, 2000, 3000, and 4000 years. Particles released at the point of containment and after 1000 years have very similar travel times from the repository edge to the St. Bees Sandstone aquifer of c.1300 years. After this initial thermal peak, the travel time for particles released after 2000 and 3000 years increases by c.100 years for every 1000 years. However, particles released after 3000 years and 4000 years show a smaller increase in travel time of 40 years because the temperatures in the thermal pulse are closer to that of the surrounding conditions. The travel time for these particles is still over 350 years shorter than natural conditions.

Conclusions

Using the excellent dataset at Sellafield, UK, as an example, a coupled thermo-hydrogeological model has been developed, which is calibrated to reproduce the chloride concentrations and temperature profiles observed in field data. Post-closure thermal disturbances are then investigated through a heat source representing a hypothetical repository.

A thermal pulse around the repository develops within decades, which increases the effective hydraulic conductivity and velocity of groundwater, thus increasing the transport rate of potential leachate from the repository. Flow pathways are controlled by the interaction between the natural flow field and the steep heterogeneities of the fault damage zones.

A rise in groundwater temperature of over 1.0 °C is predicted at the surface after just 1100 years. Advective water flow along fault damage zones, at 1.5× natural velocity, transports heat and solute towards the surface. If failure of the engineered multi-barrier system were to occur at the point of confinement/closure, transport of conservative chemical species representing repository leachate from the western edge of the repository could contaminate an aquifer within 1300 years. The long-lived thermal pulse causes travel times of conservative chemical species, released more than 4000 years after closure, to remain reduced by over 350 years compared to natural conditions.

Further modelling would be required to fully quantify the post-closure thermal effects of high-level waste disposal at this site by extending this section into 3D and including the effects of both repository excavation and engineered barrier emplacement in a coupled thermo-hydro-mechanical-chemical model. However, this preliminary study illustrates the requirement to include advective transport at this site. Alteration of the natural flow regime due to heat from HLW and SF and flow pathway interactions with heterogeneities in the geological barrier (such as permeable fault damage zones) are important considerations in the development of a post-closure safety case.