Bentonite clay barriers in nuclear waste repositories

Deep disposal of high-level radioactive waste is the preferred solution worldwide for the long-term disposal of nuclear waste. This concept involves a series of geological and engineered barriers that provide isolation of the waste from the biosphere. Most designs involve bentonite clays as seals in different forms. During the operation of the repository, the bentonite will be subjected to a series of complex thermo-hydro-mechanical phenomena that will interact with each other. Predicting the long-term safety of geological repositories thus involve a rigorous analysis of these multi-physical processes. This paper presents a review of recent numerical approaches and analyses that have aimed to improve the understanding of processes that will take place in clay barriers over the lifetime of nuclear waste repositories. The understanding of bentonite behavior from laboratory experiments under relevant conditions is analyzed. Constitutive models that attempt to predict such behavior are presented, focusing on the stress-strain model ACMEG-TS. These models are implemented in the finite element code Lagamine which allows for the study of real scale tests. Two application cases are presented: the performance of a clay barrier according to the Swiss design, and a model of the FEEBX in situ experiment, which was modelled after a real repository under natural conditions. Overall, the relevant processes are well captured quantitatively by the models, allowing for the establishment of sound basis for future prediction and long-term design of the final underground repositories.


Introduction
According to the World Nuclear Association, in 2017, 10.5% of total electricity worldwide was produced by nuclear energy. This source of energy requires the management of spent fuel with long-lived radionuclides that might be radioactive for 100k to 1M years. It is therefore necessary to envisage solutions for the longterm isolation from the biosphere. Deep geological disposal is currently the most feasible and safest option for the long-term disposal of high-level spent fuel. Different designs of geological repositories have been proposed by different countries (see e.g. [22] for a review of designs). Because of its favorable properties, such as low permeability and self-healing capacity, most of these designs involve the use of bentonite clays as a sealing material.
Predictive analyses of the long-term behavior of such a disposal method need to be based on a sound understanding of the physical processes that are involved throughout the lifetime of the repository. In this context, the highly coupled thermo-hydro-mechanical (THM) phenomena that occur in engineered and geological barriers must be captured adequately by predictive models. Both laboratory-(elementary scale) and fieldscale tests (repository scale) have been extensively performed during the last decades. These experiments have greatly improved the knowledge of the interplay between phenomena of different physical natures and the bentonite behavior under these conditions. These observations have allowed for the development and validation of numerical tools with the ultimate goal of performing predictions about the long-term performance of deep geological repositories.
This paper presents a review of numerical analyses that aimed to improve the understanding of bentonite behavior in clay barriers. The analyses are based on a coupled THM constitutive framework, based on the behavior of bentonite evaluated at the laboratory scale (elementary scale) that, implemented in a finite element code, allowed for the identification of coupling between different physical phenomena at a large scale. The paper starts with a short introduction to the main properties of bentonite and the experimental evidence upon which the constitutive framework is based. Afterwards, the constitutive model ACMEG-TS, that provides a framework for the THM behavior of clays, is introduced. The model is implemented in the finite element code Lagamine [3,5,6], and has been used to analyze and interpret relevant experiments at both small and large scales. First, the simulation results of a small-scale laboratory test that had the objectives of providing a benchmark for the numerical modelling of bentonite swelling and homogenization are presented. Two models representative of real repositories are then introduced. The first case corresponds to the prediction of canister movement in the Swiss design of clay barrier. The second case considered the analysis of the large-scale test FEBEX, which provided a wealth of experimental data, and involved a blind prediction of the final state of the barrier after 18 years of operation.

