## 1. Introduction

Disposal of radioactive waste is currently considered as the preferred option to store and isolate this waste from biosphere over hundreds of thousands of years for waste management worldwide. In the geological disposal, the temperature in the near-field rock mass will increase due to the heat generated by the radioactive waste and in turn affect the time-dependent behavior of the host rock and consequently change the overall long-term performance of the disposal. Therefore, a good understanding of the time-dependent behavior of the host rock at elevated temperatures is of great significance to evaluate the stability of nuclear waste disposal.

Extensive effort has been made to investigate the time-dependent behavior of rocks like sedimentary rock and rock salt [1–4] and hard rocks like granite [5–10] and basalt [11]. Granite is a widely recognized potential host rock for the disposal of radioactive waste due to its low permeability and high strength [12]. Moreover, a good understanding of time-dependent deformation in granite will also assist in our understanding of crustal deformation and natural hazards since granite is a major constituent of the continental crust. The distribution of stress and permeability evolution in granite is of great importance when we analyze natural hazards and rock engineering projects. Despite the importance of granite in these high-temperature environments, very few studies of time-dependent brittle creep on granite have been made at elevated temperature [2, 8, 9, 13]. During a brittle creep experiment of rock at a constant differential stress for an extended period of time, the strain against time first decelerates before accelerating as macroscopic sample failure is approached [14]. The onset of the acceleration to failure in brittle creep experiments has been ascribed as the result of the sample reaching a microcrack density at which microcracks can interact and coalesce, sometimes referred to as the critical damage threshold [2, 11, 13, 15]. Time-dependent brittle deformation is often attributed to a mechanism of subcritical crack growth, namely stress corrosion cracking [14, 16]. Stress corrosion describes fluid-rock reactions that occur preferentially between a chemically active pore fluid and strained atomic bonds at the crack tips and therefore is sensitive to environmental factors such as stress, temperature, and pore fluid reactivity [16]. Temperature influences the crack growth rate through the Arrhenius temperature dependence of crack growth rate and because temperature affects the stress dependency of the rate of crack growth. Experimental tests have shown that speed of crack growth during experiments increases with an increase in temperature [14, 17, 18].

Due to the aforementioned importance of understanding the long-term mechanical behavior of rock for disposal of radioactive waste, here we present a two-dimensional constitutive creep model to describe the time-dependent deformation of granite at different temperatures and confining pressures. We first formulate the coupled thermomechanical time-dependent model. We then validate the model and determine the required thermo-physico-mechanical parameters with published experimental data [10, 19]. Finally, the results of brittle creep simulations at different temperatures, differential stresses, confining pressures, and sample homogeneities were presented and discussed.

## 2. Description of time-dependent constitutive model

In this chapter, a quantitative model is proposed to describe the coupled heat transfer and rock failure problems associated with rock exposed to elevated temperatures. For a model used to investigate time-dependent creep deformation at high temperature, the coupled effect of the medium deformation and heat transfer must be important. Three components including a heat transfer description, a stress description, and a failure description must be accounted for. The descriptions of heat transfer, stress, and failure in the model are, respectively, presented.

### 2.1. Material heterogeneity description

It is important to introduce material heterogeneity to reveal a collective macroscopic behavior different from that of the individual elements. The mechanical parameters such as uniaxial compressive strength and Young’s modulus of the mesoscopic material elements are assumed to be homogeneous and isotropic and are randomly assigned with a Weibull statistical distribution (Eq. (1)) [19] to reflect the material heterogeneity at a mesoscale:

in which *u* is the scale parameter of an individual element at a mesoscale and the scale parameter

in which *i*.

**Figure 1** shows two numerical standard rock specimens with a homogeneity index **Figure 1** correspond to different values of element strength and Young’s modulus. Therefore, the homogeneity index is an important controlling parameter for the macroscopic mechanical response of a numerical sample [20].

### 2.2. Heat conduction description

It is assumed that the conductive flux is saturated at a value comparable to the enthalpy per unit volume. The heat conduction equation in a medium is known as Fourier’s law:

in which ^{2}, *k _{ij}* is the thermal conductivity tensor of the medium in W/(m °C), and

It is assumed that thermal equilibrium between the phases and heat conduction as the dominant mechanism of heat transfer and the energy balance equation in an anisotropic medium can be expressed as

in which *q* denotes the heat source generated in the medium in W/m^{3}, *ρ* is the bulk density of medium in kg/m^{3}, *c* is the specific heat of the medium in J/(kg °C), and *t* is the time. The energy balance equation for an isotropic material can be reduced:

### 2.3. Stress description

