## Abstract

A key component of the site comparison planned for the deep geological disposal of spent fuel and high-level waste (SF/HLW) in Switzerland is the assessment of the evolution of repository-induced perturbations in the repository nearfield associated with thermal effects from heat production due to radioactive decay of radionuclides, as well as gas pressures developing in the backfilled underground structures from the anaerobic corrosion of the steel waste canisters and tunnel support materials. The assessment of such effects is integrated in the site comparison through safety indicators used to evaluate repository performance. In this context, probabilistic assessments need to integrate the uncertainty of the entire ensemble of input parameters, and estimate the propagation to these indicators in a reliable and computationally efficient manner. This paper presents the development of a methodology for an indicator-based assessment of heat- and gas-induced effects in a SF/HLW repository in Opalinus Clay integrating a probabilistic treatment of parametric uncertainty. The workflow is demonstrated using preliminary data, repository configurations and indicators. Complementary simulations are presented to demonstrate the feedback to the optimization of repository design in order to mitigate repository-induced effects that can potentially compromise the safety function of the engineered and natural barriers.

The safe disposal of spent fuel and high-level waste (SF/HLW) is a major challenge for the nuclear energy industry. In Switzerland, the National Co-Operative for the Disposal of Nuclear Waste (Nagra) proposed three candidate siting regions for the SF/HLW repository at the end of Stage 2 of the Swiss site-selection plan, the so-called Sectoral Plan for Deep Geological Repositories (SGT). In each siting region, the Opalinus Clay Formation was proposed as the host rock formation. As a decision basis, preliminary analyses have been elaborated, addressing the proposed repository configurations in the three siting regions: Jura Ost (JO), Nördlich Lägern (NL) and Zürich Nordost (ZNO). Among other things, the analyses addressed the thermal effects from heat production due to the radioactive decay of radionuclides in SF/HLW (Senger *et al.* 2014; Papafotiou & Senger 2016*a*, *b*), as well as gas pressures developing in the backfilled underground structures from the anaerobic corrosion of the steel waste canisters and tunnel support materials (Papafotiou & Senger 2014, 2016*a*, *b*). A key component in the site comparison planned in the SGT Stage 3 is the assessment of the evolution of these perturbations in the repository nearfield – hereafter referred to as repository-induced effects – using the site-specific information gained from field investigations. The assessment of repository-induced effects feeds into the site comparison through safety indicators used to evaluate repository performance for different repository configurations. In this context, probabilistic assessments need to integrate the uncertainty of the entire ensemble of input parameters, and estimate the propagation to these indicators in a reliable and computationally efficient manner.

This paper presents the development of a methodology for an indicator-based assessment of heat- and gas-induced effects in a SF/HLW repository in Opalinus Clay integrating the probabilistic treatment of parametric uncertainty. The workflow is demonstrated in this paper using preliminary data, repository configurations and indicators. The workflow is part of a criteria-based indicator approach that is being currently developed in support of the site-selection process in SGT Stage 3 and will integrate component-scale with repository-scale model-based assessments using site-specific data gained from the site investigations. The simulations presented in this paper therefore aim to demonstrate the probabilistic workflow and do not constitute quantitative assessments of repository performance. Moreover, complementary simulations are presented here to demonstrate the feedback to the optimization of repository design in order to mitigate repository-induced effects that can potentially compromise the safety function of the engineered and natural barriers. The probabilistic assessments were conducted based on simplified thermo-hydraulic (TH) problem statements to reduce the computational effort associated with the solution of the fully coupled thermo-hydraulic-mechanical (THM) problem. The simulations that address these issues comprise:

Sensitivity analyses using the Morris Screening methodology (Morris 1991) to identify the most sensitive parameters with regard to pressure and temperature in the repository and surrounding host rock.

Monte Carlo analyses to estimate the envelope of simulated temperatures and pressures, as well as the risk of uncertain parameters resulting in failure of performance indicators, defined as temperatures and pressures exceeding certain criteria.

Repository configurations with regard to different repository depths and distances between the SF/HLW emplacement tunnels.

Repository configurations with regard to different emplacement densities of SF and HLW in the emplacement tunnels.

The simulations for (a), (b) and (c) are developed based on a 2D model (Senger *et al.* 2014) that describes the generic configuration of a single SF/HLW emplacement tunnel, taking into consideration the combined effect of heat emission and gas accumulation (component scale). To account for the range of parametric uncertainties and repository configurations, the model input parameters and geometries are varied accordingly. The simulations for (d) are developed based on a 3D model (Papafotiou & Senger 2016*a*) that represents a generic configuration of a SF/HLW repository, taking into consideration the combined effect of heat emission and gas accumulation in the different components of the repository and during the entire repository cycle including the operation and post-closure phases (repository scale). To account for different emplacement densities, the same loading of waste (i.e. the total number of SF and HLW canisters) is emplaced in different portions of the emplacement tunnel, while the remaining tunnel contains no waste. The entire modelling framework has been developed with the equation-of-state module EOS5 for water and hydrogen of TOUGH2 (Pruess *et al.* 1999) and the inverse version iTOUGH (Finsterle 2007).

## Methodology

### Basic design of the SF/HLW repository

The generic concept for the HLW repository included the following main components: access ramp and shafts, operations and construction tunnels, the pilot repository, emplacement tunnels for long-lived intermediate-level waste (ILW), the underground rock laboratory (URL), seal sections, and emplacement tunnels for the SF/HLW (Fig. 1). The generic repository layout was for a waste inventory estimated based on 50 years of operation of the existing nuclear power plants. The reference distance between emplacement tunnels was 40 m.

The waste inventory included steel canisters with spent fuel elements from the operation of nuclear power plants and vitrified HLW from the reprocessing of spent fuel. Reference dimensions of steel canisters assumed a canister length of 4.6 m and a radius of 0.525 m. Several possible emplacement variants were considered for the SF/HLW repository. In this report the design variant that was used corresponded to the reference concept for a deep SF/HLW repository, as described in Nagra (2010). Waste containers were emplaced on a pedestal of compacted bentonite blocks, whereas the rest of the tunnel was backfilled with compacted bentonite pellets (Fig. 2). The tunnel walls were reinforced with shotcrete lining, steel mesh and rock anchors. These ensured that the tunnels sustained the required geometry during the emplacement of waste canisters and backfilling with bentonite.

