Abstract
The Prototype Repository (PR) tunnel is located at the Äspö Hard Rock Laboratory near Oskarshamn in the southeast of Sweden. In the PR tunnel, six full-sized deposition holes (8.37 m deep and 1.75 m in diameter) have been constructed. Each deposition hole is designed to mimic the Swedish reference system for the disposal of nuclear fuel, KBS-3V. The PR experiment is designed to provide a full-scale simulation of the emplacement of heat-generating waste. There are three phases to the experiment: (1) the open tunnel phase following construction, where both the tunnel and deposition holes are open to atmospheric conditions; (2) the emplacement of canisters (containing heaters), backfill and seal in the first section of the tunnel; and (3) the emplacement of canisters, backfill and seal in the second section of the tunnel. This work describes the numerical modelling, performed as part of the engineered barrier systems (EBS) Task Force, to understand the thermo-hydraulic (TH) evolution of the PR experiment and to provide a better understanding of the interaction between the fractured rock and bentonite surrounding the canister at the scale of a single deposition tunnel. A coupled integrated TH model for predicting the wetting and the temperature of bentonite emplaced in fractured rock was developed, accounting for the heterogeneity of the fractured rock. In this model, geometrical uncertainties of fracture locations are modelled by using several stochastic realizations of the fracture network. The modelling methodology utilized information available at early stages of site characterization and included site statistics for fracture occurrence and properties, as well as proposed installation properties of the bentonite. The adopted approach provides an evaluation of the predictive capability of models, it gives an insight of the uncertainties to data and demonstrates that a simplified equivalent homogeneous description of the fractured host rock is insufficient to represent the bentonite resaturation.
A geological disposal facility (GDF) is considered by waste management organizations as a solution for long-term management of higher activity radioactive waste. A GDF is a multi-barrier system, isolating the radioactive waste from the biosphere. It consists of a natural barrier (i.e. host rock) and engineered barriers. Engineered barriers include the wasteform, the waste disposal container and the buffer material to protect the container.
An example of an engineered barrier system (EBS) is the Swedish (SKB 2006, 2010) KBS-3V concept, which has been developed for the final disposal of higher activity radioactive waste. It comprises a copper canister containing four spent fuel assemblies emplaced in deposition holes, within disposal tunnels, and surrounded by a protective buffer material. The buffer material consists of highly compacted sodium bentonite, which is an absorbent aluminium phyllosilicate clay rich in montmorillonite. It is the montmorillonite that allows bentonite to swell in the presence of water (Börgesson 1985). Its low permeability and its capacity to swell makes bentonite an important material for use in the geological disposal of radioactive waste. When saturated, it protects the disposal containers from detrimental thermal, hydraulic, mechanical and chemical processes and retards the release of radionuclides from any failed containers. Waste management organizations have identified bentonite as an ideal buffer material and in-situ experimental programmes have therefore been implemented to enhance understanding of the interaction of the emplaced bentonite with a host rock. There are several examples of such in-situ experiments, these include: The Full-Scale Engineered Barriers Experiment FEBEX (ENRESA 2000); The Engineered Barrier Emplacement Experiment EB (Mayor et al. 2005); The Bentonite-Rock Interaction Experiment BRIE (Vidstrand et al. 2017); and the Prototype Repository (PR) experiment, data from which are used in the modelling presented in this work.
To increase understanding of such barrier systems, a numerical model that describes the wetting behaviour of bentonite in response to interactions with a host rock and disposal container has been developed. The model was informed by and validated against data, collected from the PR experiment. Its purpose is to provide a predictive and informative tool, able to make predictions of saturation and temperature for this long-term experiment and to give insights into the impact introduced on these predictions through uncertainties in the fracture network.
PR experiment
The Äspö EBS Task Force, which is organized and led by the Swedish Nuclear Fuel and Waste Management Company (SKB), provides a forum for international specialists and modelling groups to collaborate in furthering the understanding of processes relevant to bentonite within the context of final disposal of radioactive waste. As part of this Task Force, modelling the PR experiment (Goudarzi & Johannesson 2009) is undertaken at the Äspö Hard Rock Laboratory (HRL) near Oskarshamn (see Fig. 1) in the southeast of Sweden. The Äspö HRL has a 3.6 km long tunnel, which was excavated to a depth of approximately 460 m in the early 1990s.
The location of the PR, TASG and the BRIE tunnels at the Äspö HRL. This figure is reproduced from SKB (2016).
The PR tunnel consists of two sections, described as Section 1 and Section 2, with four full-sized deposition holes (8.37 m deep and 1.75 m in diameter) located in Section 1 and two in Section 2 (see Fig. 2). Each deposition hole contains a copper canister equipped with electrical heaters to represent the thermal loading of high-heat generating wastes. The two Sections were backfilled with a mixture of bentonite (30% by weight) and crushed rock (70% by weight), with concrete plugs installed at the end of both sections. This design is based on the Swedish reference system for the disposal of nuclear fuel (KBS-3V). During the latter part of 2001, installation of the canisters, buffer, backfill and seal took place in Section 1, whereas Section 2 was installed in the early part of 2003.
Sections of the PR tunnel. Sections 1 and 2 are backfilled with a mixture of bentonite and crushed rock.
The TASG tunnel running sub-parallel to the PR tunnel was bored to allow easy monitoring and control of certain aspects of the PR, as well as to provide power to the heaters installed within each canister. The PR was heavily instrumented with sensors (a large number of which failed after the tunnel drainage was closed) in both the tunnels and deposition holes (363 in Section 1 and 394 in Section 2). The sensors recorded temperature, groundwater pressure and mechanical loading.
Inflow rates were also sampled at several locations of the PR tunnel and in the six deposition holes. A characterization investigation of the fractured rock at the Äspö site local to the PR tunnel (Stigsson et al. 2001) was also conducted through an exploratory drilling programme, involving flow logging, borehole pressure build-up tests and interference tests. The main goals of the PR modelling related to the development of a methodology that:
furthered understanding of bentonite resaturation in the presence of a fractured host rock;
modelled the influence of heaters installed in deposition holes to simulate the emplacement of spent fuel;
established the effect of interactions between deposition holes and the influence of tunnel backfilling on bentonite wetting.
Modelling methodology
The adopted methodology is comprised of the following steps:
Modelling inflow. This modelling activity aims to describe the heterogeneous flow of groundwater into the PR experiment from the rock. This is achieved using discrete fracture network models (DFN), which are based on a statistical description of fractures in the rock. The DFN concept is a powerful approach, able to model stochastically and deterministically (where sufficient characterization has taken place) the individual flow conduits in fractured rock. Based on this concept, several purely hydraulic DFN calculations are performed for a variety of stochastic realizations of the fracture network. The aim is to simulate the inflow rate for each deposition hole in the two sections of the PR tunnel.
Upscaling. In this activity the heterogeneity of the host rock, as represented explicitly in the DFN concept, is captured in a continuum representation of the rock. This is achieved through an upscaling process in which the properties of the DFN are described by equivalent continuous porous medium (ECPM) block properties on a suitable scale (i.e. the hydraulic conductivity of the two concepts are equivalent on this scale). The purpose of this is to allow additional capability by modelling not only inflow rates from the heterogeneous rock, but also bentonite resaturation and temperature at sensor locations, as required in post-installation phases of the PR experiment.
Modelling the PR experiment. The PR experiment is divided into three phases. The first one corresponds to pre-installation conditions (i.e. open tunnel and deposition holes), whereas the last two phases correspond to post-installation conditions (i.e. introduction of tunnel backfill and emplacement of bentonite buffer). These phases are described below.
Phase 1: Excavation of tunnel and deposition holes.
Phase 2: Emplacement of bentonite buffer and canisters in Section 1 of the tunnel, followed by backfilling of Section 1. A plug was positioned between Sections 1 and 2, and the power to canisters 1–4 was turned on.
Phase 3: Section 2 of the tunnel was backfilled and plugged, after the emplacement of canisters 5 and 6 and associated bentonite buffers. Power to canisters 5 and 6 was switched on.
To model the PR experiment simulation tools were selected that had both DFN and CPM capabilities. These tools had to be numerically efficient to enable the main features of the deposition tunnel to be captured without making the calculations numerically intractable. The ConnectFlow (Joyce et al. 2010; Amec Foster Wheeler 2015) modelling platform is an efficient numerical tool suitable for DFN and CPM modelling, with parallelized solvers, providing thermo-hydraulic-chemical capabilities (THC). It has an interface to FLAC3D (ITASCA Consulting Group), which allows for mechanical (M) coupling. Therefore, the modelling of the PR experiment was carried out using DFN and CPM (TH) capabilities of ConnectFlow.
Geometric model description
The 3D-model domain used an explicit representation of the PR tunnel (approximated by a cuboidal tunnel of equal cross-sectional area), the six deposition holes, the TASG instrumentation tunnel and the main access tunnel. Each deposition hole was discretized with a fine enough grid to accurately capture the groundwater in the fractured rock and the wetting behaviour of the bentonite. The mesh (see Fig. 3) generated consists of c. 270 000 elements with the majority of the refinement concentrated close to the PR tunnel volume.
(a) The 3D grid; (b) the explicit representation of the Prototype Repository, TASG and Main Access tunnels; (c) slices through the computational grid for xz, yz and xy planes, and a plan view of the refinement of a deposition hole.
DFN modelling of the PR-tunnel deposition hole inflow
The PR is located in granitic bedrock, with groundwater flow occurring through a series of interconnected fractures. A DFN modelling approach is adopted for the bedrock that explicitly represents the flow conducting fractures, providing representation of the flow and transport at a local level.
A characterization of the fractured rock at the Äspö site local to the PR Tunnel, conducted by (Stigsson et al. 2001), provides fracture statistics consistent with the bulk properties of the rock mass. Based on exploratory drilling campaigns of the PR, a DFN parametrization was derived and is utilized here; a detailed presentation is given in (Tsitsopoulos et al. 2018). Fracture statistics, given as distributions for intensity, size, orientation and transmissivity, have been inferred from fracture traces on tunnel walls and intersections with boreholes. The stochastic network was generated using the collated statistics. Realizations of these stochastic networks are able to capture the variability of the flow in the rock. In addition, deterministic features have been characterized through field measurement at global (repository) and local (deposition hole) scale. In the analysis that follows two DFN models are generated.
Global refinement DFN model: Stochastic fracture statistics aimed at capturing the variability of the flow in the rock surrounding the PR. In addition, two large-scale deterministic structures that intersected the PR tunnel were included.
Local refinement DFN model: Uses the same stochastic fracture network statistics from the global refinement DFN with the addition of minor deterministic fractures conditioned to accurately represent the spatial distribution of inflows within the instrumented deposition holes (Tsitsopoulos et al. 2018).
The deterministic features explicitly represented in the DFN models are shown in Figure 4. Deposition holes 5 and 6 can be seen to be intersected by deterministic features in the local refinement DFN model.
Positioning of the deterministic fractures as implemented in the two DFN scenarios.
For both DFN models, 40 stochastic realizations are generated using ConnectFlow and DFN simulations were performed under open tunnel and deposition holes’ conditions (i.e. atmospheric pressure and with no backfill emplaced). These simulations provide a prediction of the total inflow to the sections of the PR tunnel and each of the deposition holes. Selected results from both models are compared to measured values in Figure 5.
Calculated inflows in deposition hole 6 for 40 realizations, generated using the locally- and globally-refined DFN models, are compared with the measured inflows (central red dashed line).
Evidently, there is variability among the simulated inflows, obtained from different stochastic realizations. This behaviour is exaggerated for realizations corresponding to the globally-refined model which has not constrained the occurrence of higher flowing fractures. Overall, introduction of deterministic features in the locally-refined model results, as expected, in better constrained realizations.
In Figure 6, simulated inflows obtained for realization 14 of the locally-refined DFN model are compared to measured inflows of the deposition holes and sections of the PR tunnel. A more detailed presentation of this information is given in the work of (Tsitsopoulos et al. 2018). There is a tendency to overestimate the inflow except in the case of deposition hole 6. The reason for this is the locally-refined DFN model has used fracture characterization information to locate the positions of fractures in deposition hole 6. Such features tend to dominate the inflow and from Figure 6 we can infer that the inclusion of such features is essential in obtaining a more reliable description of the inflow.
Computed inflows for realization 14 of the locally-refined DFN model compared with the measured inflows at deposition holes (DH) and tunnel sections.
Upscaling the DFN model to an ECPM model
DFN modelling allows for a more natural description of the heterogeneity of the host rock. However, a continuum model is more appropriate if other processes, such as the heat conduction or the wetting of emplaced bentonite are required. It is necessary that the continuum model represents the inherent heterogeneity of the fracture network and so a methodology referred to as upscaling, implemented in ConnectFlow, has been developed. It calculates both a directional permeability tensor and fracture kinematic porosity for grid cells, based on the fracture network (Jackson et al. 2000). The upscaling methodology implements a flux-based approach, requiring several flow simulations to be performed through the DFN model on the scale of individual grid cells. For grid cells that do not correspond to fractures in the DFN, permeabilities and porosities are assigned corresponding to a granitic rock matrix. The resulting continuum model is termed the ECPM.
An illustration of the heterogeneity introduced in the ECPM model is shown in Figure 7. Even though the upscaling methodology is independent of element size, it can be observed that as the element size reduces, the upscaled permeabilities have a greater contrast due to the effect of reducing the volume that permeability is averaged over.
Illustration of the upscaling process. On the top figure stochastic and deterministic fractures are coloured by transmisssvity (DFN representation). The figure at the bottom shows the permeability tensor in the z (vertical) direction kzz, derived through the process of upscaling.
As the upscaling is a computationally demanding process, not all realizations of the DFN models are upscaled to an ECPM. Five realizations of the globallly-refined model and four of the locally-refined model are upscaled. Realizations demonstrating a favourable comparison with measured inflows are selected, in particular with deposition hole 6, as this is where most of the high-quality sensor data is located.
Modelling of phases 1–3
Having obtained several heterogeneous ECPM models, which correspond to the selected DFN realizations, focus is directed towards testing the model by simulating temperature and saturation at several locations in the PR experiment. The adopted approach is designed to model the experimental phases described earlier and repeated here for convenience:
Phase 1 is a purely hydraulic calculation. It considers open tunnels where the pressure is set to be atmospheric. A steady state calculation is performed to set up an initial pressure field.
Phase 2 is a coupled TH calculation. It considers the time between installation of Section 1 (deposition holes 1–4) and the installation of Section 2 (deposition holes 5–6).
Phase 3 is a coupled TH calculation using a heterogeneous representation of the rock. It considers the time post-installation of Section 2. The repository tunnel for Section 2 is backfilled, the canister and bentonite buffer installed to each deposition hole and the plug between Section 2 and the rest of the PR tunnel is added. A numerical simulation is performed for 2500 days.
Initial and boundary conditions (ICs and BCs) for the modelling phases are as follows. For Phase 1 the hydrostatic pressure is provided as an IC. The steady state pressure field obtained in Phase 1 is then passed as an IC in Phase 2 and the pressure field obtained at the end of Phase 2 is passed as an IC for the Phase 3 simulation. In terms of temperature, Phase 1 is isothermal, therefore an IC of ambient temperature (15°C) is provided at the initiation of Phase 2. The temperature field obtained at the end of Phase 2 is provided as an IC for Phase 3. The top and bottom boundaries of the domain are sufficiently far from the tunnel excavations that a Dirichlet boundary condition of hydrostatic pressure and ambient temperature can be assumed. Neumann pressure and temperature BCs are applied to the lateral extent of the domain with a ‘no flow’ condition prescribed (i.e. a symmetry condition). This assumption is considered appropriate due to the number of similar additional tunnels of the HRL that reside outside the modelling domain. Finally, atmospheric pressure and ambient temperature BCs are set on the TASG tunnel and the main access tunnel connected to the PR Tunnel. The BCs and ICs applied in the three modelling phases are illustrated in Figure 8.
Schematic diagram of the PR tunnel, illustrating the initial and boundary conditions used in the three-phase modelling process.
Heat is introduced to the PR experiment through heaters installed in the emplaced canisters. The power decay of the canisters has been approximated as a piecewise constant function (Tsitsopoulos et al. 2018). The power input methodology adopted in Phases 2 and 3 as well as the power input are presented in Figure 9. The canisters were affected by two significant instances of power outages, the details of which are shown in Figure 9.
Power evolution used in the modelling for deposition holes 1 through 4 for Phase 2 (top) and all six deposition holes for Phase 3 (bottom) accompanied by a schematic illustration. A common power input is applied to deposition holes 1 through 4 (blue line), whilst unique power profiles are assigned to deposition hole 5 (green) and 6 (red). The first power outage effected all canisters and was caused by damage due to excessive water pressure in the buffer and backfill. In addition, problems with canister 6 resulted in a second power outage lasting around 2 months.
Material specifications for the host rock, bentonite and tunnel backfill are presented in some detail in Tsitsopoulos et al. (2018). For Phases 1, 2 and 3 the H and TH simulations are carried out in ConnectFlow, where it assumed that there is direct contact between the rock wall and bentonite blocks, i.e. no air-filled gaps. The governing equations, numerically solved, consist of the conservation of mass:
In the equation above, the residual pressure, PR, is related to the total pressure, PT and given by:
Here g [m s−2] is the gravitational acceleration, ρ0 [kg m−3] is a reference density, k is the permeability tensor, kr is the relative permeability and μ [kg ms−1] is the dynamic viscosity.
The heat transport equation as implemented in ConnectFlow is given as:
In the above equation, the term H [W m−3] corresponds to a heat source/sink, while information regarding the average heat capacity (ρc)a and the average thermal conductivity Γa, is presented in Table 1, accompanied by the functional forms of density and viscosity.
Functional forms of liquid density, viscosity, average thermal conductivity, average heat capacity and liquid suction
The default van Genuchten retention curve has been adopted to model the rock and backfill of the tunnel (Van Genuchten 1980), whereas for the bentonite a modified van Genuchten model is introduced accompanied by a cubic law for the relative permeability (Åkesson et al. 2010). The formulae used for saturation as a function of suction and relative permeability as a function of saturation are summarized in Table 2, and a detailed analysis can be found in (Tsitsopoulos et al. 2018).
Retention models used in the modelling
Results obtained from the transient TH simulation of Phase 3 of the PR experiment are presented in Figures 10 and 11. A maximum temperature of just below 60°C is observed at the outer boundaries of the buffer boxes of deposition holes 1–4, at this point temperature in deposition holes 5 and 6 is ambient. The simulation suggests that it takes 100 days for the maximum temperature on the external boundaries of deposition holes 5 and 6 to increase by approximately 30°C. The temperature for deposition holes 5 and 6 increases until the second year of the simulation, where the peak calculated temperature is approximately 55°C. The computed temperature evolution for deposition holes 5 and 6 is at steady state between years two and four, whereas the variations observed after the fourth year of the simulation are related to the power input of Canister 6. In contrast, the temperatures simulated in deposition holes 1–4 appear to slowly, but monotonically, increase with time, reaching a maximum of 60°C at the end of the simulation. These higher temperatures simulated in Section 1 of the PR tunnel are a consequence of the close proximity of the four canister heat sources, whereas Section 2 of the tunnel only contains two deposition holes, with the power output from canister 6 significantly reduced due to difficulties with intermittent power failure of the heater.
Temperature distribution of the PR through the modelling of Phase 3 of the PR experiment. Figures (a–h) show the temperature evolution at different times, as computed at the external boundaries. Here [s] indicates seconds, [d] days and [a] years.
The first nine surface plots present the computed pressure evolution in the bentonite on a slice taken in the x–y plane. The final slice shows a representation of the hydraulic conductivity of the heterogeneous rock (outside the bentonite ring) and homogeneous bentonite. This illustrates why the resaturation of the bentonite ring is strongly influenced by the fracture network.
Figure 11 shows the pressure evolution along a horizontal slice through the canister and buffer location. The suction within the bentonite-filled ring surrounding a canister is decreasing with time, but not uniformly. The lack of homogeneity in the resaturation process can be explained by the presence of a fracture network in the immediate vicinity of the buffer. At the bottom of Figure 11, the horizontal slice shown focuses on materials that are distinguished by their hydraulic conductivity. The fracture zones are coloured green, whereas elements belonging to the buffer ring or the rock are dark blue. By direct comparison, it can be observed that the resaturation is directly driven by fractures and the assumptions made in the model about how groundwater moves through them. On this basis, the performance of a heterogeneous numerical model in representing the resaturation at a sensor level is strongly dependent on the accurate representation of the fracture network in its proximity.
Model evaluation
As most of the available sensor data covering the full duration of Phase 3 of the PR experiment was collected in deposition hole 6, the numerical model is evaluated against these measurements. The positioning of the ten sensors recording pressure and temperature is shown in Figure 12. Numerical simulations of Phase 3 were performed for five realizations of the globally-refined DFN model and for four realizations of the locally-refined DFN model.
Positioning of the 10 sensors located within the bentonite buffer of deposition hole 6. Each sensor measures temperature and liquid suction over time.
Computed temperature profiles are found to be in a good agreement with the measured data for all sensors (Tsitsopoulos et al. 2018). Deviations between measurements and simulated temperature evolutions, presented in Figure 13, are less than a few degrees at all sensors considered. Discrepancies can be attributed to features not incorporated into the model, such as the heterogeneity of thermal properties of the fractured rock and changes in thermal conductivity of bentonite as it resaturates. The larger deviations between measurements and simulated temperatures observed at earlier times in the experiment are likely to be a consequence of the timescales required for the bentonite to swell sufficiently to fill the outer gap (which is not modelled) between the bentonite and deposition hole wall. At early times it is expected that the gap acts as a thermally insulating layer. At later times, larger deviations are attributed to power outages that have not been explicitly modelled.
Computed temperature profiles obtained for five different realizations are directly compared with the experimental data gathered at four sensors. The variation in the computed temperature profiles among realizations is shown to be negligible.
As shown Figure 13, the temperature evolutions in the bentonite buffer is independent of the discrete fracture network (i.e. realization number or DFN model). This is mainly a direct consequence of the dominance of conduction over advection as the heat transfer mechanism. The model assumes the thermal conductivity of the bentonite is not dependent on saturation. It is recognized that this is a simplification; however, it is expected that, because the buffer is relatively thin and the rock is the major thermal conductor, this would have a relatively minor impact on temperature (a few degrees).
The sensors provide data for liquid suction, which can be used to inform rates of bentonite wetting as resaturation of bentonite is a suction driven process. A comparison of the modelled and observed liquid suction values is given in Figure 14. Here, several realizations of the two DFN models are considered in order to determine the impact different realizations have on the simulations. For sensors WB657 and WB658 it can be observed that the overall variability between realizations has been significantly reduced for the locally-refined DFN model.
Computed suction profiles obtained at three sensors for five realizations of the globally-refined DFN and four realizations of the locally-refined DFN model. Equivalent experimental data recorded at the saturation sensors are also shown.
This trend is also applicable for sensor WB663, but the reduction is now less prominent. Overall, the deterministic specification of fractures that intersect deposition holes 5 and 6, with properties based on mappings and hydraulic tests, has brought the range in predictions closer to the observed data. Further, it has reduced the variability between realizations. The variability still present is a consequence of the stochastic far-field fracture network and its ability to supply groundwater to the conditioned deterministic fractures intersecting the deposition hole.
Resaturation of the bentonite buffer in the heterogeneous PR model is dominated by the location of the fractures intersecting the deposition hole, which provide a pathway for groundwater ingress. Hence, an accurate reproduction of the observed data would require an accurate representation of the DFN, at least near the deposition holes. Gathering such data would be very difficult, particularly since a GDF would typically involve thousands of depositions holes. It is therefore worth emphasizing the versatility of the DFN approach as:
one of the few modelling approaches that can provide insight into the variability of flow at a resolution level of an individual fracture;
a viable method for providing an understanding of the resaturation rates that are likely to be observed across a GDF, by utilizing site characterization statistics.
Sensitivity analysis
This section focuses on the sensitivity of the predictions to key thermal and hydraulic parameters. The analysis consists of five numerical experiments, investigating the following sensitivity cases:
the host rock intrinsic permeability;
the bentonite buffer intrinsic permeability;
the tunnel backfill intrinsic permeability;
the thermal conductivity of the host rock;
power input into the system of canisters.
To illustrate the sensitivity of those results on material properties, a single realization (realization 14) of the locally-refined DFN model was selected to act as the base-case. This realization was identified as optimal from the resaturation predictions as presented in Tsitsopoulos et al. (2018). Due to the dominance of conduction over advection and the assumption of constant thermal properties, thermal effects are found to be unaffected by changes in hydraulic parameters. Therefore, for the first three sensitivity studies, attention is directed to the hydraulic effects, whereas for the last two cases only thermal effects are presented. The variants from the base-case scenario for each sensitivity study are summarized in Table 3. Simulation results are compared to measured data for sensors WB658 and WB611 in Figures 15 and 16. In addition, simulated saturation (for the first three sensitivity cases) and temperature (for the last two sensitivity cases) of the bentonite in deposition hole 6 at the end of Phase 3 (2500 days) are presented in Figures 17 and 18, respectively.
Suction profiles obtained for the base-case scenario (realization 14 of the locally-refined DFN), and the sensitivity cases considering the rock, bentonite and tunnel backfill intrinsic permeability. Results are directly compared with data measured at the sensors.
Temperature profiles obtained for the base-case scenario (realization 14 of the locally-refined DFN), and the sensitivity cases considering the rock conductivity and power input. Results are directly compared with data measured at the sensors.
Saturation (S) as computed at the end of the numerical simulation of Phase 3 on the bentonite in deposition hole 6, for the three cases investigated as part of the rock matrix, bentonite and tunnel backfill permeability sensitivity studies. The points used in the 3D scatter plot are the nodes of the computational grid responsible for the representation of the bentonite.
Temperature as computed at the end of the numerical simulation of Phase 3 on the bentonite in deposition hole 6 for the three cases investigated as part of the thermal conductivity and power input sensitivity studies. The points used in the 3D scatter plot are the nodes of the computational grid responsible for the representation of the bentonite.
Sensitivity studies considered and their variants
In the sensitivity study that considers host rock-matrix intrinsic permeability, the base-case value was set to kR = 10−21 [m2]. Now two additional scenarios are considered i.e. kR = 10−20 and kR = 10−22. As anticipated, an increase/decrease of the rock matrix permeability results in higher/lower rates of saturation for sensor WB658 (top of Fig. 17), and therefore full saturation is achieved at earlier/later time. This is not the case for sensor WB611, where a change in the rock permeability leaves the suction evolution unaffected. This behaviour results from the presence of a fracture in the immediate vicinity of the sensor. Since fractures drive the resaturation process, the influence of the rock-matrix permeability is less apparent in this area. Therefore, we can conclude that for the range of variants considered, the intact rock permeability has little influence at sensor level in the presence of a fracture zone. It is worth noting that varying the permeability of the rock matrix by an order of magnitude around the base-case value does not have a linear impact on the computed saturation rates. For sensors not in the vicinity of a fracture (i.e. WB658) it can clearly be observed that the impact on increasing the intact rock permeability on wetting of the bentonite becomes more significant. This is because the rock permeability can become comparable with the bentonite permeability, and hence is able to supply sufficient groundwater for the wetting process.
Figure 15 shows the sensitivity case where the permeability of the bentonite is varied. In this figure the evolution of suction for the base-case (sensors WB658 and WB611) of kR = 10−20 are compared with the variant cases of kB = 6.410−21 and kB = 2.210−21. The selected values capture the range that have been experimentally determined (see Olsson et al. 2013) for the buffer in deposition hole 6. Here, suction profiles computed at sensors WB658 and WB611, are dependent on the bentonite permeability. This is a contrast to the rock-permeability sensitivity study, where computed suction profiles for sensor WB611 (positioned close to a fracture) were found to be independent on the rock permeability. Even though the resaturation process is driven by the neighbouring fractures and is quicker, the measured hydraulic impact is weighted by the ratio of the hydraulic permeability of the bentonite and that of the fracture zone.
The final hydraulic related sensitivity study considers the effect of the tunnel backfill permeability. In the base-case model, the backfill is a mixture of bentonite and crushed rock with a permeability of kTB = 310−18. In the first of two variants, a backfill of crushed rock is assumed, which is a higher conductive material with permeability kTB = 10−16; for the second variant the same hydraulic properties of the bentonite buffer used in the base-case model (i.e. kTB = 10−20) are assigned to the entire backfill material. Suction evolutions computed at the sensors for both variants, as shown at the bottom of Figure 15, do not show any variation with backfill permeability. It can be seen in Figure 17, for the low permeability case, that there is an effect on saturation at the top of the bentonite ring that sensors, due to their positioning (located further down the hole), do not capture. This serves as a reminder that the sampled data, although providing detailed information regarding the temperature and pressure evolutions, only cover discrete locations in the bentonite and cannot necessarily be used to imply bentonite behaviour away from the sensors.
A further sensitivity study considers the effect of changing the thermal conductivity from the base-case value of λR = 2.63 [W mK−1], to λR = 2.3 and λR = 3.6. As expected, and shown in Figure 16 (top), reduction of the thermal conductivity by 12.5% resulted in an increase of the simulated temperature (of the order of 5%), whilst an increase of the thermal conductivity by 36.9% had the opposite effect and resulted in a drop of the temperature profiles on the order of 14%. Further, it can be observed that for the variant where the thermal conductivity is increased, the temperature profiles are slightly shifted to earlier times, suggesting that the response of the system to the power input is more rapid.
The temperature distribution in the bentonite as obtained at the end of simulation is shown in Figure 18 with the maximum temperature varying from 56°C (base-case) to 48°C when λR = 3.6 and 60°C when λR = 2.3.
The final sensitivity analysis considers changes in power input to the system. Two scenarios are considered, the first assumes a 10% increase in power input for canister 6, whereas the second adopts a 10% power increase for all canisters. The system responds linearly, as observed in Figure 16, so that an increase of 10% in the input energy results in an approximate 10% increase of the computed temperature ΔT = T − Tamb. It is worth noting that the temperature profiles observed for the two variants display minor differences of the order of 1–2°C, in particular for sensors in the vicinity of deposition hole 5. Profiles with the most elevated temperature correspond to the variant with increased power input for all canisters. This is to be expected, as more energy is introduced into the system. This sensitivity study demonstrates the localized nature of the thermal effects. There is only a marginal difference in the simulated temperature of the bentonite surrounding canister 6 for the variant where the power input is increased for both canister 5 and canister 6 and for the variant where the power input is increased for canister 6 only.
Conclusions
In this work a methodology suitable for large-scale modelling of coupled TH processes was presented and applied to the PR experiment. Numerical simulations were performed using ConnectFlow's DFN and CPM capabilities.
Application of the methodology in the modelling of the PR allowed for a heterogeneous representation of the host rock surrounding the PR through DFN modelling. Two DFN models were developed, termed as global refinement and local refinement models. Both utilized fracture statistics local to the PR tunnel and large-scale structures. The latter also included a deterministic specification of fractures intersecting deposition holes 5 and 6, reflecting fracture mappings taken from the PR experiment. For both DFN models, numerous stochastic realizations were generated and used to predict inflows that had been measured along sections of the PR tunnel and its open deposition holes. A number of realizations were then selected, according to a criterion of prediction consistency, with the measured inflows. These realizations were then upscaled to an equivalent continuous porous media representation, suitable for modelling the three phases of the PR experiment. The vast majority of available sensor data in terms of temperature and pressure corresponded to Phase 3 (i.e. bentonite buffer installed in deposition hole 6) and numerical modelling focused on the ability to reproduce the sensor measurements. Interpretation of simulation results emphasized that:
Temperature evolutions computed at the sensor locations are in a good agreement with the experimental data (see Fig. 13). Deviations by a few degrees observed in few sensors can be attributed to the assumption of homogeneous thermal properties for the rock/bentonite.
Resaturation of bentonite, emplaced in low permeability fractured host rock, is dominated by water flow through the fracture of the adjacent host rock. To fully understand the spatial distribution of groundwater flow in the vicinity of the hole, a detailed characterization of the local fracture network would be required. The Locally-refined DFN model is able to incorporate this additional information gathered from the experiment, while still accounting for uncertainty caused by the less well characterized regions of the fractured rock. Additional characterization data improved the prediction of suction evolution (see Fig. 14), as compared to the globally-refined model and also constrained the variability among stochastic realizations.
The adopted modelling methodology has a broader versatility. Through the use of site characterization fracture statistics, stochastic simulations can provide an understanding of the range of resaturation rates that are likely to be observed, not just in these detailed experiments, but more broadly at a GDF which may comprise of many thousands of deposition holes.
Analyses also included a sensitivity investigation examining the uncertainty introduced in the modelling through assumptions made in the physical parameters characterizing thermal and hydraulic properties. The key findings are:
The host rock permeability can have a significant impact on the resaturation of the bentonite (see Fig. 17), however only in the case where resaturation is not dominated by adjacent fractures. For the cases considered this impact is not linear; an increase in rock permeability is found to have a more significant effect than a decrease in permeability.
Uncertainty in bentonite permeability can have a significant impact on the computed saturation profiles, with variations in bentonite permeability values causing a larger impact on bentonite wetting than variations in rock permeability (see Fig. 17).
A backfill material of higher permeability could be used without noticeably altering the resaturation process in the bentonite (see Fig. 17). In this study an increase of two orders of magnitude caused no observable change to the wetting of bentonite.
The thermal conductivity of the host rock has an important effect on the computed temperature profiles (see Fig. 18).
Acknowledgments
The authors would like to thank Radioactive Waste Management Limited for supporting our participation in the SKB Task Force on Engineered Barrier Systems. The authors would also like to thank the reviewers whose comments have been incorporated into this paper.
Funding
The authors would like to thank Radioactive Waste Management Ltd for providing financial support for this study.
- © 2018 Nuclear Decommissioning Authority. Published by The Geological Society of London. All rights reserved
