Heating-induced creep and potential creep rupture of clay liners for nuclear waste repository

The aim of this study is to assess the potential of encountering a heating-induced creep rapture of clay liners in nuclear waste repository. Groundwater and soil contaminations may occur if the elevated temperatures, expected in the vicinity of nuclear waste repository, trigger creep rapture of the clay liners. In this study, we utilize simulations based on the discrete element method (DEM) to understand the conditions under which heating-induced creep rupture can take place. In lieu of the conventional local/non-local damping mechanism usually utilized in DEM simulations to dissipate energy, the DEM simulations presented in this study incorporate the rate process theory as a damping mechanism to model soil creep. The results of a base anisotropic model at 70 °C show a dramatic increase in the creep rate at high temperatures showing creep rupture. Such undesired behavior can be mitigated by engineering clay liner materials to sustain and resist the expected high temperatures expected around nuclear waste repository.


Introduction
Interests to study temperature effects on the deformations and strength of soils have gained a wide attention in the last few decades within the engineering community. Many natural phenomena (e.g., extreme temperature cycles) and energy geo-structures (e.g., heat exchanger piles) involve thermal loading of soils, which in turn impacts soil behavior [1][2][3][4][22][23]. Adequate design of such geo-structures is, therefore, crucial to ensure their safety and resiliency.
Owing to their small permeability, clays are usually used as liners in nuclear waste repositories to prevent leakage. The overburden pressures experienced by these clay liners may be accompanied with elevated temperatures over the lifespan of the radioactive buried wastes. Clay liners may, therefore, be subjected to a sustained thermo-mechanical loading. Under this loading condition, clay liners experience accelerated creep. The creep magnitude and rate are functions of both the applied pressure and temperature [5]. Under a certain pressure and temperature, the continuous creep may lead to a creep rupture [6]. Owing to the particulate mechanism of soil creep, the conditions upon which a creep rupture under elevated temperature occurs are not fully understood or constrained.
Therefore, the aim of this paper is to investigate the potential occurrence of creep rupture of clay liners at elevated temperatures. A discrete element model based on the rate process theory is employed here to model the creep phenomena. The rate process theory is first introduced briefly and factors affecting creep rate and its magnitudes are presented thereafter. Results of two simulations at two different temperatures are then presented and discussed in light of the rate process theory and creep phenomena.