Heat generation from the radioactive decay of the waste were estimated for different waste types based on Nagra's inventory developed for the preliminary assessments of SGT Stage 2 (Nagra 2013). The present analysis adopted the reference heat generation rates for a MOX/UO_{2} canister, which had the highest output in the inventory (Fig. 2). Gas generation from the corrosion of the waste canisters can also have a significant effect on the long-term evolution of the SF/HLW nearfield (Diomidis 2014). In the framework of SGT Stage 2, comprehensive sensitivity studies were performed to estimate the gas generation rates in a SF/HLW repository. A range of gas generation rates was provided as input for the sensitivity analyses of gas release (Papafotiou & Senger 2014), representing a reference value, as well as the upper and lower bounding values. In this paper, the gas source term drew on the reference gas generation rates from these analyses. The rate of gas generation from canister corrosion corresponded to 0.01 m^{3} of hydrogen gas (SATP)/m/year for 60 000 years. The rate of gas generation from the corrosion of tunnel installations corresponded to 0.06 m^{3} of hydrogen gas (SATP)/m/year for 2000 years (Fig. 2).

Material properties were adopted from the preliminary database developed in SGT Stage 2. The TH properties of the Opalinus Clay were required for the simulation of heat transport through the host rock and its impact on pore pressure. Thermal properties of the Opalinus Clay were taken from earlier studies (Johnson *et al.* 2002) and TH calibrations of temperature profiles from deep boreholes in the candidate siting areas (Papafotiou & Senger 2013*a*, *b*, *c*). Gas-related properties of the Opalinus Clay had been derived in Senger *et al.* (2013). Two-phase flow and associated parameterizations were described based on the assumption that the Opalinus Clay can be described as an effective homogeneous medium.

### Development of the SF/HLW emplacement tunnel model

On the component scale, the simulations were developed using the 2D model from Senger *et al.* (2014) that described the generic configuration of a SF/HLW emplacement tunnel, taking into consideration the combined effect of heat emission and gas accumulation. Senger *et al.* (2014) examined deterministically the early thermal period associated with the re-saturation of the bentonite and heating of the emplaced SF/HLW canisters, taking into account potential regional geothermal conditions associated with basal heat flow, as well as the long-term effects associated with gas generation from corrosion. The 2D model was implemented using integrated finite differences (IFD) with nested gridding near the emplacement tunnel, allowing a detailed representation of the SF/HLW canister and bentonite buffer materials. The lateral extent of the model was determined by the symmetry axis through the SF/HLW tunnel and the half-distance between two adjacent SF/HLW tunnels. The model grid was extended by 1D elements connected to the land surface and to the Crystalline basement at the bottom, representing the local hydrogeological and geothermal conditions in the candidate siting regions ZNO, JO and NL. The TH properties of the hydrogeological units above and below the host rock formation were taken from Papafotiou & Senger (2013*a*, *b*, *c*).

To account for different distances between emplacement tunnels, the mesh and lateral extent of the 2D model of a SF/HLW emplacement tunnel (Senger *et al.* 2014) were modified based on the principle of superposition. Six alternative models were implemented to account for a distance between the emplacement tunnels equal to 65, 75, 90, 110, 125 and 135% of the base configuration. It is therein assumed that the ensemble of emplacement tunnels could be represented with a single half-tunnel due to symmetry. The implication from this approach is that the thermal effect was overestimated as the half-tunnel model does not consider lateral dissipation of heat through the rock surrounding the emplacement field. Moreover, the 2D representation does not consider gas transport along the excavation damaged zone (EDZ) and repository tunnels, which is a key mechanism for gas release in the repository (Papafotiou & Senger 2016*a*, 2017). In the current context, the 2D model was used for the development and demonstration of the probabilistic workflow, and not for the quantitative assessment of repository performance. Gas generation was included in the workflow development and demonstration but gas-induced effects need to be further evaluated with the integration of analyses on the repository scale to account for the effects of effective gas-release pathways.

The evaluation of the TH evolution in the emplacement tunnel and nearfield was carried out using four different indicators:

temperature evolution at the canister surface (indicator

*T1*);pressure evolution at the bentonite–canister interface (indicator

*P1*);pressure in the Opalinus Clay at a distance of 20 m above the tunnel (indicator

*P2*);pressure in the Opalinus Clay at the lateral boundary, the respective midpoint between adjacent emplacement tunnels (indicator

*P3*).

These were based on preliminary safety indicators used to evaluate the barrier function of the bentonite and surrounding rock in the performance assessment of the geological repository. Temperature and pressure inside the tunnel are related to chemical and mechanical alterations of the bentonite buffer; and pressure in the Opalinus Clay nearfield is related to the mechanical integrity of the host rock. The performance of each investigated repository configuration was evaluated based on failure criteria defined for each of the preliminary indicators. In this demonstration, it was assumed that the failure criterion for temperature (*T1*) is peak temperature exceeding 150°C, whereas the failure criterion for pressure (*P1*, *P2* and *P3*) is peak pressure exceeding the lithostatic value at the corresponding depth below ground level. The SF/HLW emplacement tunnel model geometry representation with the performance indicator locations are shown in Figure 3. The indicators selected here aimed to demonstrate the methodology. The indicators used for model-based assessments of repository-induced effects would need to be revised depending on the aims of the analysis (i.e. integrating repository-scale analyses for gas-induced effects), and may also depend on the repository configurations and site-specific conditions and input parameters determined with the site investigations.

### Development of the SF/HLW repository model