Multi-physical behaviour of bentonite
Bentonites are clays with large amounts of smectite minerals. The smectite mineral is characterized by a basic structure that by itself is not electrically neutral and attracts exchangeable cations such as Na, Ca or Mg, that ensure electro-neutrality. Depending on the environmental conditions, several smectite particles may be stacked together [18]. The interlayer, defined as the space between two consecutive smectite layers, is arguably the differential factor with other clays, being able to adsorb up to three or four layers of water molecules due to the interaction with the exchangeable cations and the smectite surface forces. It has been recognized that up to 60% of pore-water in bentonite might be adsorbed water [16].
The water adsorption process involves either large swelling strains or large swelling pressures, depending on the confinement conditions. The main safety properties of bentonite for a deep geological repository are its high swelling pressure, which provides self-healing capacities, and its low hydraulic conductivity combined with a high solute retention potential. Figure 1 shows the dependency of swelling pressure on the dry density (equivalently void ratio). The target dry density of bentonite in a clay barrier is therefore a key value of the design, such that a given range of hydraulic conductivity and swelling pressure upon saturation is achieved.
Several forms of bentonite have been proposed in order to fill the space between the canisters and the tunnel, the most studied being compacted blocks and bentonite pellets of high density. For instance, the current Swiss design envisages the use of bentonite blocks to support the steel canisters containing the radioactive waste, and the use of pellets to fill the remaining space inside the tunnel. A large scale demonstration project is currently in operation in the underground laboratory of Mont Terri (northern Switzerland) in order to demonstrate the feasibility of such a concept [13].  Figure 2 presents a section type representative of the Nagra concept of nuclear waste repository. The figure is used to illustrate the main phenomena that are expected in a repository, in particular on the multi-physical actions to which the bentonite will be subjected to. In the contact with the canister, the bentonite will be subjected to the heat emitted by the canister that might reach 150 degrees Celsius in early stages of disposal. Simultaneously, the bentonite will be progressively hydrated due to the inflow from the host-rock. Consequently, the outer part will tend to swell, compressing the inner part that will contract due to the increase of suction induced by the increase of temperature. Therefore, in order to predict how the bentonite will behave in these conditions, the characterization of the material has to include its response under relevant changes of stress, water saturation and temperature.
One of the key properties that need to be evaluated is the water retention behavior. A precise determination of the water retention behavior is a prerequisite in order to perform quantitative predictions of the evolution of the state of the barrier. It significantly affects the hydraulic conductivity, thermal conductivity and the effective stress. Therefore, the time needed to reach equilibrium and the final distribution of dry density are dependent on the water retention curve as well. Seiphoori et al. [20] investigated the water retention behavior of granular MX80 bentonite. Their results are shown in Figure 3, showing that, after saturation was reached, a subsequent drying showed significant deviation with respect to the wetting path. Moreover, the second wetting path did not coincide with the first one. These findings were related to the modification of the fabric after the saturation. Indeed, MIP results showed that a significant change of the pore size distribution took place due to the homogenization of the initial structure formed by the highly compacted grains. . Water retention behaviour of MX80 granular bentonite subjected to a wetting-drying cycles (data from [20]).
The coupling between the hydraulic and mechanical response is apparent on the development of high swelling pressure when the bentonite is wetted under confined conditions, or the development of high swelling strains when it is wetted without kinematical restrains. Water retention properties strongly depend on the dry density. Likewise, the compressibility is dependent on the suction to which bentonite is subjected. Figure 4a shows the oedometric response upon compression of the granular MX80 bentonite after equilibration at different values of suction [21].  The compressibility is shown to be highly dependent on the suction value. Figure 4b shows the high void ratio that is attained upon saturation under low stress conditions [21]. An advanced triaxial cell allowing for simultaneous control of temperature and suction was used in [21] to study the response of bentonite to thermo-hydromechanical actions. Figure 5 shows the effect of a temperature increase on the deformation of granular bentonite samples subjected to a controlled total suction of 20 MPa under different confining stress. As a response to temperature increase, the bentonite contracted under high stress and it expanded under low stress.