Assuming the total strain for a stressed medium is made up of elastic, creep, and thermal components. The total strain can be decomposed:

in which the subscripts *e*, *c*, and *T* refer to the elastic strain, creep strain, and thermal strain, respectively. The elastic strain

in which *G* is the shear modulus,

in which *t*, i.e.,

where

For creep problems, a Norton-Bailey equation [4] known as a constitutive law of the creep strain rate was adopted to characterize time-dependent creep deformation based on the approach of the equation of state theory:

where *A*, *m*, and *n* are constants dependent on temperature. The constant *n* usually denotes stress component of greater than one, *m* is usually a fraction, *U* is the creep activation energy determined empirically as proportional to the slope of a plot of log *R* is the universal gas constant, and *T* is the absolute temperature in Kelvin.

Equation (10) can also be expressed in a strain rate form because the strain rate is of great interest for our study:

Creep flow rule under multiaxial stress conditions can be expressed in tensor form:

in which

where *S _{ij}* is the deviatoric part of

*σ*, and

_{ij}*σ*is effective stress defined as

_{e}This creep model can describe the decelerating creep commonly seen in laboratory [14], but it fails to capture the acceleration in strain rate in the approach to macroscopic sample failure. Thus, a damage evolution law for the accelerating creep of rock is incorporated at this stage.

The static stress equilibrium equation can be expressed in tensor form:

in which *σ _{ij}* is the total stress tensor in MPa and

*f*is the body force per unit volume in MPa. The continuity equation can be expressed in terms of the displacement gradient:

_{i}in which *u* is the displacement of the medium. The constitutive equation for elastic isotropic medium can be expressed as

in which *G* is shear modulus,

Thus, the final governing equation of heat transfer in a medium can be written in displacement form:

### 2.4. Damage evolution description

An elastic damage constitutive law is incorporated in the model to describe its stress-strain relationship. The elastic modulus for an isotropic elastic medium is expressed:

where *E* and *E*_{0} are the Young’s moduli of the damaged and the undamaged material, respectively, and

According to damage mechanics theory, the Young’s modulus of an element will degrade according to Eq. (19) when the stress on the element exceeds a damage threshold. The stress-strain curve of the element is considered linear elastic with a constant residual strength until the given damage threshold is reached. The elastic damage constitutive law of each element under uniaxial stress condition is shown in **Figure 2**. The maximum tensile stress criterion and modified Mohr-Coulomb criterion are chosen as the damage thresholds to determine whether any elements are damaged in tension or shear, respectively:

where *f*_{t0} are the uniaxial compressive and tensile strength, respectively, and *ϕ* is the internal friction angle of the element. The tensile strain criterion is always preferential because the uniaxial tensile strength is far below uniaxial compressive strength of rock material.

The damage variable D according to the constitutive law as shown in **Figure 2** can be described.

where *ε*_{1} is the maximum principal strain and *ε*_{3} is the minimum principal strain. *ε*_{t0} is the maximum principal strain in tension, *ε*_{c0} is the maximum principal strain in compression, and r is a constitutive coefficient with a value of 2. The element behaves elastically during loading and unloading when the major tensile or compressive strain of the element is less than the maximum principal strain in tension or compression, respectively. The element will be damaged by Eq. (22) when either of these strain thresholds is exceeded.

## 3. Model validation

### 3.1. Constant strain rate experiments

In this section, the proposed model was calibrated with experimental data [10, 18, 21] to determine appropriate input parameters. Mechanical data from uniaxial and triaxial constant displacement rate experiments at room temperature were used to obtain the physico-mechanical input parameters at the mesoscale. The geometry of the granite samples in laboratory tests was 100 mm in length and 50 mm in width. The compressive strength of air-dried Beishan granite under confining pressures of 0, 1, and 5 MPa is 165.2, 174, and 216.8 MPa, respectively. The Young’s modulus and Poisson’s ratio of air-dried Beishan granite in uniaxial compression are 43 GPa and 0.25, respectively. The randomly generated numerical samples were discretized into 20,000 square elements. The size of the modeled sample was kept the same for all numerical simulations in the present paper. A suite of unconfined and confined compressive tests on the modeled rock samples were performed at constant room temperature and constant confining pressures of 0, 1, and 5 MPa. The physico-mechanical input parameters of the individual elements at a mesoscale used in the simulations were listed in **Table 1**.