On the repository scale, a 3D model was developed with a nested IFD mesh which resulted in more than 300 000 model cells representing 25 emplacement tunnels for SF (including a pilot repository), three emplacement tunnels for HLW, each with a length of 700 m (Fig. 4), and the remaining repository components in the Opalinus Clay host rock and confining units. Based on the preliminary emplacement schedule (Papafotiou & Senger 2016*a*), the three outermost emplacement tunnels in the main repository were the HLW emplacement tunnels, whereas the remaining 24 tunnels were considered as SF emplacement tunnels. The total storage volume and length of the three emplacement tunnels in the pilot repository corresponded to one SF emplacement tunnel. Heat generation rates were prescribed in the model accordingly per emplacement tunnel type, and took into account the distribution of compartments and intermediate seals in each case. The IFD grid along the emplacement tunnels was discretized with 10 m cells. The heat generation rates corresponded to those used in the 2D model of a SF/HLW emplacement tunnel (Papafotiou & Senger 2016*a*).

To assess alternative repository designs, the source terms used for heat and gas were adjusted to account for emplacement of the same volume of waste in either 100 or 75% of the total tunnel length in each emplacement tunnel. The configuration with only 75% of the tunnel length used for waste emplacement corresponded to approximately two compartments less in each of the SF (two compartments comprised 25% of the heat generation cells) and HLW tunnels (two compartments comprised 20% of heat generation cells), and one compartment less in each of the pilot emplacement tunnels.

The evaluation of the TH evolution in the emplacement tunnel and nearfield was carried out with the following metrics:

Time histories of temperature and pressure (which corresponded to gas pressure when gas was present and liquid pressure when not) at selected locations in the repository, as shown in Figure 5.

Time histories of pressure (which corresponded to gas pressure when gas was present and liquid pressure when not) at selected locations in the Opalinus Clay, as shown in Figure 5.

### Methodology for probabilistic analysis

The probabilistic analyses were performed using the 2D model of a SF/HLW emplacement tunnel. The aim of the probabilistic analysis was to identify the most sensitive parameter datasets associated with the highest risk of failure of selected performance indicators. These datasets were subsequently used to optimize repository configurations with regard to performance indicators, wherein alternative configurations were sought so that the risk of failure was reduced or even eliminated for the assumed parameter distributions.

Assessment of the relative importance of different input parameters in models and understanding their respective impact on model predictions is applicable to different fields of science and engineering. Determining the impact of parameter variations is known as sensitivity analysis and can broadly be divided into local and global analyses (Sobol 1990; Saltelli *et al.* 2008; Sudret 2008). In most cases, a global sensitivity analysis is performed to analyse a model in general or to investigate, quantify and rank the effects of parameter variation or parameter uncertainty on the overall model uncertainty (Anderson & Burt 1985; Winter *et al.* 2006; Oladyshkin *et al.* 2012). Global sensitivity analysis also helps engineers to produce more robust designs, and decision-makers to optimize uncertainty in problems related to risk assessment (de Barros & Rubin 2008; de Barros *et al.* 2009).

The first component of the probabilistic framework developed in this study was a computationally efficient sensitivity analysis technique, called Morris Screening, to identify the most sensitive input parameters, as well as input–output non-linearities, in the model (Morris 1991). The basic principles of the Morris Screening methodology and the implementation in the simulation framework using the SF/HLW emplacement tunnel model are described below. The most sensitive screened parameters were subsequently combined with the estimated uncertainty in the heat source in a Monte Carlo framework to assess envelopes of pressures and temperatures, and the risk associated with failure of the selected performance indicators at three repository depths at the respective candidate sites.

#### Morris Screening

The Morris Screening method (Morris 1991) is a computationally efficient heuristic method for performing ‘global sensitivity’ analyses without resorting to exhaustive sampling across the parameter space. The Morris method may be summarized as below:

Run the model with an initial set of (best guess) inputs.

Randomly perturb an input variable by a small factor (Δ), while keeping others constant.

Calculate the ‘elementary effect’ for the perturbation, by taking the absolute difference in the two outputs |

*F*(*x*+ Δ) −*F*(*x*)| divided by Δ.Repeat steps (2) and (3) for all

*n*variables, selecting each at random while keeping other variables constant (at their previous values).Repeat steps (1), (2), (3) and (4)

*R*number of times, with different starting points. Step (1) is repeated with a randomly selected initial set. This leads to a total of*R*(*n*+ 1) model evaluations. Keep track of all elementary effects for a given variable perturbation.Calculate the mean (

*μ*) and standard deviation (*σ*) for the elementary effects for each variable.Plot the mean and standard deviations on a scatter plot.

The mean and standard deviation plot from the Morris Screening method reveals sensitive and insensitive parameters in addition to important relationships between the input and output variables. Variables that have low *μ* and *σ* are classified as relatively unimportant. Variables that have high *μ* but low *σ* are important but have linear input–output dependencies, whereas variables with high *μ* and *σ* are important with non-linear impacts on the output.

The selection of parameters to be subjected to Morris Screening was based on previous deterministic sensitivity studies that indicated the range of expected temperatures and pressures for the different repository configurations (Senger *et al.* 2013, 2014). In that context, deterministic scenarios were simulated using alternative values for the TH properties of the emplaced materials and Opalinus Clay. Additional insights on the sensitivity of the thermal response to input parameters were gained by the prediction–evaluation analyses performed for the full-scale (FE) experiment in Mont Terri (Papafotiou *et al.* 2015). Temperatures at the heater surface in the initial heating phase (i.e. the first 2 months of heating) indicated sensitivity to the canister and bentonite thermal properties. Temperatures in the tunnel and nearfield during the early heating phase (i.e. the first 3 years of heating) showed pronounced sensitivity to the thermal conductivity of the bentonite, and the thermal conductivity and heat capacity of the surrounding rock. Preliminary predictions of pressure indicated sensitivity to the thermal properties, as well as the permeability and pore compressibility of the Opalinus Clay.