Stress-strain constitutive model ACMEG-TS
Francois and Laloui [10] presented a constitutive model for the THM behavior of clays, ACMEG-TS, which accounts for the aforementioned coupling between the mechanical, hydraulic and thermal responses. The model is formulated using an elasto-plastic framework based on the Hujeux's critical state model [11]. The strain tensor is divided into an elastic part (reversible) and a plastic part (irreversible) as: where is the strain tensor and superscripts e and p stand for elastic and plastic strains respectively. A generalized effective stress and suction are used as constitutive variables [14]. The generalized effective stress is based on Bishop's expression (Bishop 1959) where the pore pressure averaging coefficient is taken equal to the degree of saturation [19]: where is the effective stress tensor, is the total stress tensor, is the degree of saturation, is water pressure; is air pressure; = − is suction and = − is the net stress tensor. The use of such effective stress implies several hydro-mechanical coupling mechanisms, a detailed description of which can 10 100 be found in Nuth and Laloui [14]. Changes in effective stress are solely related to changes in elastic strains as where and is the elastic stiffness tensor, which has the following form: where is the shear modulus and the bulk modulus. These moduli depend on the confining stress = tr( )/3 as where and are the shear and bulk modulus respectively at a reference pressure and is a parameter controlling the elastic nonlinear behaviour. In non-isothermal conditions, additional elastic strains are induced by temperature changes according to where is given by where is the thermal expansion coefficient at a reference temperature and is the initial critical state pressure at the reference temperature. Effective stress paths on the yield surface will result in plastic strains. These might be induced by two mechanisms, namely isotropic d and deviatoric d plastic strains. Thus each mechanism corresponds to a different yield surface whose expressions are, respectively, where is the preconsolidation pressure, is the deviatoric stress, is the slope of the critical state line in the ( − ) plane; is a material parameter defining the shape of the deviatoric yield surface; = with the critical state pressure; and correspond to the radius of pure elastic domain for each mechanism, allowing for plastic strains within the yield limits. It can be readily observed that both mechanisms are coupled by means of the preconsolidation pressure, . depends on the volumetric plastic strain, temperature and suction as: where is the plastic compressibility modulus (slope of the − ln relation); is the air entry suction value; and and are material parameters accounting for the increase/decrease of elastic domain with suction and temperature. Figure 6 shows graphically how temperature and suction influence the size of the yield surface, decreasing when temperature increases and increasing when suction increases above the air entry suction.
The flow rule of the isotropic mechanism is associated whereas that for the deviatoric mechanism can be nonassociated. Thus, the plastic potentials have the following form: where is a non-associativity parameter. The magnitude of plastic strains depends on the derivatives of these potentials as follows where and are the plastic multipliers which are determined following Prager's consistency condition extended to multi-dissipative materials [17]. In order to account for unsaturated conditions a relationship between degree of saturation and suction is introduced [10]. An elastoplastic approach is used to model water retention. This approach includes hysteretic behavior controlled by two yield surfaces corresponding to drying and wetting paths. The conceptual model is represented in Figure 7. A change in is induced when suction reaches either the wetting limit or the drying limit, satisfying respectively: where and stand for the yield activated during a wetting or a drying process respectively, is the drying yield limit which, during a desaturation/saturation process (i.e. = 0 or = 0), remains equal to the actual value of ; and a material parameter controlling the water retention hysteresis.
If the initial state is saturated, is equal to the air entry suction and increases when suction exceeds this value according to the following hardening law

= exp(− ∆ )
( 1 7 ) where is the slope of the retention curve in the ( − ln ) plane (see figure 7). The same process is activated for a drying path in the opposite way (expression (16)) Thus, expression (17) describes the hardening process that leads to changes in . To account for the dependency of water retention on dry density and temperature, is taken as a function that depends on the material state as where is the initial air entry suction, and and describe the evolution of air entry suction with temperature and volumetric strain respectively.

Constitutive equations for water, air and heat flow
The balance equations of mass, energy and momentum are based on the compositional approach and are described in detail in Collin et al. [4]. For the sake of conciseness only the most relevant constitutive relationships will be described here.
Water flow is modelled by means of Darcy's law neglecting the gravitational forces: where is the vector of water flux, is the tensor of intrinsic permeability, is the relative permeability, is water viscosity, is water pressure. Relative permeability evolves with the degree of saturation, following an exponential law = (20) where is the permeability at saturated state, and is a material parameter. In the present case it will be considered that the permeability tensor is isotropic, i.e.: = The influence of deformation on the intrinsic permeability is taken into account by means of the Kozeny-Karman formula: where , is the initial intrinsic permeability, stands for porosity, is the initial porosity and and are material parameters.
The effect of temperature on water is important when considering vaporization. Vapor in the porous medium is supposed to be in thermodynamic equilibrium with liquid water, and, using Kelvin-Laplace's law as the definition of relative humidity h, the following relationship is obtained: where is the molar mass of vapor, is the gas constant of water vapour, and , is the saturated vapor density, that is itself dependent on temperature. The overall air density is thus: This relationship is used in the vapor diffusion law that is based on Fick's law in a porous medium: where is the vapor flow, is the diffusion coefficient and τ the tortuosity. Heat transport is governed by both conduction and convection: = −Γgrad( ) where λ is the thermal conductivity of the mixture. Depending on the properties of each material, this physical characteristic is either considered as a function of the volume ratios of solid, liquid water and gas phases, or a specific function for the material.

