Coupled THM modelling of engineered barriers for the final disposal of spent nuclear fuel isolation

Abstract This paper describes the thermohydromechanical (THM) simulation of engineered barrier systems (EBS) for the final disposal of nuclear spent fuel in Finland. The bentonite barriers were simulated with the Barcelona Basic Model and the model was calibrated from laboratory tests. The evolution of gap closure and the presence of a fracture intersecting the disposal were analysed. The simulations were performed in 2D axisymmetrical geometries. Full 3D simulations were carried out in order to check the effect of the third dimension. The time required for the barriers to reach full saturation, the maximum temperature, deformations and displacements at the buffer–backfill interface and the homogenization of components both locally and globally are the main interests. The effect of rock fracture and the hydraulic conductivity of the rock are subjected to 2D sensitivity analyses.

The final disposal of nuclear spent fuel in crystalline bedrock is developing in Finland. Olkiluoto bedrock was proposed as the site for the repository (Posiva 2013). The disposal is based on the use of multiple release barriers, which ensures that the nuclear spent fuel cannot be released into the biosphere or become accessible to humans. The release barriers include the physical state of the fuel, the disposal canister, the bentonite buffer, the backfilling of the tunnels and the surrounding rock.
The KBS-3V concept (Posiva 2013) consists of deposition holes spaced between 7.5 and 11 m of each other in backfilling tunnels excavated at a depth of approximately 400 -450 m and located at approximately 25 m from each other (Fig. 1). The spent nuclear fuel will be encapsulated in final disposal canisters made of cast iron and enclosed in a copper shell. These canisters will be placed in the deposition holes and surrounded with bentonite clay, which protects the canister from any potential jolt in the bedrock and slows down the movement of water in the proximity of the canister.
In this paper, buffer (MX-80 bentonite), backfill (Friedland clay) and pellets are modelled using the Barcelona Basic Model (BBM) originally described in Alonso et al. (1990). This model is briefly described here. Elastic, isotropic, non-isothermal volumetric strains are defined by: where e is the void ratio, k i , k s and a are parameters, p ′ is net mean stress, s is suction and T is temperature. For deviatoric elastic strains, a constant Poisson ratio is used. Plasticity is accounted for using a Modified Cam-Clay yield surface (and plastic potential) with the following equation: where p s ¼ ks and p 0 correspond to the intersection of the ellipse with the p ′ -axis, q is the equivalent shear stress, and M and k are parameters. The preconsolidation mean net stress (p o ) is a function of suction in the BBM with the following form: where p c is the reference stress of the loadingcollapse curve, p * o is the initial yield mean net stress, l(0) is the slope of the virgin elastoplastic compressibility for saturated conditions of the soil and k is the slope of the unload -reload line (elastic response). The function l(s) is the slope of the virgin elastoplastic compressibility for a given suction (volumetric compressibility index), written as: where r and b are parameters. Finally, the hardening law for the BBM is given as: where d1 p v is the plastic volumetric strain increment.