The initial set of uncertain parameters comprised a large set of hydraulic and thermal properties assigned to the tunnel materials, the EDZ and geological units (Opalinus Clay, and over- and underlying formations) that were explicitly described in the model. Several of these parameters could be screened out based on the current understanding of the TH evolution in the emplacement tunnel and the nearfield. For example, the earlier analyses showed that thermal effects are, in general, expected to dissipate within the upper and lower confining units (UCU and LCU), so that parametric uncertainty in the respective geological units above and below these could be excluded from the analysis. Furthermore, two-phase conditions existed in the bentonite following emplacement but developed in the surrounding rock only at later times when the thermal pulse had propagated through the surrounding rock. Gas saturations in the rock, in general, were low, so they did not impact the thermal properties. It was thus assumed that the uncertainty in thermal conductivity was sufficiently described by the range of values used for saturated conditions. Similarly, two-phase properties (i.e. van Genuchten parameters) were not included in the perturbed parameters as it had been shown that the impact on thermal effects was relatively insignificant compared to that of absolute permeability and the thermal properties of the tunnel materials and surrounding rock. Although gas generation was included in the current workflow development, the uncertainty in the gas-related parameters (including gas generation and release pathways through the repository) is not addressed in this study and will be evaluated in the next stage of the phenomenological analyses for gas-induced effects. The set of parameters considered for Morris Screening thus included:

thermal conductivity (

*λ*), heat capacity (HC), absolute permeability (*K*), pore compressibility (*C*_{p}), heat expansion coefficient (*C*_{t}) and porosity (*ϕ*) of the Opalinus Clay;*λ*, HC,*K*,*C*_{p}and*C*_{t}of the EDZ;*λ*, HC,*K*,*C*_{p}and*C*_{t}of the bentonite pellets;*λ*, HC,*K*,*C*_{p}and*C*_{t}of the bentonite blocks;*λ*and HC of the waste canister;*λ*, HC,*K*and*C*_{p}of the UCU;*λ*, HC,*K*and*C*_{p}of the LCU.

The absolute permeability of the Opalinus Clay was assumed anisotropic, with *K _{xy}* equal to 5 times

*K*according to the considerations from Senger

_{z}*et al.*(2013). The absolute permeability of all other formations and materials was assumed to be isotropic. The permeability of the EDZ was perturbed as an independent parameter; however, using a mean that was always larger than that of the undisturbed Opalinus Clay. Thermal conductivity in this demonstration was assumed to be isotropic. The uncertainty of the anisotropy in thermal and hydraulic properties will be addressed in subsequent assessments incorporating site-specific data.

Based on these considerations, a set of 31 parameters was used for the Morris Screening, as shown in Table 1. Normal and log-normal distributions were assumed for the uncertain parameters. With the initial parameter run and the 31 independent parameter perturbations, each Morris Screening set consisted of 32 simulations. A total of 10 Morris sets were conducted to calculate the mean and standard deviation of the system response for each independent parameter perturbation. The number of Morris sets was determined based on the convergence of the mean and standard deviation, and therefore could vary depending on the uncertain parameters and the purpose of the application. The elementary effect of interest is the response of the peak value of the four performance indicators. The sensitivity of each perturbed parameter was therefore evaluated based on the elementary effect on peak temperature at the canister surface, peak pressure at the canister–bentonite interface, peak pressure at 20 m distance above the tunnel and peak pressure between the tunnels.

#### Monte Carlo analyses

For the SF/HLW emplacement tunnel model, the Monte Carlo analysis was implemented using the most sensitive parameters selected from the Morris Screening. Monte Carlo uncertainty analysis consists of an exhaustive sampling across the parameter space to characterize the sensitivity and uncertainty in the model output. While extremely powerful, Monte Carlo analysis is computationally expensive because it typically requires several hundred or thousands model calculations. Given the high computational expense of Monte Carlo simulations, more efficient approaches, generally referred to as stratified sampling, have been developed that divide the sampled population into smaller groups (strata). In stratified random sampling, the strata are formed based on members’ shared attributes or characteristics. A random sample from each stratum is taken in a number proportional to the stratum's size when compared to the population. These subsets of the strata are then pooled to form a random sample. The most commonly used form of stratified sampling is the Latin Hypercube sampling procedure, where the range of each parameter is divided into several intervals of equal probability and a value is selected at random from each interval (Helton 1993; Helton & Davis 2003).

A total of eight material parameters were selected as uncertain parameters for the Monte Carlo analysis, namely: *λ*, HC, *K* and *C*_{p} of Opalinus Clay; *K* of the EDZ; HC and *K* of the bentonite pellets; and *K* of the UCU (see the results from the Morris Screening). The parameters were sampled from their probability density functions using the mean and standard deviations given in Table 1. Parameters were sampled randomly, as well as with Latin Hypercube sampling.

The repository depths and associated geothermal and hydrogeological conditions for ZNO, JO and NL were used for the Monte Carlo analysis. The reference hydraulic properties for deep Opalinus Clay (Senger *et al.* 2013) were used as the best estimate values (mean of probability density function) in simulation variants with repository depths of 600 and 750 m bgl (below ground level), whereas the reference hydraulic properties for shallow Opalinus Clay (Senger *et al.* 2013) were used for the best estimate values at a depth of 450 m bgl. Thermal conductivities of the Opalinus Clay and of the over- and underlying formations corresponding to the respective candidate sites and repository depths were taken from Papafotiou & Senger (2013*a*, *b*, *c*) and Senger *et al.* (2014). The standard deviations used for sampling the uncertain parameters were selected to capture the estimated range of values and were assumed to be the same for all repository depths.

The uncertainty in the heat source term was also included in the analyses by sampling the heat source rate between the high and low bands of the average heat source rate (Nagra 2013). For this, an additional uncertain parameter is introduced that corresponds to a scaling factor of the heat source rate.

Based on these considerations, Monte Carlo simulations were performed in six different variants (Table 2):

MC1: an initial set of 200 simulations randomly sampling all of the 31 uncertain parameters (repository depth of 600 m bgl);

MC2: a set of 1000 simulations randomly sampling the eight most uncertain parameters selected with Morris Screening (repository depth of 600 m bgl);

MC3: a set of 1000 simulations randomly sampling the eight most uncertain parameters selected with Morris Screening, and the scaling factor for the heat source (repository depth of 600 m bgl);

MC4: a set of 1000 simulations using Latin Hypercube sampling to sample the eight most uncertain parameters selected with Morris Screening, and the scaling factor for the heat source (repository depth of 600 m bgl); an additional set of 1000 simulations was run using a distance of 50 m between the tunnels (MC4b);