Models of laboratory experiments: homogenisation after gap filling
Dueck et al. [7] performed a series of tests to analyze the response of bentonite upon wetting under different scenarios. The aim of these tests was to provide a benchmark for numerical modelling and validation. One of the tests was modelled with AMCEG-TS and is reported here. The test consisted in three phases: wetting under constant volume conditions, swelling to a height of 25 mm and equilibrium with external water under the final volume. The phases are shown graphically in Figure 8. The swelling pressure simulated for the test 1a01 is shown in Figure 9 together with the experimental results. The calibration was performed to fit the swelling pressure both radially and axially at the equilibrium state of the first stage (saturation stage). The remaining stages, that is, unloading, swelling, and subsequent swelling pressure, are predictions, since no further calibration of parameters was attempted. Overall model performance is satisfactory, with a qualitative reproduction of the measured evolution of swelling pressure. After the saturation phase, subsequent unloading and swelling stages are well reproduced in terms of axial stress.  The results of the simulations for the test 1a02 are shown in Figure 10 in terms of swelling pressure. The particularity of this test is that swelling pressure was monitored in the radial direction at three different levels. Both initial and final values of the position at which stresses are obtained numerically are shown in the inserts of Figure 10.
Due to the large strains involved in the problem, a direct comparison between the measurements (obtained at a fixed point) and the numerical results (obtained at points that evolve in height) is not straightforward. In Figure 10 the measurements obtained at a height of z=30 mm are compared to two integration points that were located at z=30mm at a different time. One was located at z=30 at the beginning, while the other was located at z=30 mm at the end of the simulation (after swelling). Thus, these two integration points provide the envelope of stresses developed at z=30 throughout the test.
In terms of swelling pressure, the results are in good agreement both qualitatively and quantitatively with the experimental values. This demonstrates the predictive capabilities of the model provided the results are obtained with the very same set of parameters calibrated with the saturation stage of test 1a01. Radial swelling pressures are well predicted in the bottom parts of the sample. The only noticeable difference is in the rate of pressure development at the beginning at z = 15 mm which is slower in the model compared to that measured in the experiments.  Figure 11 shows the distribution of radial stress in the sample as the gap fills and homogenizes. Once equilibrium is reached, between (1000-1600 hours) a nonhomogeneous stress state is found in the sample with higher radial stresses located at the bottom of the sample.

Analysis of canister movements
One of the questions related to the EBS (Engineering Barrier System) design is linked to the positioning of the canister in the center of the tunnel and the evolution of this position, especially during the saturation process. A change in position could lead to an inhomogeneous distribution of stresses as well as temperature, leading to a less efficient buffer and possible damage of the canister. Dupray et al. [8] showed how a non-symmetrical and heterogeneous (granular and block bentonite) engineering barrier system could possibly behave. The key issues are related to the possible different reactions of the granular bentonite versus the blocks when submitted to hydration and heating/cooling cycles. Dupray et al. [8] thus investigated the ability of ACMEG-TS to simulate the expected behavior of an EBS design in terms of temperature, hydraulic and mechanical evolution.  The setup of the model was taken from the current Swiss concept of nuclear repository which entails a block pedestal that ensures enough strength to support the canisters and a backfilling of bentonite pellets. Figure 12 shows the finite element mesh that was used to perform the analyses.
Comparisons with previous studies of the same case using different sets of physical processes highlighted the validity of the assumptions adopted as well as the differences due to the mechanical couplings of thermal and hydraulic behavior. The simulated swelling pressures corresponded well to those measured in the laboratory at the same densities, although they can be treated as maximum values due to the construction gaps that will exist in real cases and allow a certain amount of free swelling. The generalized effective stress framework associated with the advanced water retention behavior seems well adapted to predictions of in situ swelling behavior. This study shows that the movements of the canister are primarily governed by the swelling pressures of the two materials and, in this case, are dominated by the highest swelling pressure created by the bentonite blocks. The thermal dilation of the materials also creates a slight uplift. As a result, a heave motion is obtained. This heave was also observed in a real-scale scale experiment EBS [12]. The maximal displacement of 15 mm does not affect the heat and water flow processes and only reduces the thickness of the buffer by 2 %.