Background
Soils undergo continuous deformations under sustained loads and stresses (e.g., mechanical, thermal, or both), a phenomenon commonly known as creep. The creep process in any soil can be divided into three stages ( Fig.  1). High creep rates occur at early times defining the primary creep stage. As time passes, the rate of creep slows down to a constant rate introducing the secondary creep stage. Under certain boundary conditions, the tertiary stage of the creep starts when the creep rate increases dramatically causing a special type of failure known as creep rupture [7]. These three creep stages were observed in several experiments on different soil types [5,[7][8][14][15][16].
These three creep stages were attributed to continuous particles' rearrangements under the same conditions [7,14]. These structural rearrangements were explained theoretically by the continuous attempts of every one of the contacting pair of particles to reach a pair-wise equilibrium state. However, the presence of multiple contacts per particle forces each particle to keep reorienting itself targeting a global equilibrium state with all its neighbor particles. Triggered by the targeted pair-wise and global equilibrium states, the creep process in soils, its rate and dependency on time and effective stresses were conceptualized using the rate process theory [5,8]. This theory constrains the time-dependent relative deformations between pairs of particles by energy barriers separating adjacent equilibrium positions [7]. Accordingly, a relative movement between any two particles at a given equilibrium position occurs when they acquire enough energy to surmount an energy barrier and reach a new equilibrium position. The magnitude of the energy required to trigger the movement is known as the activation energy (ΔF in Fig. 2), which depends on material type, process in consideration. Furthermore, changes in the applied conditions (e.g., temperature and/or stress) can increase the energy within the system beyond the activation energy and, consequently, accelerating or slowing down the creep rate.
Energy required to cause displacement Fig. 2. Schematic of the effect of increased stress and temperature on energy barriers based on the rate process theory (adapted from [7]). Fig. 2 shows a schematic of the energy displacement for the same material at two different effective stresses and temperatures. From statistical mechanics, it is known that the average thermal energy per flow unit is kT, where k is Boltzmann constant (1.38x10-23 J K -1 ) and T is the absolute temperature (K). Even if the considered material is at rest under the same conditions, thermal vibrations occur at a frequency of kT/h, where h is Planck's constant (6.624x10-34 J s -1 ). A direct relation exists between the magnitude of the frequency of the thermal vibration and the creep rate as shown in Eq. 1 [5]. This relation indicates that when thermal vibrations increase the energy between contacting particles beyond the activation energy limit, particles move with respect to one another targeting the next equilibrium position.
in which R is the universal gas constant (8.314 J/mol. K), the parameter X is the proportionality material-dependent constant that depends on time and soil structure, λ is the distance between equilibrium positions, f is the force driving the creep, and ε is the strain rate.
Furthermore, increasing soil temperature and/or the applied stress introduce additional energy to the particulate system. Such additional energy triggers soil particles to move between equilibrium positions, accelerating the creep process. This behavior is incorporated in the rate of creep strain Eq. 1 using the sinh-based term. This term normalizes the energy due to stress changes (fλ/2) by the thermal energy (kT). The latter energy was described earlier. The former energy due to stress changes can be explained using Fig. 2. Increasing the applied stresses beyond those corresponding to the conditions for curve A will distort the heights of energy barriers. The new conditions will correspond to curve B of the distorted energy barriers. If f represents the force acting on a flow unit, then the barrier height will be reduced by an amount fλ/2 in the direction of the applied force, while it will be increased by a similar amount in the opposite direction. λ represents the distance between successive equilibrium positions as shown in Fig. 2. Thus, the distorted energy curve B will be displaced a distance δ from that of corresponding to the initial condition (i.e., curve A).
According to Eq. 1, the creep rate is expected to increase when the soil temperature and/or applied stresses increase. This suggests that clay liners used for nuclear waste repository may experience increasing creep rates due to the elevated temperatures in their vicinity. More critically, higher risks will develop if failures in the clay liners occur due to creep ruptures. Thus, we investigate creep rupture potentials of nuclear waste clay liners as triggered by elevated liner temperatures.

Overview of discrete element method (DEM)
In this study, we utilize numerical simulations based on the discrete element method (DEM). The DEM offers a unique framework to model the behavior of geomaterials at the particle scale [9]. It was originally developed to study the behavior of rocks and granular soils [10]; it was then extended to model the mechanical behavior of cohesive soils [11,[17][18]. In this method, soils are modeled through tracing the evolutions of interparticle forces between multiple moving "discrete" particles at different time steps. Typically, linear contact models are used in DEM simulations where particles' interactions are simulated using contact springs with defined normal (kn) and tangential (ks) stiffness. Thus, the normal and tangential forces between each contacting particles are calculated at each time step from the respective displacements at the contact plane. These displacements are obtained by numerically solving the equation of motion using a finite difference scheme.
While they offer unique abilities to explain basic soil responses using microscale particles interactions, DEM simulations are computationally expensive. These high computational costs come from the need to trace movements of individual particles in a soil assembly and the corresponding interparticle forces on each contacting pair. Large relative movements between particles will typically correspond to significant changes in the interparticle forces causing instabilities in numerical solvers. In order to avoid such conditions, DEM simulations use small time increments (or time steps) for the numerical solvers. Over each time step, particle interactions remain unchanged including contact forces and the particles' accelerations and velocities. More information about the DEM formulation is provided in [9].

Rate process theory in DEM
The first implementation of the rate process theory in a discrete element framework was reported in [12]. In that work, a time-dependent damping was incorporated in DEM instead of the typical artificial damping used in DEM simulations to dissipate energy. Since creep deformations correspond to particles sliding at their contacts, the utilized contact model made use of a sinhtype viscous dashpot. This dashpot acts in the tangential direction of the contact by defining the tangential sliding velocity ( s  ) based on the rate process theory as given in Eq. 2 [12].
where n1 is the number of bonds between opposing atoms at the contact, and ft and fn are the tangential and normal forces acting on the contact, respectively.

Temperature effects
The damping mechanism introduced by Eq. 2 was implemented in the 2D particle flow code (PFC 2D ) [13].
To assess the applicability of the adopted contact model in capturing the expected increase in creep rate with temperature, an assembly of two 2 mm diameter disks pressed together by a 10 N normal force was ran at different temperatures. A simple linear elastic contact model with kn = ks = 1x10 7 N/m was used. Parameters for the rate process equation (Eq. 2) were selected from the literature [12] and are listed in Table 1. The two disks were given an initial relative tangential velocity of 20 μm/sec (each particle moved at 10 20 μm/sec in opposite directions) until the maximum tangential force was developed. Then, the velocity of each of the particles was set to zero and the model was allowed to creep. The ratio between the tangential and normal forces (ft/fn) was then plotted versus time at different temperatures (Fig. 3).  As shown in Fig. 3, there is a smooth decrease of ft/fn over time due to the sliding movements associated with the creep process. Additionally, it is noticed that higher sliding rates, as indicated by the decrease of ft/fn over time, occurred at higher temperatures. This response conforms the expected temperature-dependent creep behavior as appeared in the original rate process theory (Eq. 1).

Generation of numerical specimen
To determine the potential of heating-induced creep rupture in nuclear waste clay liners, the model discussed in the previous section was applied in a DEM simulation consisting of 150 particles. While the sizes of the particles in clay liners are in the submicron range, the particle sizes considered in this study are in the 2-3 mm range. Scaling the particle sizes up for creep models is commonly adopted for adequate computational times [19][20]. Particles with larger mass result in larger time steps [21] reducing the computation time for a given creep time. The results of the models with scaled-up particle sizes, however, still represent actual responses since the selected stiffness of the contact springs are also scaled up.
The 150 particles considered in DEM simulations in this study were placed randomly in a 10x10 mm container. A linear contact model incorporating the rate process theory, with kn = ks = 1x10 7 N/m, was applied at all contacts between soil particles. To minimize boundary effects, the external walls of the container were modelled as smooth boundaries with 0.0 friction coefficient. Moreover, the contact stiffness between soil particles and the external walls were assigned the same values as the stiffness of particle-particle contacts (i.e., 1x10 7 N/m). The initial assembly was then isotropically consolidated to an all-around effective stress of σ'3 = 400 kPa. During simulation stages, the assembly was relaxed by solving for a minimum average unbalanced force ratio (R), as defined in Eq. 3, of 1x10 -5 [13].
where Up is the vector summation of the unbalanced forces acting on a soil particle, and Fp is the sum of the magnitudes of the interparticle forces acting on a soil particle.

Creep simulation
Upon stabilizing at the desired isotropic confining stress (i.e., 400 kPa), an axial deviatoric stress of 200 kPa was applied to the top boundary of the model. Thus, the specimen experienced anisotropic stress conditions corresponding to Ko of 0.67. Such condition mimics the typical stress state applied in clay liners. Additionally, the time-dependent creep damping defined earlier in Eq. 2 was incorporate into the model when the deviatoric stresses was applied. Once applied, the soil assembly was cycled for an additional 50,000 cycles to allow the model to creep. The developed model was allowed to creep at two different temperatures: 20 °C and 50 °C. The first simulation (T = 20 ℃) corresponds to the typical creep response observed in laboratory environments. While the second simulation at T = 50 ℃ mimics the response at moderate temperatures in the vicinity of nuclear waste repositories.

Results and Discussion
The volumetric strain evolutions of the two specimens during the creep stage at the two considered temperatures are shown in Fig. 4. The fundamental aspects of creep in soils are be observed in this figure. The first 1500 seconds represented the primary creep over which creep rates initially increased then decreased over time. However, the final volumetric strain magnitudes at the end of the primary creep stage in the two specimens were very close with slightly more creep in the specimen at 50 °C. Thus, it appears that the effect of temperature is less pronounced during the primary creep process.
The effects of temperature are, however, more pronounced in the secondary and tertiary creep stages. The specimen at 20 °C experienced a secondary creep stage, with almost constant creep rate, the extended until the end of the simulation (i.e., 4500 seconds). In other words, no creep rupture was predicted under the applied stress condition at 20 °C. On the other hand, the second specimen at 50 °C experienced a short ~1000 seconds secondary creep stage after which the tertiary stage started after about 2500 seconds (Fig. 4). In this stage, a sudden increase in the volumetric strain due to creep was predicted as shown in Fig. 4. Since all other model parameters were kept the same, it appears that the considered moderate temperatures caused the predicted changes in the creep rates. The results of the numerical simulations were replotted in terms of the strain rate (i.e., strain/time) as shown in Fig. 5. This figure shows that creep rate generally decreases with time during the primary and secondary creep stages. This response was found correct for the two considered temperatures (i.e., 20 and 50 °C). However, the strain rate increased during the tertiary creep stage of the specimen at 50 °C leading to a final creep rupture of the modeled specimen.
Since the main difference between the two simulations was the temperature at which creep occurred in each of the specimens, it can be stated that the main driving potential causing creep rupture was the elevated temperature. Thus, geo-structures subjected to moderate to high temperatures in the vicinity of nuclear waste repositories may experience failures due to creep ruptures. While the repositories DEM simulation at elevated temperatures reported in this study predicted creep rupture after about one hour, it should be noticed that this simulation was idealized to understand the overall behavior and does not represent the exact time to rupture for field conditions. Actual clay liners, however, may take several years to experience creep rupture under the field conditions. Thus, it is very important to constrain the time to creep rupture under expected field conditions (i.e., stresses and temperatures) in the future. Once known, clay liners will need to be designed to minimize the strain rate (Eq. 2) to prevent such creep rupture. This may be done by engineering the host material (i.e., clay liners) to increase their activation energy (ΔF) needed to initiate the creep process leading eventually to creep rupture. The authors are currently investigating the potential of engineering the clay liners in the vicinity of nuclear waste disposal facilities to mitigate the potential risks associated with creep rupture.

Conclusions
The heating-induced creep rupture of geomaterials is studied through numerical simulations based on the discrete element method (DEM). Since the excess or unbalanced energy is dissipated through sliding at the particle contact, the rate process theory was implemented in the numerical code as a damping mechanism instead of the artificial damping usually used to dissipate energy in DEM. The basics of the rate process theory was presented, after which two numerical simulations at two different temperature were carried out to study the creep behavior at higher temperature. It was shown that creep rupture can occur at temperatures comparable to those found in certain geotechnical application, posing a need to consider the possibility of encountering such behavior in actual geo-structures if not designed adequately to mitigate creep rupture phenomenon.