MC5: a set of 1000 simulations using Latin Hypercube sampling to sample the eight most uncertain parameters selected with Morris Screening, and the scaling factor for the heat source, assuming a repository depth of 750 m bgl;

MC6: a set of 1000 simulations using Latin Hypercube sampling to sample the eight most uncertain parameters selected with Morris Screening, and the scaling factor for the heat source, assuming a repository depth of 450 m bgl; an additional set of 1000 simulations was run using a distance of 50 m between the tunnels (MC6b).

The Monte Carlo simulations were evaluated using the four performance indicators *T1*, *P1*, *P2* and *P3* that were also used for evaluating parameter sensitivities with Morris Screening. For each of these performance indicators, failure criteria were defined as peak temperature and/or pressure exceeding certain values, given in Table 3.

#### First-order reliability method (FORM)

To provide approximations of the risk of parameter combinations resulting in failure, the Monte Carlo samples resulting in failure of any of the four criteria were screened out and further evaluated using a level II first-order reliability method (FORM) (Ang & Tang 1984; Madsen *et al.* 1986).

For each Monte Carlo simulation that resulted in failure, the parameter vector *X* comprising the eight parameters (or nine parameters when uncertainty in the heat source is included) was transformed to a set of independent standardized normal variables *U* (Hohenbichler & Rackwitz 1981; Ditlevsen 1981) by computing the difference between *X* and the best estimates normalized by the standard deviation. The next step was to find the set of parameters closest to their mean values (with highest probability) that would lead to the system response being at the prescribed target (e.g. failure criteria for *T1*–*P3*). This ‘limit-state’ response can be given by *g*(*U*) = 0, which is the system state based on the transformed parameter set, corresponding to the prescribed target. This can be described as a minimization problem, wherein the reliability index *β* is estimated as the minimum distance from the origin to the failure surface in the *U*-space. Once *β* has been found, the exceedance probability (risk) can be given as the standard normal cummulative probability function for the critical system state *β*.

### Repository layout optimization

Alternative configurations were also used to assess the potential for optimizing the repository layout with input from the probabilistic analyses. The alternative configurations comprised variants of distance between emplacement tunnels (2D model) and variants of emplacement density within tunnels (3D model).

In the 2D model of a SF/HLW tunnel, the alternative configurations of distance between the emplacement tunnels corresponded to 65, 75, 90, 110, 125 and 135% of the base-case scenario (40 m distance). For each repository depth, the following variants were simulated using the alternative distances between the tunnels:

the base case corresponding to the best estimates of input parameters;

the input parameters resulting in the maximum probability of failure in

*T1*,*P1*,*P2*and*P3*(input from the Monte Carlo simulations);the input parameters resulting in the minimum peak pressure in

*P1*–*P3*(input from the Monte Carlo simulations).

Two configurations of waste emplacement were implemented in the 3D model of the SF/HLW repository corresponding to waste emplacement in the entire length of the SF/HLW emplacement tunnels and emplacement in 75% of the length of the tunnels. The following simulation cases were performed using both configurations for waste emplacement:

the base case corresponding to the best estimates of input parameters (repository depth of 700 m bgl);

cases with the input parameters resulting in the maximum probability of failure in

*P1*,*P2*and*P3*in the 2D emplacement tunnel model.

The parameter sets resulting in failure were extracted from the Monte Carlo variants with Latin Hypercube sampling. For this, the Monte Carlo variant MC5 was used, which corresponded to the deepest of the repository configurations (750 m bgl) and in accordance with the TH parameter set used in the base case of the 3D SF/HLW repository model (Papafotiou & Senger 2016*a*).

## Results

### Sensitivity of parameters

For the first Morris Screening set, parameter values in the first simulation run are those used in the base-case simulation (TH00) of Senger *et al.* (2014). One independent parameter is then perturbed by a small factor for the remaining 31 simulations in the set. The perturbation factor is randomly selected from a normal distribution between −2 and 2%. The order of parameters perturbed in the set is randomly selected and all previous parameter changes in the set are maintained. In the second to the 10th Morris Screening set runs, the initial parameter values (run 1) for the perturbed parameters are selected randomly from the mean, corresponding to the base case parameters, and the standard deviation for each variable shown in Table 1. The computation time on a single machine (Intel^{®} Xeon^{®} CPU E5-1620 v3 @ 3.50 GHz) for each set varied between 5 and 12 h depending on the input parameters.

The time histories of temperature and pressure at the waste canister surface (performance indicators *T1* and *P1*) and in the Opalinus Clay (performance indicators *P2* and *P3*) simulated in Morris sets 1–10 are shown in Figure 6. While some parameters result in changes in temperature or pressure that are not graphically visible (i.e. pressure curve overlaps with the previous one), other parameters result in more distinct changes. The synthesis of these elementary effects will result in the parameter sensitivity with regard to peak values of the selected performance indicators *T1*, *P1*, *P2* and *P3*.

The overall TH evolution is similar to that observed in the deterministic sensitivity cases discussed in Senger *et al.* (2014). Temperatures at the waste canister reach their maximum values between 3 and 10 years (Fig. 6). Peak temperatures vary between approximately 108°C (Set 10) and 127°C (Set 9) across the 10 Morris sets. After approximately 800 years, temperature at the canister has decreased below 100°C in all simulations. It is indicated that temperature may increase temporarily despite the continuous decrease in heat output at the canister (Fig. 2). This is due to the TH evolution associated with the 2D model considerations (i.e. no-flow lateral boundaries, buffer re-saturation, EDZ representation). Peak pressure at the waste canister is reached between 200 years (Set 9) and 2000 years (Set 5), and ranges between 10.5 MPa in Set 10 and 14.2 MPa in sets 2 and 3. Pressure at a distance of 20 m above the tunnel peaks at between 200 years (Set 9) and 700 years (Set 5), with maximum values of between 10 MPa in Set 10 and 14.2 MPa in Set 5. Two (or more) peaks of pressure between the tunnels are observed in most cases. Maximum pressure between the tunnels is reached between 200 years (sets 5 and 10) and 800 years (Set 2), whereas the lowest and highest peak pressure values correspond to 10.2 MPa in Set 10 and 14 MPa in Set 3, respectively. Figure 7 shows the time histories of gas saturation in the bentonite next to the waste canister (location of *P1*), indicating the evaporation of porewater near the canister during the early heating phase followed by re-saturation of the buffer. Spatial plots of pressure, temperature and gas saturation after 60 and 60 000 years are given in Figure 8. Gas saturation in the locations of indicators *P2* and *P3* remain zero during the early heating phase when overpressures are induced by the differential thermal expansion between porewater and the rock framework. In the late post-closure phase (i.e. 60 000 years), gas saturation in the host rock reaches values of the order of 1 × 10^{−3}.