Hydraulic and thermal constitutive models
The fluid flow is governed by Darcy's law, given as: where q l is the volumetric flux of liquid, k is the intrinsic permeability tensor, k rl is the phase relative permeability, m l is the viscosity of the fluid, P l is the pressure of the fluid and r l is the density of the fluid. Relative permeability is calculated as a power of the degree of saturation: where S l is the degree of saturation and m is a parameter. The degree of saturation for the liquid phase is calculated using the retention curve with the relationship of Van Genuchten (1980); this is written as: where S e is the effective degree of saturation of porous media, S l is the degree of saturation of liquid, S rl is the residual degree of saturation, S ls is the maximum degree of saturation, P g is the gas pressure, P l is the liquid pressure, l is the shape function coefficient for the retention curve and P is a parameter that can be interpreted as the pressure of air entrance. The molecular diffusion of vapour is governed by Fick's law: where i w g is the non-advective mass flux vector, f is the porosity of porous media, r g is the density of the gas phase, S g is the degree of saturation of the gas phase, D w g is the molecular diffusion coefficient for vapour in the gas phase, v w g is the mass fraction of vapour in the gas phase, I is the identity tensor and t is a tortuosity coefficient.
The diffusion coefficient of vapour is given by: where D v is the coefficient of diffusion (5.9 × 10 26 m 2 s 21 K 2n Pa) and n is a parameter (n ¼ 2.3). For the air gap element, it is convenient that retention properties evolve as the gap closes. In order to represent this response using the available functions in CODE_BRIGHT, the dependence of the retention curve on porosity has been taken into account when considering the following equation: where a is a parameter. This function allows a low retention capacity for the gap to be used when it is open (very small saturation), which evolves into a higher retention capacity for the gap as it closes (closer or similar to what occurs in the surrounding clay). The heat transfer process is governed by Fourieŕs law, given as the heat flux vector: where thermal conductivity depends on the degree of saturation in the following way: where i c is the conductive flux vector of heat, T is the temperature, l is the thermal conductivity, l sat is the thermal conductivity of the water-saturated porous medium, l dry is the thermal conductivity of the dry porous medium and S l is the degree of saturation.
For the thermal analysis, the required parameters for the materials are the thermal conductivity (l), the specific heat (c s ) and the solid density (r s ).