Analysis of the FEBEX test
The FEBEX in situ test was a near-to-real experiment of nuclear waste disposal in a granitic formation. The objectives of the experiment were to study the technical feasibility of such a design and to obtain data from real scale test in order to better understand the processes. A layout of the experiment is shown in Figure 14.
Dupray et al. [9], modelled the THM behavior of the buffer, which was made of FEBEX bentonite, and the surrounding host rock (granite), by means of the finite element method with Lagamine. The ACMEG-TS was used to model the constitutive behavior of bentonite. The material parameters of FEBEX bentonite were obtained from numerical simulations of a series of oedometric tests under controlled temperature and suction. The parameters for the thermal and hydraulic diffusion models were extracted from the literature.  Despite the non-simulated construction gaps, the results of the simulations showed good agreement with the experiment. The confined swelling behavior of highly compacted bentonite, including three distinct regimes, was well reproduced. It is an important feature of the ACMEG-TS model that results directly from the introduction of a wetting collapse mechanism in the framework of generalized effective stress. Moreover, the results obtained were interpreted in the light of elastothermoplasticity of unsaturated materials. In particular, the effects of heating and wetting were highlighted, and the spatial and time distributions of their respective influence on the volumetric behavior were shown. The influence of the water retention behavior, which includes the effects of hysteresis, temperature, and mechanical strain, was shown in terms of suction and mechanical effects.
On the basis of these results, two major features of the behavior of compacted bentonite were highlighted. The first aspect is the hysteretic water retention behavior of bentonite, which, coupled with adequate modelling, allows for a proper reproduction of the evolution of relative humidity in the buffer. The second aspect is the complexity of the swelling behavior of bentonite, which was analyzed using the framework of suction-induced plasticity. This allows for a better understanding of experimentally-observed behavior while being a step toward more precise modelling. In such a way, the use of an advanced thermoplastic constitutive model using an unsaturated formulation significantly advances the knowledge of the highly coupled processes occurring in a clayey formation, initially unsaturated, in the near field of a heat-emitting radioactive waste.
A further analysis of the FEBEX test was reported in ref. [15]. In that case a blind prediction of the final state of the barrier was made before the post-mortem analysis. The model featured analyses up to the dismantling phase, including 11 stages of construction and corresponding dismantling activities. Figure 16 shows the results for two sections of temperature and relative humidity. The cooling stage, resulting from switching off the first heater is well captured in the simulation. Likewise, the vaporization of water close to the heater, leading to a reduction in relative humidity in that zone, is well reproduced by the model, closely following the measurements.   Figure 17 shows the predicted values of dry density and water content compared to the values that were obtained from the post mortem analysis in the section between the two heaters [23]. Overall, good agreement was achieved in most of the sections analyzed.

Conclusions
The long-term safety analyses of deep geological repositories for nuclear waste require advanced models for the geomechanical processes that will take place. These require in-depth understanding of multi-physical couplings. This paper provided a summary of thermohydro-mechanical numerical models for clay barriers in nuclear waste disposal focusing on the role of the thermohydro-mechanical constitutive model. It is concluded that most of the thermo-hydro-mechanical processes that are involved in the performance of clay barriers are well understood, providing confidence in long-term evaluations of safety performance. Nevertheless, despite the considerable advances that have taken place in the basic understanding of the relevant physical phenomena in clay barriers, this is still a topic of active research. Current developments are under way in order to enhance the capabilities of constitutive models, such as including adsorption and capillary processes into the water retention formulation [1].