The numerical and experimental stress-strain curves for the granite simulations are plotted in **Figure 3**. It can be seen from **Figure 3** that the simulated stress-strain curves are in good agreement with the experimental stress-strain curves. It needs to be noted that the nonlinear behavior at the early beginning of the experimental stress-strain curves is a result of the closure of preexisting compliant microcracks. This nonlinearity is not replicated in the present model because the stress-strain behavior of an element is considered linear elastic until the given damage threshold is attained. We highlight that the model input parameters were the same for each of the simulations in **Figure 3**, adding confidence that the model is capable of accurately capturing the short-term mechanical behavior of Beishan granite under different confining pressures.

### 3.2. Single-step brittle creep experiments

The numerical model was validated via a successful attempt to replicate published uniaxial and triaxial experimental data for Beishan granite. We will now run a suite of conventional brittle creep simulations on Beishan granite under different constant temperatures of 23, 50, and 90°C in order to find the required thermo-physico-mechanical properties of the granite. These thermophysical parameters were determined from the experimental data, and material constants *A*, *m*, and *n* for the numerical model are listed in **Table 1**. It is noted that the same physico-mechanical input parameters listed in **Table 1** were used for the simulations. The temperatures of 23, 50, and 90°C represent room temperature, the in situ rock temperature, and the maximum temperature on the canister surface, respectively.

The numerically simulated creep curves together with the experimental creep curves are plotted in **Figure 4**. It can be seen that the simulated creep curves agree well with the experimental curves. The numerical creep curves clearly capture the phenomenology of brittle creep: the strain rate first decelerates followed by an acceleration in strain rate prior to final macroscopic failure.

**Figure 5** presents the snapshots of the damage evolution of the numerical specimens deformed at constant temperatures of 23, 50, and 90°C. **Figure 5** shows when and where damage and failure occur in the numerical specimen. The value of Young’s modulus is randomly distributed since the numerical rock samples are heterogeneous. Therefore, the elements with a low Young’s modulus act as nucleation sites for damage. As time goes on, these damaged elements grow and form localized damage zones. The localized damaged zones alter the stress field in their surrounding region, and these alterations further trigger the dynamic extension of the damage zone. Eventually, a thoroughgoing macroscopic fracture forms that signals the macroscopic failure of the sample.

### 3.3. Multistep brittle creep experiments

Chen et al. [18] also performed a series of multistep creep experiments to investigate the influence of temperature and stress on the time-dependent behavior of Beishan granite. Multistep creep tests were performed at confining pressures of 0, 1, and 5 MPa and at temperatures of 23°C (room temperature) and 90°C (the maximum temperature on the canister surface according to the current disposal conceptual design in China), respectively. The stress applied on the sample during the experiments was increased stepwise to predetermined percentages of the average peak stress (20, 40, 60, and 80%). The sample was kept at each level of stress for 1 week. We also performed numerical multistep creep simulations under the same conditions to further validate our model. It is noted that the numerical multistep creep simulations were all performed using the determined physico-mechanical and thermophysical parameters listed in **Table 1**. During the simulations, the rock samples are fixed in the vertical direction but can move freely in the horizontal direction.

The numerically simulated multistep creep curves together with the experimental multistep creep curves are plotted in **Figure 6**. It can be seen from **Figure 6** that the simulated multistep creep curves are in good agreement with the experimental curves. In detail, the model captures the influence of both confining pressure and temperature on the mechanical behavior, and the strain at the different stress steps and the time to failure are very similar between the experiment and the model (**Figure 6**).

We therefore conclude that, based on the above validations, our model can be used to investigate time-dependent creep of low-porosity granite at different temperatures. We will now use our model to further explore the influence of temperature, differential stress, confining pressure, and sample heterogeneity on brittle creep in low-porosity granite.

## 4. Numerical simulations

### 4.1. Model setup

In this section, the influence of temperature, differential stress, confining pressure, and sample heterogeneity on brittle creep of granite was investigated with the proposed model. The numerical samples as shown in **Figure 7** are the same as the samples modeled in the validation described above. The modeled samples were discretized into 20,000 square elements. We applied various axial stresses of 140, 145, 150, 155, and 160 MPa and various constant temperatures 23, 40, 50, 75, and 90°C on numerical samples with a inhomogeneity index 4, 5, and 6, respectively, to study the influence of temperature, differential stress, and sample heterogeneity on brittle creep of granite. The loading conditions are also shown in **Figure 7**. The relevant model parameters in the simulations are the same as the parameters listed in **Table 1**.

### 4.2. Effect of temperature on brittle creep

In addition to the validations above, two additional simulations of brittle creep were conducted under uniaxial compressive stress of 150 MPa at constant temperatures of 40 and 75°C. The creep curves and creep strain rate curves for various constant temperatures are shown in **Figure 8**, and the curves clearly show the accelerating-decelerating phenomenology of brittle creep observed in laboratory tests. We can see that there is a clear temperature effect on brittle creep in granite from these simulations.