Each parameter perturbation thus results in an elementary effect as shown in Figure 6, which is expressed in this demonstration as the absolute change in peak temperature or pressure observed in *T1*, *P1*, *P2* and *P3*, respectively. This result is illustrated in Figure 9, which shows the absolute change in peak temperature and pressure for each perturbed parameter and set of Morris screenings.

An equal weighting of the response in all four performance indicators resulted in the Morris plot shown in Figure 10, indicating that the most sensitive parameter is the absolute permeability of the Opalinus Clay followed by the thermal conductivity, the heat capacity and the pore compressibility of the Opalinus Clay, the absolute permeability of the EDZ, the thermal conductivity and absolute permeability of the bentonite pellets, and the absolute permeability of the UCU. Alternatively, different weights can be used for performance indicators *T1*–*P3* to increase sensitivity with regard to temperature or pressure at specific locations. Figure 10 shows the parameter sensitivity resulting when only indicators *P1*–*P3* are used in the analysis. It is indicated that sensitivity to parameters of the emplacement materials and the EDZ is reduced, whereas sensitivity to parameters of the UCU and LCU is increased.

### Monte Carlo simulations

The Monte Carlo simulations are performed in eight sets to account for different variants related to sampling the uncertain parameters and the repository configuration (Table 2). The computation time for 1000 runs on a single machine (Intel^{®} Xeon^{®} CPU E5-1620 v3 @ 3.50 GHz) varied between 2 and 4 days, depending on the input parameters. The time histories of temperature and pressure simulated with variant MC4 are shown in Figure 11 (only 200 simulations are shown in the figure). The mean of peak temperature and pressure is indicated in each case, together with the bandwidth of 1 SD. An overview of the statistics of the peak temperatures and pressures is given in Table 4. The standard deviation indicates an average change of 6% from the base case in *T1* and 10–12% change in *P1*–*P3*, respectively.

Table 4 gives a summary of the mean peak temperatures and pressures computed for *T1* and *P1*–*P3* across all Monte Carlo runs, respectively. The standard deviations are given in brackets. It is indicated that the addition of the heat source to the uncertain parameters results in the same averaged effect in terms of the selected performance indicators (MC2 and MC3). Similarly, the random sampling and Latin Hypercube sampling of parameters result in the same means and standard deviations (MC3 and MC4). MC5 and MC6 result in overall higher and lower temperatures and pressures, respectively, related to the repository depth. The standard deviations indicating the variation around the mean are, however, similar for all repository depths. Variants MC4b and MC6b, which are based on 25% greater distance between the tunnels, result in lower means and standard deviations of peak pressures compared to MC4 and MC6, respectively.

### Risk with regard to failure of selected criteria