Material properties
Friedland clay was chosen as the reference material for manufacturing the backfill blocks (Keto et al. 2013). Although the backfill consists of different components (i.e. a foundation bed, blocks and pellets) (Fig. 2), it is treated as a single material.
It was planned that the pellets be installed between the buffer blocks and rock (Juvankoski et al. 2012), and that the pellets would occupy a volume of 10 -20% of the excavation volume associated with the deposition holes. It is also anticipated that the pre-compacted buffer blocks would swell and compress the pellet component. As a result of this, it was anticipated that the buffer blockpellet material in the tunnel should meet the performance requirements of the buffer described in Juvankoski et al. (2012).
In many aspects, it can be said that the pellet component is, in fact, a simple filler of the space between the buffer and the rock. There are several reasons for using fillers in gap. Fillers can secure better thermal conductivity and prevent spalling of the surrounding rock into the deposition hole. Also fillers might be needed in order to retain an adequate average density (Marjavaara & Kivikoski 2011).
In order to determine the hydromechanical parameters of the buffer (MX-80) and backfill (Friedland clay), an experimental programme was followed. Table 1 summarizes the initial properties of the materials for each test considered. For each material, an oedometer and an infiltration test was performed. Figure 3 shows the calibration results from the oedometer test data for Friedland clay (backfill) (Fig. 3a) and MX-80 bentonite (buffer) (Fig. 3b). As the samples were initially unsaturated, there was a stabilization (flooding) process at the beginning of the oedometer tests. It can be seen that both materials showed considerable expansive behaviour before loading. Once the saturation took place, samples were exposed to loading and unloading phases. In this study, the Barcelona Basic Model (Alonso et al. 1990) parameters used in Toprak et al. (2013) for MX-80 bentonite clay have been updated. The BBM parameters for MX-80 bentonite and Friedland clay were obtained by calibration of the oedometer tests. For the pellets, a preliminary set of parameters that take into account non-linear elasticity was established. The pellets represent a small volume fraction of the total volume of the clay barrier. The calibrated BBM parameters that were used in the 2D and 3D modelling of MX-80 bentonite, Friedland clay and the pellets are given in Table 2. Figure 4 shows the final distribution of the water content and the infiltration tests set up for Friedland clay (Fig. 4a) and MX-80 bentonite ( Fig. 4b) at the end of the experiment. The infiltration cell was made of stainless steel. The inner diameter of the infiltration cell was 50 mm and the height of the sample was approximately 63 mm. There were two sintered porous frits made of stainless-steel balls, with a diameter of 10 mm, at the top and bottom of the sample. The materials (MX-80 bentonite and Friedland clay) were placed directly inside the cells. During the tests, a target water pressure was imposed at the bottom of the sample using a pressure -volume controller manufactured by GDS (http://www.gdsinstruments.com). On the other side, the water pressure was zero (the atmospheric pressure was fixed as the origin). The tests were performed at laboratory temperature (approximately 248C). The suction of the materials was measured before the tests using a chilled mirror psychrometer WP-4 (http://www.decagon.com). The sample volume was kept constant during the test. The water inlet flow was measured continuously. At the end of the test, the sample was cut into slices with a saw. After cutting the sample, the water content was measured.
The swelling near the injection zone is possible because other parts of the sample undergo compression. Some level of suction is still measured at the end of the test. This could be due to the fact that the sample was not fully saturated, and dismantling of the sample with an associated unloading induces some suction. The hydraulic parameters of MX-80 bentonite and Friedland clay were calibrated from these infiltration tests and are listed in Table 3.
The host rock (Hagros et al. 2003) and the canister (Raiko 2012) are considered to be linear elastic with parameters E, n and a, respectively, for Young's (elastic) modulus, Poisson's ratio and the coefficient of thermal expansion. Hydraulic and mechanical properties of the host rock and canister are given in Tables 3 and 4, together with the properties of the gap.
Mechanical parameters for the air gap were proposed in Toprak et al. (2013). The 10 mm air gap

Governing equations
The equation of equilibrium of forces (or stresses), the equation of mass balance of water and the equations of mass balance of energy were solved   simultaneously to obtain displacements, pressure and temperature at each point and time of the model. The formulation can be found elsewhere: for instance, in Olivella et al. (1994Olivella et al. ( , 1996.

Modelling of vertical disposal schemes
The model geometry and materials are shown in Figure 5a. This figure shows boundary conditions that are both hydraulic (Fig. 5b) and mechanical (Fig. 5c). The analysis assumed under axisymmetrical conditions. The initial water pressure for all materials in the deposition hole was 241 MPa. The initial temperature was 10.58C throughout the domain modelled. The rock had a hydrostatic water pressure (with P l ¼ 4.36 MPa at z ¼ 436 m and P l ¼ 3.98 MPa at z ¼ 398 m, where z is the depth below the surface), which was maintained along with corresponding boundary conditions. The temperature at the boundaries is seen to evolve with time (Fig. 6a). This temperature can be calculated from an analytical solution that takes into account the presence of all canisters (Ikonen 2003). An initial confining stress of 10.6 MPa was considered for the host rock. This confining pressure was also used as a boundary condition (applied on the top boundary).The lateral and bottom boundaries had prescribed normal displacements. During the excavation process (assumed to be 1 year), a prescribed liquid pressure (25 MPa) was applied to the excavation surface to represent the process of ventilation. This boundary condition was removed when the buffer and backfill materials were emplaced. During the simulated time period, the hydrostatic liquid pressure was imposed on the top and bottom boundaries.
Regarding the power evolution of the canister, there are two main parameters playing a fundamental role, these are the residual power at the time of deposition and the decay rate. The work by Hökmark et al. (2009) was taken as a reference for the calculation of the power and the decay heat rate. The power as a function of time for an individual canister can be expressed as: In this expression, P(0) is the canister power at the time of deposition, and a i, and t i are parameters. Two parameter sets from SKB of power data relative to the time of cooling prior to disposal are given in Table 6. A burnup of 38 MWd/kg U (megawatt days per kg of uranium) is used as a reference. The coefficients given in Table 6 are valid for an initial power level of 1837.3 W (in the case of 30 year-old fuel) and an initial power level of 1545.3 W (in the case of 40 year-old fuel). The work presently being performed is targeting a 1700 W initial power level at the time of deposition. The power function and the prescribed temperatures used in the models presented here are shown in Figure 6. Values of 1700 W of canister power, 25 m of tunnel spacing and 11 m of canister spacing, and variable temperatures at the boundaries (Fig. 6a, b) were used to perform the thermohydromechanical (THM) calculations (see Toprak et al. 2013 for detailed information on thermal calculations under axisymmetrical conditions).
With the formulation and parameters decided, a series of THM calculations were carried out using a 2D axisymmetrical configuration, as well as using a 3D configuration. For the 2D axisymmetrical configuration, three cases described in Keto et al. (2013) were been considered. Each one has been given a name depending on the hydraulic properties chosen for the host rock, and one case includes a higher permeability zone representing a fracture. The cases are: † Normal deposition hole (Normal_Case). The rock hydraulic conductivity is 1.52 × 10 212 m s 21 (equivalent to 1.52 × 10 219 m 2 of intrinsic permeability). Assuming water density of r ¼ 1000 kg m 23 , water viscosity m ¼ 0.001 Pa s and gravity g ¼ 10 N kg 21 , gives the following conversion: K ¼ (rg/m)k ¼ 1000 × 10/0.001 k: that is, K (m s 21 ) ¼ 10 7 k (m 2 ), where K is hydraulic conductivity and k is intrinsic permeability. The rock is considered  homogeneous: that is, no fracture is considered. This value of intrinsic permeability is considered to be a standard value for the rock at the repository conditions. † Wet deposition hole (Wet_Case). In this case, the rock hydraulic conductivity is 1.52 × 10 212 m s 21 (equivalent to 1.52 × 10 219 m 2 of intrinsic permeability) and a predefined fracture was considered with a transmissivity of 1.1 × 10 29 m 2 s 21 . The fracture (high conductive zone) is horizontal because for the axisymmetrical geometry, it is not possible to represent inclined surfaces. The fracture thickness considered is 8 cm in order to avoid elements that are too  small. The fracture is modelled with continuum elements. The equivalent hydraulic conductivity of the fracture is 1.37 × 10 28 m s 21 (1.37 × 10 215 m 2 of intrinsic permeability). † Dry deposition hole (Dry_Case). The rock hydraulic conductivity is 3.53 × 10 213 m s 21 (equivalent to 3.53 × 10 220 m 2 of intrinsic permeability). In order to get drier conditions, the intrinsic permeability of the rock is reduced by approximately one order of magnitude, compared to the so-called normal case. Figure 7 shows the hydraulic boundary conditions and intrinsic permeability for the three cases. In the case of the wet deposition tunnel, there is a zone with higher conductivity representing a fracture, and an appropriate boundary condition is applied at the end of the fracture. The temperature of water inflow through the fracture is the temperature in the lateral boundary (no heat-flow boundary condition). As indicated in Figure 7, the rock is less permeable in the dry deposition case compared to the wet and normal cases. There are three critical zones that are considered relevant in this THM model, these are: (1) 'the buffer-backfill interface' (Å kesson et al. 2010;Leoni 2013), where the vertical displacements show the largest values; (2) 'canister wall in contact with the buffer', the maximum temperature of the engineered barrier is reached on this surface (Ikonen 2003), and desaturation of the buffer takes place because of intense heating from the canister; and (3) 'buffer blocks under the canister', the maximum stresses are achieved here owing to the combination of the canister weight and swelling of the clay materials.
Results for 2D axisymmetrical modelling Figure 8 shows the temperature evolution in the canister. The temperature increases due to heat generation in the canister, which, as shown in Figure 6, is characterized by decay to half-value in about 70 -80 years. For the geometry considered, a peak value above 808C takes place and the temperature decreases afterwards. However, the shape of the temperature evolution in the canister is influenced by the water saturation of the different components of the barrier and the presence of the gap. Both heat capacity and thermal conductivity depend on the degree of saturation.
The maximum temperature is not significantly influenced by the hydration conditions in general, but the evolution at earlier times is different due to the effect of the gap and saturation conditions of the barrier components. The gap under dry conditions produces thermal isolation as it remains open (air has lower heat transport capacity than water or soil). When the rock intrinsic permeability is very low (Dry_Case), the gap closure is delayed and the temperature in the canister increases due to the isolation effect. For the Dry_Case, the maximum temperature in the canister is higher, by a few degrees compared to the other cases and occurs in a different peak that occurs earlier. The peak that occurs around 1 year is caused by the presence of the gap, so the temperature decrease after that peak is motivated by the increase in conductivity induced by the combined effect of saturation and closure. The peak caused by the gap closure occurs earlier than the peak caused by the decaying power function. Figure 9 shows the evolution of the liquid pressure for the three cases in the bentonite ring adjacent to the canister. In the case of the wet deposition hole, as the global intrinsic permeability (including the fracture) is higher compared with the other cases, the saturation of materials takes place faster. Desaturation of the bentonite ring is much greater in the Dry_Case because the rock intrinsic permeability is lower than in the other cases, and this implies a lower thermal conductivity of the dried-up buffer (Börgesson et al. 1994;Tang et al. 2007;. With regard to stress development, the most critical zone is below the canister where the stresses become greater. Figure 10 shows the swelling pressure of buffer blocks under the canister for the three cases. Stress increases as liquid pressure evolves from extremely negative values caused by heating towards zero and positive values (i.e. hydration produced by a water supply from the host rock). The effect of the fracture is clearly observed in Figure 10,   as the section where stress is plotted is close to the fracture. In the case of high intrinsic permeability of rock, water flows from the rock to the buffer in a much faster way. The time to reach a stabilized stress and the values achieved are different between the three cases. A lower intrinsic permeability of rock delays the process of swelling and, hence, the effective stress development. A larger effective stress is developed below the canister as saturation takes place faster. This can be explained by the fact that backfill is more rigid during buffer swelling in the Wet_Case compared to the Normal_Case. The same behaviour is observed if the Normal_Case is compared with the Dry_Case. Figure 11 shows the swelling pressure of the backfill (central zone) for the three cases. As this zone is far away from the fracture, the time to achieve the maximum swelling pressure of the central backfill mainly depends on the intrinsic permeability of the rock and not on the water entry from the fracture. In the Wet_Case, the swelling pressure of the central backfill reaches a value of 7.5 MPa and stabilizes at 6 MPa. Except for the Dry_Case, the other two cases reach a peak value of mean effective stress at the same time. The reason for this is that the intrinsic permeability of the rock is of greater importance with respect to the time required to reach full saturation of the backfill. However, the swelling of the buffer blocks and the ring is greater in the Wet_Case. As a result, the achieved stress in the backfill is greater in the Wet_Case, despite it taking the same time to reach full saturation as the Normal_Case.
The peak that appears in the evolution of effective stress is caused by the effect of unloading after saturation. Once saturation is reached, the maximum effective pressure is achieved. Then, pore pressure becomes higher than atmospheric pressure and, as the material is saturated, effective stresses come into play. The increase in pore pressure implies a decrease in effective stresses. This is represented in Figure 12, which shows the suction -mean effective stress path. Although suction is a positive definite variable, the negative values have been represented in order to highlight the fact that effective stress decreases after full saturation of the medium. Figure 13 shows the evolution of vertical displacements at the buffer-backfill interface for the three cases. The vertical displacements are generated by the swelling of the buffer blocks. Displacements are mainly concentrated at the buffer-backfill interface. Owing to the swelling, the buffer blocks push the backfill up. In the case of the dry deposition hole (Dry_Case), as the swelling of buffer blocks takes place slowly, the achieved displacements at the interface show a lower value. In other words, the buffer and backfill swell in a more simultaneous way for the Dry_Case compared to the Normal_Case and Wet_Case for which the buffer hydrates more rapidly than the backfill.
As a conclusion, it can be said that owing to the delay in water supply in the Dry_Case, displacements take place later compared to the other two cases. The maximum buffer-backfill interface displacement peak is about 5 -6 cm in the Dry_Case. In the Wet_Case, vertical displacements reach to 9-10 cm, and for the Normal_Case they are the range of 7-8 cm. Table 7 summarizes the results in terms of the maximum temperatures achieved, times to reach saturation, saturated density of the buffer and backfill, stress development, and maximum displacements. The Dry_Case curve for temperature shows an early peak (absolute maximum) that is caused by the lower thermal conductivity of the open gap, which remains open for a longer time compared to the other cases. Earlier saturation is motivated by the presence of the highly conductive zone that simulates a fracture. Final densities show little variation. As the global volumetric deformation of the engineered barrier is nearly zero (except for the deformation of the rock, which is very small), mean effective stress development is of the same order of magnitude as the swelling pressure. There are differences, however, as the saturation of the engineered barrier takes place in different ways depending on the conditions of hydration (intrinsic permeability and fractures, and the shape of the buffer-backfill system). A greater development of stresses is calculated with the Wet_Case. Finally, displacement at the interface is influenced by the earlier or later expansion of the buffer.

Results for 3D modelling
In this section, some results for 3D modelling are presented in order to check the influence of the simplification performed in the axisymmetrical models. For simplicity, the 3D geometry does not take into account the presence of either the pellets or the air-gap element. Besides this, backfill tunnel has a greater volume compared to the 2D calculations. The BBM parameters used for the buffer and backfill are the same as in the previous 2D Fig. 12 Suction v. mean effective stress path for a point below the canister and a point in the backfill. Negative suctions are included to show that after saturation the mean effective stress decreases due to the development positive pore pressure. axisymmetrical calculations. The value of the intrinsic permeability for the rock corresponds roughly to the Dry_Case (1.52 × 10 220 m 2 ), the intrinsic permeability of the backfill is the same as in the previous cases, and the intrinsic permeability of the buffer is higher but is in the range of typical values for bentonite. So, the values chosen fall in the same order of magnitude of intrinsic permeability and this is favourable for the numerical solution, especially for the 3D model with a relatively simple mesh. Further work on mesh optimization is required in order to obtain 3D solutions of the THM problem, with barriers having greater contrasts in permeability. Figure 14 shows a comparison between the 2D and the 3D case in terms of the total mean stress. As the backfill tunnel has a greater volume for the 3D geometry, the whole system needs more time to reach full saturation. Therefore, the generated stresses reach the steady-state conditions in faster times in the case of 2D calculations. The values achieved for the total mean stresses are in the same range (i.e. from 8 to 12 MPa).

Conclusions
The laboratory tests (oedometer and infiltration) performed by B + Tech Oy (Finland) were modelled using the finite-element code Code_Bright. The Barcelona Basic Model (BBM) was used to model the mechanical constitutive behaviour of the MX-80 bentonite and Friedland clay, together with Darcy's law for the flow. The parameters obtained from the modelling of laboratory tests were used to perform 2D and 3D sensitivity analyses for the final disposal of nuclear spent fuel. The values in the Table 7 correspond to representative points for the axisymmetrical cases. In the Dry_Case, the gap between the canister and the buffer closes later compared to the other cases. The maximum temperature in the canister reaches a value of 848C for the Dry_Case (compared to 82.58C), which is a clear effect of the air-gap element. As discussed throughout this paper, the fractures and hydraulic conductivity of rock have an impact on the results obtained in terms of maximum temperatures, swelling stress and displacements.
The buffer-backfill interface has been analysed and the maximum vertical displacement has the same order of magnitude as the one calculated previously by Leoni (2013).
With respect to the 3D calculations, the system reaches steady-state conditions later than in the 2D axisymmetrical models, mainly because the backfill volume is larger. However, the values obtained for stresses in the 3D models are in the same range as in the 2D calculations.