It is noted that the creep strain rate strongly depends on temperature as observed in brittle creep experiments. The evolution of creep strain rate at various constant temperatures is shown in **Figure 8**. The strain rate first decreases, reaches the minimum creep strain rate, and finally increases as the sample approaches macroscopic failure. The simulations show that there are several orders of magnitude difference in the minimum creep strain rate between a lower temperature (23°C) and a higher temperature (90°C). A large increase in strain rate at a higher temperature results in a large decrease in the time to failure (**Figure 9**) as observed in laboratory experiments on granite [8].

### 4.3. Effect of differential stress on brittle creep

Here, the results of single-step brittle creep experiments under different constant-applied differential stresses were presented to study the effect of differential stress on brittle creep in granite. Uniaxial creep tests at a constant temperature of 50°C and constant axial stresses of 140, 145, 150, 155, and 160 MPa are numerically performed, and the numerically obtained creep curves for the five simulations are shown in **Figure 10**. Likewise, **Figure 10** clearly shows the accelerating-decelerating phenomenology of brittle creep seen in laboratory tests. The simulations indicate that the differential stress has a great influence on brittle creep in granite, as observed in brittle creep experiments [14]. First, the creep strain rate is higher when the differential stress is higher. A change in differential stress from 140 to 160 MPa leads to an increase in the minimum strain rate by over an order of magnitude. Similar observations have been experimentally obtained for many rock types [14]. As a result of the higher strain rate at higher differential stress, the time-to-failure is reduced as differential stress is increased (**Figure 11**).

Moreover, brittle creep tests under different constant confining pressures of 0, 2, and 10 MPa, but the same constant temperature and applied axial stress of 50°C and 150 MPa, respectively, were also performed, and the brittle creep curves were presented in **Figure 12**. The simulations capture the decelerating-accelerating phenomenology in laboratory experiments and reveal that confining pressure plays a great influence on brittle creep in granite. The creep strain rate reduces with an increase in the confining pressure (**Figure 12**): an increase in confining pressure from 0 to 10 MPa leads to a reduction in the minimum strain rate by about an order of magnitude; the results are consistent with those from brittle creep experiments at different confining pressures. **Figure 12** also shows snapshots of each of the failed samples. It can be seen that more localized shear damage occurred in rock samples at higher confining pressures.

### 4.4. Effect of material heterogeneity on brittle creep

As we know, rock is heterogeneous, and thus we use a statistical Weibull distribution to reproduce material heterogeneity in rock. A set of simulations were performed to study the effect of material heterogeneity on brittle creep in granite with different homogeneity indices but at a temperature of 50°C, confining pressure of 0 MPa, and constant-applied axial stress of 150 MPa. The mean Young’s modulus and mean UCS of the elements were kept the same (43 GPa and 350 MPa, respectively), and only the homogeneity indices were changed in these simulations. As described above, a larger homogeneity index implies that the elements within the sample are closer to the mean Young’s modulus and the mean UCS. Therefore, a sample with a larger homogeneity index indicates that there are fewer low-strength elements in rock sample, and it will become stronger and more brittle.

The simulated creep curves and the evolution of creep strain rate presented in **Figure 13** show that an increase in material homogeneity leads to a decrease in creep strain rate and an increase in time to failure. For example, the minimum creep strain rate decreases by about an order of magnitude, and the time to failure increases from about 111,600 to 358,200 s when the homogeneity index increases from 4 to 6. On the contrary, a decrease in material homogeneity results in an increase in the creep strain rate and a decrease in the time to failure.

## 5. Concluding remarks

A numerical time-dependent thermomechanical model was proposed to simulate brittle creep in granite under different loading conditions and different temperatures. The mechanical parameters such as uniaxial compressive strength and Young’s modulus of the mesoscopic material elements assumed to be homogeneous and isotropic are randomly assigned with a Weibull statistic distribution to reflect the material heterogeneity at a mesoscale. It is noted that the model can well capture the cooperative interaction between microcracks in the transition from distributed to localized damage. The model is validated against published experimental data, and then conventional brittle creep experiments at various constant temperatures, applied differential stresses, confining pressures, and sample homogeneities were simulated. The simulations accurately capture the short- and long-term mechanical behavior of the experimental brittle creep tests using unique thermo-physico-mechanical properties. The simulations further show that an increase in temperature and differential stress leads to an increase in the creep strain rate and therefore a decrease in time to failure and an increase in confining pressure and sample homogeneity leads to a decrease in creep strain rate and an increase in time to failure. Thus, the model proposed in the present paper will help the management and optimization of rock engineering projects in granite.