The risk of failure in the Monte Carlo was calculated using a FORM with failure criteria corresponding to 150°C for *T1*, and lithostatic pressure at the corresponding depth for *P1*–*P3*. A risk calculation is not performed for MC1 due to the large number of parameters (31) combined with the relatively small number of samples (200). The estimated risks for all parameter samples leading to failure in MC4 is illustrated in Figure 12. Criterion *T1* results in failure in six out of 1000 simulations in MC4. The maximum risk of failure corresponds to 0.07% (simulation #190). Criterion *P1* results in failure in 163 simulations, with a maximum failure risk of 6.6% (simulation #1). Criterion *P2* results in failure in 175 simulations, with a maximum failure risk of 6.98% (simulation #22). Criterion *P3* results in failure in 105 simulations, with a maximum failure risk of 2.24% (simulation #462).

Figure 13 shows a synthesis of maximum and minimum risk calculated with each of the Monte Carlo variants for criteria *T1*–*P3*. It is indicated that all variants result in a maximum failure risk below 10%. Variants MC2, MC3 and MC4 (600 m bgl) result in an overall similar risk of failure; variant MC5 (750 m bgl) results in a lower risk of failure due to the higher pressure values used as failure criteria. The comparison between MC4 and MC4b (40 and 50 m distance between the tunnels, respectively) shows that the maximum risk of failure is significantly reduced (i.e. from 7% to below 1%) when the distance between the tunnels is increased by 25%.

### Additional sensitivity cases for the optimization of distance between tunnels

In addition to the Monte Carlo variants MC4b and MC6b, optimization of the repository design was assessed using deterministic sensitivity cases of the 2D emplacement tunnel model based on alternative distances between the emplacement tunnels corresponding to 65, 75, 90, 110, 125 and 135% of the base-case scenario (40 m distance). For each repository depth, simulation variants were performed using the best estimates of input parameters, as well as parameter samples resulting in the maximum failure probabilities for *T1*–*P3* and parameter samples resulting in minimum pressures in *P1*–*P3*. The time histories simulated with the base-case scenario and the alternative distances between tunnels at a repository depth of 600 m bgl are shown in Figure 14. Time histories from simulation variants with decreased tunnel distance (distance scaled by 0.90, 0.75 and 0.65) are shown in green, whereas variants with increased tunnel distance (distance scaled by 1.10, 1.25 and 1.35) are shown in blue. It is indicated that decreased tunnel distances result in overall higher temperatures and a shift of peak temperature times to approximately 20 years after closure. On the other hand, increased tunnel distances result in lower temperatures at later times (i.e. after 50 years), while peak temperatures after approximately 10 years are similar to the base case. Peak pressures show a significant dependency on tunnel distance. Whereas the base case does not result in failure of any of the criteria, reducing the distance to 90% results in failure of *P1*, and near failure in *P2* and *P3*. On the other hand, increasing the tunnel distance can further decrease pressures in the emplacement tunnel and the nearfield. Peak pressure at *P1* corresponds to 92% of the lithostatic pressure in the base case, and 88, 82 and 79% of lithostatic pressure when the tunnel distance is increased by 10, 25 and 35%, respectively.

A synthesis of the results is provided in Figure 15; it shows the simulated peak pressure in *P1*–*P3* relative to lithostatic (in %) for different distances between the emplacement tunnels. It is indicated that the base-case scenarios for 600, 450 and 750 m bgl result in failure only when the distance between the tunnels is decreased below the reference distance of 40 m. The simulation variants with parameter samples that result in maximum risk of failure for 40 m distance indicate that failure is mitigated when the distance between the tunnels is increased to approximately 45–50 m. Simulation variants with parameter samples that result in the lowest peak pressure among the Monte Carlo samples allow for a significant reduction in the distance between the tunnels, even to 25 m, without causing failure of criteria *P1*–*P3*. It is thus indicated that all configurations of possible parameter samples (i.e. between the samples resulting in minimum pressure and maximum risk of failure) can be optimized using a distance between the tunnels of between 25 and 50 m to eliminate the risk of failure.

### Additional sensitivity cases for the optimization of emplacement density

The same parameter samples (i.e. base case and samples resulting in maximum failure probability) are also used in simulations with the 3D repository model to assess the impact of emplacing waste at a lower density over the entire length of emplacement tunnels, or at a higher density in 75% of the tunnel length. The following simulation cases are performed using both configurations for waste emplacement:

the base case corresponding to the best estimates of input parameters (repository depth of 700 m bgl);

a case with the input parameters resulting in the maximum probability of failure in

*P1*in the 2D emplacement tunnel model (input from Monte Carlo MC5);a case with the input parameters resulting in the maximum probability of failure in

*P2*in the 2D emplacement tunnel model (input from Monte Carlo MC5);a case with the input parameters resulting in the maximum probability of failure in

*P3*in the 2D emplacement tunnel model (input from Monte Carlo MC5).

Depending on the input parameters, the simulation of the post-closure phase using the 3D repository model requires between 24 and 48 h run in parallel on 8 Xeon^{®} E5-2630v3 @3.2 GHz processors. The time histories of temperature and pressure simulated with the base case (reference parameters) at different locations in the SF/HLW repository are shown in Figure 16. Time histories of pressure in the Opalinus Clay from the same simulation are shown in Figure 17. Thick lines show the results from the base case (emplacement in 100% of tunnel length), whereas dashed lines show the results when waste is emplaced in 75% of the tunnel length.

Time histories from the base case show peak temperatures of 124, 106 and 126°C in the SF, HLW and pilot emplacement tunnels, respectively (Fig. 16). The first peak of pressure due to thermal expansion occurs after approximately 500 years, and corresponds to approximately 9.8 and 8.0 MPa in the SF and HLW tunnels, respectively. At the same time, pressure in the pilot repository only increases to approximately hydrostatic pressure due to the faster dissipation of heat through the nearfield. The time histories of pressure in the Opalinus Clay (Fig. 17) show a peak pressure of 9.9 MPa at 30 m above the repository. The time histories of temperature and pressure in the variant with waste emplaced in 75% of the emplacement tunnels show peak temperatures of 150, 122 and 155°C in the SF, HLW and pilot emplacement tunnels, respectively (Fig. 16). This corresponds to 21, 15 and 23% increase in peak temperatures at the corresponding locations compared to the base case, respectively. The first peak of pressure due to thermal expansion corresponds to 11.6, 8.8 and 7.4 MPa in the SF, HLW and pilot emplacement tunnels, respectively. This corresponds to an increase of 18, 10 and 8% compared to the base case in the selected locations, respectively. The time histories of pressure in the Opalinus Clay (Fig. 17) indicate a peak pressure of 11.3 MPa at a distance of 30 m above the repository. This indicates an increase of 13–14% in pressures in the surrounding host rock above the emplacement tunnels compared to the base case.

The simulation results thus show that the increase in waste emplacement density in the emplacement tunnels by 33% results in an increase in peak pressures and temperatures between 10 and 20% in the base case. Figures 15 and 16 also show the failure criteria *P1*–*P3* which are similar to those used for the probabilistic analysis, adjusted to the repository depth. It is indicated that, also for increased waste emplacement density, pressures in the repository and surrounding host rock remain below the failure criteria throughout the simulation.

In Monte Carlo variant MC5, the maximum failure probability corresponds to simulation runs #222 (indicators *P1* and *P3*) and #697 (indicator *P2*) from variant MC5. A synthesis of the results obtained from the 3D repository model for these parameter samples is given in Figure 18 and shows the simulated peak pressure in *P1*–*P3* relative to the lithostatic value (in %) for different emplacement densities’ respective emplacement tunnel lengths with waste. The base case corresponds to emplacement in a 700 m length in each tunnel, whereas the alternative scenario corresponds to emplacement in a 525 m length in each tunnel.

It is indicated that parameter sets that resulted in failure of *P1*–*P3* in the Monte Carlo simulations result in peak pressures below the lithostatic pressure in the 3D repository model. This effect is discussed in detail in Papafotiou & Senger (2016*b*), and is related to the limited representation of heat- and gas-release mechanisms in the 2D representation. The 2D model is based on the assumption of superposition, which means that it does not account for heat dissipation laterally through the host rock surrounding the emplacement field. Moreover, in the 2D representation, the gas that is generated at early times from the corrosion of materials in the tunnel (liner, anchors) is assumed to accumulate in the emplacement tunnel, whereas in the 3D repository model, gas is transported along the tunnels, through the tunnel plugs and the EDZ into the operation and access tunnels. In the analyses presented here, the pressures and temperatures simulated with the 2D model are not used as quantitative evaluation measures of repository performance, but only as a comparison metric for the relative performance of the different sites. In the site selection of SGT Stage 3 additional model enhancements will be undertaken to improve the quantitative predictions of the 2D models (i.e. connect additional volumes to the emplacement tunnels, scale heat and gas generation rates).

## Conclusions

This paper presents the development of a methodology for future site-specific assessments of repository-induced effects integrating the probabilistic treatment of parametric uncertainty. The methodology emphasizes the impact of uncertainties of the model input on the response of evaluation metrics related to preliminary safety indicators, namely simulated pressures and temperatures at selected locations in the emplacement tunnels and surrounding host rock. The probabilistic workflow is demonstrated here using preliminary data, repository configurations and indicators in preparation of a criteria-based indicator approach that will support the site selection process in SGT Stage 3 based on site-specific data. The probabilistic modelling is supplemented by deterministic thermo-hydraulic (TH) calculations that demonstrate the feedback to repository optimization. The analysis considers three potential repository depths, namely 450, 600 and 750 m bgl, corresponding to the candidate sites. The reference properties of the emplacement tunnel materials, Opalinus Clay, and the over- and underlying formations were taken from existing studies corresponding to the generic repository design and the preliminary data available for the candidate siting regions.

The workflow for the probabilistic assessment of repository-induced effects is developed using three key elements:

a probabilistic sensitivity analysis to identify the most sensitive parameters with regard to temperature and pressure in the repository and surrounding host rock using the Morris Screening method;

Monte Carlo analyses used to estimate the envelope of simulated temperatures and pressures (uncertainty assessment);

estimation of risk of failure with regard to four selected indicators for temperature and pressure in the tunnel and surrounding rock.

For the 2D model used in this study with approximately 30 uncertain parameters, it has been demonstrated that Morris Screening and Monte Carlo simulations can be conducted at the cost of moderate CPU times. For large 3D models, it is recommended to evaluate alternative workflows for the uncertainty assessment (e.g. meta-models/response surfaces).

The semi-generic simulations presented in this paper demonstrate the workflow in preparation of the site-specific model analyses of repository-induced effects foreseen in SGT Stage 3. The probabilistic simulations are limited to a single-component model (emplacement tunnel) and therefore do not constitute quantitative assessments of repository performance. Heat dissipation and porewater displacement laterally into the rock surrounding the emplacement field are not taken into account. Moreover, the 2D representation does not consider gas transport along the EDZ and repository tunnels, which is a key mechanism for gas release in the repository (Papafotiou & Senger 2016*a*, 2017). Gas generation is included in the current demonstration but gas-induced effects need to be evaluated with the integration of repository-scale modelling. In such an appraisal, additional focus may be placed on two-phase properties that would affect gas pressures developing in the nearfield. Similarly, reference thermal and hydraulic properties of the host rock and their anisotropy may vary between sites and result in discriminating effects. The indicators selected for the model-based evaluation would need to be revised depending on the aims of the analysis, and may also depend on the repository configurations and site-specific data and conditions determined with the site investigations.

The demonstrated probabilistic sensitivity analysis using Morris Screening showed a clear ranking of the most sensitive parameters for each of the selected indicators. Morris Screening was implemented by running 10 sets of simulations, each set initiating from a random point of the parameter space and then perturbing each of the 31 uncertain parameters. The mean response of 10 × 32 perturbations was then synthesized in Morris plots of the average mean and standard deviation in *T1*, *P1*, *P2* and *P3*. An equal weighting of the response in all four performance indicators indicated that the most sensitive parameters are the absolute permeability, the thermal conductivity, the heat capacity and the pore compressibility of the Opalinus Clay, followed by the absolute permeability of the EDZ, the absolute permeability of the bentonite pellets and the absolute permeability of the UCU. When only indicators *P1*–*P3* are considered the sensitivity of the emplacement materials is reduced whereas sensitivity to parameters of the UCU and LCU is increased. Consequently, eight sets of Monte Carlo simulations were performed to account for different configurations and input parameters associated with the candidate sites. The Monte Carlo simulations show that output uncertainties can be constrained by bracketing: (i) the parameter ranges of the rock properties (site-specific characterization of permeability, pore compressibility and thermal conductivity of the host rock); and (ii) the properties of the engineered barrier system (engineered design of permeability and thermal conductivity of the bentonite buffer). The simulations also indicate that the probabilistic appraisal can indicate potential discriminators between candidate sites in the indicator-based evaluation depending on the input parameter reference values and distributions. Risk of failure was evaluated using FORM with failure criteria corresponding to a constant temperature for *T1* and lithostatic pressure at the corresponding depth for *P1*–*P3*. Risk of failure varied between the candidate sites depending on the reference TH properties and the preliminary failure criteria. For the site selection in SGT Stage 3, failure criteria are developed for the site-specific hydrogeological, geothermal and regional stress conditions identified in the field investigations and integrated with coupled hydro-mechanical (HM) simulations for the candidate sites. In the demonstration presented in this paper, uncertain parameters and selected indicators were mostly related to the assessment of thermally-induced effects. Deterministic model runs without gas generation (Senger *et al.* 2014) indicate that, for the assumptions made in this demonstration of the workflow, the outcome in terms of sensitivities and risk estimates related to thermal effects would, overall, be similar if gas generation was not taken into account.

In addition to the probabilistic assessments of model sensitivity and model uncertainty, deterministic simulations were dedicated to the optimization of the repository configuration with respect to thermal loading and considering the uncertainty ranges identified in the probabilistic analysis. Emphasis was on the optimization of the distance between emplacement tunnels to mitigate critical temperatures and pressures in the tunnel and surrounding rock, and on the impact of emplacement density associated with different design options for the SF/HLW repository. The analyses demonstrate that an increased footprint area could significantly reduce the risk of excessive pore pressures in the repository nearfield. A methodology for the optimization of the effective repository footprint area with respect to the performance criteria was demonstrated for the preliminary site-specific conditions using the simplified 2D model representation. Additional analyses with the 3D repository model indicate that lateral dissipation of heat through the host rock and gas release through the tunnels will result in overall lower pressures so that higher emplacement densities can be considered in the repository optimization.

## Acknowledgements

The authors would like to acknowledge Nagra for funding this study in the framework of SGT-E2.

## Funding

This research received funding from Nagra under Order no. 15’161.

- © 2018 The Author(s). Published by The Geological Society of London. All rights reserved