A Coupled Thermal–Hydrological–Mechanical Damage Model and Its Numerical Simulations of Damage Evolution in APSE

This paper proposes a coupled thermal–hydrological–mechanical damage (THMD) model for the failure process of rock, in which coupling effects such as thermally induced rock deformation, water flow-induced thermal convection, and rock deformation-induced water flow are considered. The damage is considered to be the key factor that controls the THM coupling process and the heterogeneity of rock is characterized by the Weibull distribution. Next, numerical simulations on excavation-induced damage zones in Äspö pillar stability experiments (APSE) are carried out and the impact of in situ stress conditions on damage zone distribution is analysed. Then, further numerical simulations of damage evolution at the heating stage in APSE are carried out. The impacts of in situ stress state, swelling pressure and water pressure on damage evolution at the heating stage are simulated and analysed, respectively. The simulation results indicate that (1) the v-shaped notch at the sidewall of the pillar is predominantly controlled by the in situ stress trends and magnitude; (2) at the heating stage, the existence of confining pressure can suppress the occurrence of damage, including shear damage and tensile damage; and (3) the presence of water flow and water pressure can promote the occurrence of damage, especially shear damage.


Introduction
In recent years, the coupled thermal-hydrological-mechanical (THM) behaviour in porous media has been a subject of great interest in many engineering disciplines, such as geothermal energy extraction, recovery of hydrocarbons, safety assessment of waste repositories for radioactive waste, and geological sequestration of CO 2 [1][2][3][4]. Based on Biot's theory of the isothermal consolidation of elastic porous media, the phenomenological approach of poroelasticity [5,6], and the mixture theory as described by Morland [7] and Bowen [8], many coupled THM models and computer codes have been developed [3,[9][10][11]. In particular, the framework of international projects like DECOVALEX (DEvelopment of COupled models and their VALidation against Experiments) [12][13][14][15], and the tunnel sealing experiment in the Underground Research Laboratory, Canada [16] has dramatically promoted the study of the THM coupling process [17].
In addition, the creation of an excavation-damaged zone (EDZ) is expected around all man-made openings in geologic formations. Macro-and meso-fracturing and, in general, a rearrangement of rock structures, occur in this zone, resulting in drastic changes in the mechanical parameters [13].
where G is the shear modulus and v is Poisson's ratio. α (≤ 1) is the Biot's coefficient, depending on the compressibility of the constituents, and can be defined as α = 1 − K K s , where K s is the effective bulk modulus of the solid constituent (Pa) and K = 2G(1+v) 3(1−2v) is the bulk modulus of the porous medium (Pa). α T is the volumetric expansion coefficient of the bulk medium (K −1 ).
Then a modified Navier equation, in terms of displacement u i (m), water pressures p (Pa), and temperature change T (K), can be expressed as where σ ij is the components of the net body force (N·m −3 ) in the i-direction σ ij .

Water Flow Equation
According to Zhou et al. [24], if the water flow follows Darcy's law: The mass conservation equation for the water flow process in porous rock can be expressed as where where ε v is the volumetric strain, µ w is the dynamic viscosity of water (N·s·m −2 ), k is the intrinsic permeability of rock (m 2 ), ρ w is the water density (kg·m −3 ), z is the vertical coordinate, g is the gravitational acceleration (m·s −2 ), φ is the porosity of rock, a w and a s are the thermal expansion coefficient of water and rock grain under constant pore pressure and stress (1/K), respectively, and K w denotes the bulk modulus of water.

Energy Conservation Equation
Based on Fourier's law, the constitutive relation for heat diffusion can be expressed as where q T is the thermal flux (J·s −1 ·m −2 ), C w is the water-specific heat constant at a constant volume (J·kg −1 ·K −1 ), and q w is the Darcy velocity as expressed in Equation (3). λ M , λ s , and λ w are the thermal conductivities of rock medium, rock grain, and water components (J·s −1 ·m −1 ·K −1 ), respectively. The first and second terms on the right-hand side of Equation (6) are the heat conduction by fluid-solid mixture and the heat convection by water flow, respectively. Assuming thermal equilibrium between the fluid and solid phases, the energy conservation equation can be expressed as [3,6,24,25]: where (ρC) M is the specific heat capacity of the water-filled solid medium and ρ s is the mass density of the rock matrix (kg·m −3 ). C w and C s are the water and solid specific heat constants at a constant volume (J·kg −1 ·K −1 ), respectively. T ar is the absolute reference temperature in the stress-free state (K). The first term on the left-hand side of Equation (8) represents the internal heat energy change rate per unit of volume due to temperature change. The second and third terms represent the heat sink due to the thermal dilatation of water and rock, respectively [24]. Assuming constant specific heats (C w and C s ) and thermal conductivities (λ s and λ w ), substituting Equations (6), (7), and (9) into Equation (8) yields

Damage Evolution Equation
As illustrated in Figure 1, when the stress state of rock under coupled THM conditions satisfies the maximum tensile stress criterion or the Mohr-Coulomb criterion, tensile damage or shear damage will occur, respectively. The criteria can be expressed as [26][27][28]: where σ 1 and σ 3 are the first and third principal stresses (Pa), respectively; f t0 and f c0 are the uniaxial tensile and compressive strength (Pa), respectively; and ϕ is the internal frictional angle ( • ).

Effect of Damage on THM Parameters
The impact of rock damage on mechanical parameters was also considered, based on the damage theory. The elastic modulus and strength of an element degrade monotonically as damage evolves, and can be expressed as [28,29]: According to Figure 1, the damage variable D can be calculated as [26][27][28]: where ε 1 and ε 3 are the major and minor principal strain; ε t0 and ε c0 are the maximum tensile principal strain and the maximum compressive principal strain when tensile and shear damage occur, respectively; and n is a constitutive coefficient, specified as 2.0. More details can be found in previous publications [26][27][28].

Effect of Damage on THM Parameters
The impact of rock damage on mechanical parameters was also considered, based on the damage theory. The elastic modulus and strength of an element degrade monotonically as damage evolves, and can be expressed as [28,29]: where E and E 0 are the elastic modulus (Pa) of the damaged and undamaged element, respectively; and σ c and σ c0 are the compressive strength (Pa) of the damaged and undamaged element, respectively. The permeability is correlated to the damage variable according to the following exponential function [29]: where k 0 is the permeability of the undamaged element and α k is a damage-permeability effect coefficient to indicate the effect of damage on the permeability. Similarly, the effect of damage on thermal conductivity is given as: where λ s0 is the thermal conductivity of the undamaged element and α λ is a coefficient to reflect the effect of damage on the thermal conductivity. Equations (2), (4), (10), and (12) represent the fully-coupled thermal-hydrological-mechanical damage response of rock. The coupled non-linear equation is then implemented and solved directly using COMSOL Multiphysics [30], a powerful partial differential equation-based multiphysics modeling environment, together with programming based on Matlab. More details about the implementation and validation of this kind of non-linear equation with COMSOL Multiphysics have been presented in previous publications [25,28].

APSE Background
As mentioned previously, placing radioactive waste in suitable underground rock formations was one of the early options proposed as a way to remove waste from the biosphere [31]. The Swedish Nuclear Fuel and Waste Management Company (SKB) is undertaking site characterization for spent nuclear fuel. In the KBS-3 concept, the spent fuel is encapsulated in canisters and deposited in vertical 8 m deep, 1.75 m diameter deposition holes excavated in the floor of horizontal deposition tunnels at 400-700 m depths below the surface [32]. At such depths, the excavations induce stress concentrations around the holes and may induce damage in the surrounding rocks. In addition, the decay of used nuclear fuel induces temperature increase, thermal expansion, and thermal stress in surrounding rocks. To address the issue, the Äspö pillar stability experiment (APSE) was carried out in the Äspö Hard Rock Laboratory (HRL) to investigate the EDZ development and failure process in rock mass subjected to coupled excavation-induced and thermally induced stresses [20].
A drift with a width of 5 m and height of 7.5 m was excavated first for the experiment. Then, two 1.75 m diameter bore holes were excavated close to each other to create the 1 m wide pillar ( Figure 2). A swelling pressure of 0.7 MPa was applied to the confining hole to simulate the backfilling effect of the swelling clay bentonite mixture. Four vertical boreholes were then drilled and four electrical heaters were installed to heat the pillar area ( Figure 2). for spent nuclear fuel. In the KBS-3 concept, the spent fuel is encapsulated in canisters and deposited in vertical 8 m deep, 1.75 m diameter deposition holes excavated in the floor of horizontal deposition tunnels at 400-700 m depths below the surface [32]. At such depths, the excavations induce stress concentrations around the holes and may induce damage in the surrounding rocks. In addition, the decay of used nuclear fuel induces temperature increase, thermal expansion, and thermal stress in surrounding rocks. To address the issue, the Äspö pillar stability experiment (APSE) was carried out in the Äspö Hard Rock Laboratory (HRL) to investigate the EDZ development and failure process in rock mass subjected to coupled excavation-induced and thermally induced stresses [20].
A drift with a width of 5 m and height of 7.5 m was excavated first for the experiment. Then, two 1.75 m diameter bore holes were excavated close to each other to create the 1 m wide pillar ( Figure 2). A swelling pressure of 0.7 MPa was applied to the confining hole to simulate the backfilling effect of the swelling clay bentonite mixture. Four vertical boreholes were then drilled and four electrical heaters were installed to heat the pillar area ( Figure 2). : Heating hole

Determination of Meso-Mechanical Parameters
The mean value of elastic modulus and uniaxial compressive strength (UCS) of intact rock sample in APSE is 76 GPa and 211 MPa, respectively [20]. Because the impact of existing fractures on the pillar damage process is not considered in the numerical model in this paper, reduced mechanical parameters are used instead. According to Barton [33], the reduced elastic modulus and uniaxial compressive strength of rock mass range from 58.5 to 65.6 GPa and 46 to 58 MPa, respectively. Hence, the reduced elastic modulus of 62 GPa and UCS of 52 MPa are used in the present study.
According to the Finite Element Method, the numerical sample will be divided into many elements first to numerically solve the THMD model proposed in Section 2. These elements are called mesoscopic elements here and the mechanical parameters of these elements are called meso-mechanical parameters. The mechanical parameters such as Elastic Modulus and UCS obtained from uniaxial compressive test are called macro-mechanical parameters. In addition, a Weibull distribution, as defined in the following probability density function, is used here to characterize the heterogeneity of rock [29]: where c is the mechanical parameter (such as elastic modulus, UCS, and thermal conductivity) of the mesoscopic elements in the numerical specimen, 0 c is the scale parameter related to the average of the material parameters, and the shape parameter m reflects the degree of material

Determination of Meso-Mechanical Parameters
The mean value of elastic modulus and uniaxial compressive strength (UCS) of intact rock sample in APSE is 76 GPa and 211 MPa, respectively [20]. Because the impact of existing fractures on the pillar damage process is not considered in the numerical model in this paper, reduced mechanical parameters are used instead. According to Barton [33], the reduced elastic modulus and uniaxial compressive strength of rock mass range from 58.5 to 65.6 GPa and 46 to 58 MPa, respectively. Hence, the reduced elastic modulus of 62 GPa and UCS of 52 MPa are used in the present study.
According to the Finite Element Method, the numerical sample will be divided into many elements first to numerically solve the THMD model proposed in Section 2. These elements are called mesoscopic elements here and the mechanical parameters of these elements are called meso-mechanical parameters. The mechanical parameters such as Elastic Modulus and UCS obtained from uniaxial compressive test are called macro-mechanical parameters. In addition, a Weibull distribution, as defined in the following probability density function, is used here to characterize the heterogeneity of rock [29]: where c is the mechanical parameter (such as elastic modulus, UCS, and thermal conductivity) of the mesoscopic elements in the numerical specimen, c 0 is the scale parameter related to the average of the material parameters, and the shape parameter m reflects the degree of material homogeneity and is defined as a homogeneity index. For higher values of m, the mechanical parameters of more elements are concentrated closer to c 0 [28]. Hence, a series of uniaxial compressive tests on numerical samples with different meso-mechanical parameters was conducted first. When the macro-mechanical parameters obtained from the numerical sample were close to the reduced parameters from the laboratory experiments, the corresponding meso-mechanical parameters were used in later APSE simulations. Figure 3 shows the damage and failure process of a numerical sample (50 mm wide and 100 mm high) under uniaxial displacement of 0.01 mm/step at the upper boundary. Here, a homogeneity index m of 5, a mean elastic modulus of 68 GPa, and a UCS of 119 MPa are used for the mesoscopic elements of the numerical sample. Figure 3a shows the damage variable evolution, in which the tensile damage is represented as negative numbers (−1 ≤ D < 0), while the shear damage is represented as positive numbers (0 < D ≤ 1), in order to clearly distinguish the two kinds of damage modes. There is no damage in the intact rock sample before compression (Figure 3a, Step 0). As the axial displacement increases, first shear damage initiates, which is randomly distributed in the sample (Figure 3a, Step 80). Then, tensile damage occurs around the shear damage element (Figure 3a, Step 86) and coalescences gradually (Figure 3a, Step 86), leading to the formation of the main macro-fracture (Figure 3a, Step 88). The main macro-fracture located around the middle part of the sample then extends quickly (Figure 3a, Steps 92-94), leading to the final unstable failure of the sample (Figure 3a, Step 100).
According to Figure 1 and Equation (12), each mesoscopic element is damaged gradually under external loading. This means that the damage variable of each element increases monotonically, while damage modes (shear or tensile) can change with stress state. This is why the initial shear damage elements can finally develop into macro tensile fracture. Figure 3b is the elastic modulus evolution, in which the initial elastic modulus distribution is specified according to a Weibull distribution with a homogeneity index m of 5 ( Figure 3b, Step 0). The elastic modulus decreases gradually with damage developing in the mesoscopic elements, shown as black areas in the sample (Figure 3b, Steps 80-100). The reduction of elastic modulus causes stress redistribution around the damaged elements, leading to the extension and coalescence of fractures and the final failure of the sample. homogeneity and is defined as a homogeneity index. For higher values of m , the mechanical parameters of more elements are concentrated closer to 0 c [28].
Hence, a series of uniaxial compressive tests on numerical samples with different meso-mechanical parameters was conducted first. When the macro-mechanical parameters obtained from the numerical sample were close to the reduced parameters from the laboratory experiments, the corresponding meso-mechanical parameters were used in later APSE simulations. Figure 3 shows the damage and failure process of a numerical sample (50 mm wide and 100 mm high) under uniaxial displacement of 0.01 mm/step at the upper boundary. Here, a homogeneity index m of 5, a mean elastic modulus of 68 GPa, and a UCS of 119 MPa are used for the mesoscopic elements of the numerical sample. Figure 3a shows the damage variable evolution, in which the tensile damage is represented as negative numbers (−1 ≤ D < 0), while the shear damage is represented as positive numbers (0 < D ≤ 1), in order to clearly distinguish the two kinds of damage modes. There is no damage in the intact rock sample before compression (Figure 3a, Step 0). As the axial displacement increases, first shear damage initiates, which is randomly distributed in the sample (Figure 3a, Step 80). Then, tensile damage occurs around the shear damage element (Figure 3a, Step 86) and coalescences gradually (Figure 3a, Step 86), leading to the formation of the main macro-fracture ( Figure 3a, Step 88). The main macro-fracture located around the middle part of the sample then extends quickly (Figure 3a, Steps 92-94), leading to the final unstable failure of the sample (Figure 3a, Step 100).
According to Figure 1 and Equation (12), each mesoscopic element is damaged gradually under external loading. This means that the damage variable of each element increases monotonically, while damage modes (shear or tensile) can change with stress state. This is why the initial shear damage elements can finally develop into macro tensile fracture. Figure 3b is the elastic modulus evolution, in which the initial elastic modulus distribution is specified according to a Weibull distribution with a homogeneity index m of 5 ( Figure 3b, Step 0). The elastic modulus decreases gradually with damage developing in the mesoscopic elements, shown as black areas in the sample (Figure 3b, Steps 80-100). The reduction of elastic modulus causes stress redistribution around the damaged elements, leading to the extension and coalescence of fractures and the final failure of the sample.
Step 0 Step 80 Step 86 Step 88 Step 92 Step 94 Step 100 (a) Step 0 Step 80 Step 86 Step 88 Step 92 Step 94 Step 100 (b)  During numerical simulation of uniaxial compressive test, axial strain and stress at the upper boundary at each step are calculated to gain the complete stress-strain curve, from which the macro elastic modulus and UCS can be obtained. Finally, the macro elastic modulus and UCS of 62.8 GPa and 52.6 MPa were obtained ( Figure 4) based on the uniaxial compressive test above, which are very close to the reduced parameters of 62 GPa and 52 MPa from the laboratory experiments. Hence, the mean elastic modulus of 68 GPa and a UCS of 119 MPa for the mesoscopic elements are used in later simulations.
Based on Andersson and Martin [20] and the simulations above, the mechanical parameters used in later simulations are listed in Table 1.
During numerical simulation of uniaxial compressive test, axial strain and stress at the upper boundary at each step are calculated to gain the complete stress-strain curve, from which the macro elastic modulus and UCS can be obtained. Finally, the macro elastic modulus and UCS of 62.8 GPa and 52.6 MPa were obtained ( Figure 4) based on the uniaxial compressive test above, which are very close to the reduced parameters of 62 GPa and 52 MPa from the laboratory experiments. Hence, the mean elastic modulus of 68 GPa and a UCS of 119 MPa for the mesoscopic elements are used in later simulations. Based on Andersson and Martin [20] and the simulations above, the mechanical parameters used in later simulations are listed in Table 1.

Determination of In Situ Stress and Boundary Conditions
Based on hydraulic fracturing, triaxial over-coring, and back analysis, Andersson and Martin [20] give the in situ stress state for the APSE site as shown in Table 2.

Determination of In Situ Stress and Boundary Conditions
Based on hydraulic fracturing, triaxial over-coring, and back analysis, Andersson and Martin [20] give the in situ stress state for the APSE site as shown in Table 2. In the horizontal direction, the trends of major principal stress σ 1 and minor principal stress σ 3 in the Äspö 96 coordinate system are 310 • and 220 • , respectively. The APSE tunnel trend is 46 • in the Äspö 96 coordinate system. As a result, the tunnel axis was located approximately 6 • off the trend of σ 3 , i.e., the tunnel axis was not at right angles to the trend of σ 1 [34]. Hence, a numerical model as shown in Figure 5 is established in this study. One of the benefits of such a model is that only normal stress will exist at the outer boundary of the model and the tangential stress at outer boundary is zero.
Being limited by the number of mesoscopic elements and computing power for the rock damage process, the actual three-dimensional experiment was simplified to a two-dimensional model. Therefore, pillar stability modeling was performed in one plane of interest, 3.5 m below the tunnel floor. This means that only the major principal stress σ 1 and the minor principal stress σ 3 are considered. The size of the numerical model is 7.5 m × 7.5 m. The dotted red line in Figure 5a represents the tunnel trend, along which two 1.75 m diameter bore holes are excavated close to each other to create the 1 m wide pillar. The angle between the APSE tunnel trend and the minor principal stress σ 3 is represented by β (equalling 6 • for the stress trend in Table 2).

Numerical Model for Excavation Stage
As shown in Figure 5a, at the excavation stage, the vertical boundary stress σ by (corresponding to σ 3 ) and the horizontal stress σ bx (corresponding to σ 1 ) are applied in monotonically increasing mode with an increment ∆σ by of 0.2 MPa and ∆σ bx of 0.6 MPa per step to observe the stable crack propagation process. A confining pressure of 0.7 MPa is applied around the confining hole. An initial water pressure p 0 of 4.5 MPa is applied at the outer boundary and an atmospheric pressure of 0.1 MPa is applied around the holes. Figure 5b is the elastic modulus distribution with the meso-mechanical parameters listed in Table 1.
The excavation of the APSE tunnel causes the in situ stress to redistribute and concentrate within the surrounding rock. The stress state in the pillar decreases with pillar depth. As a result, the precise stress state around the experimental area is difficult to implement. Figure 6a shows the excavation-induced tangential stress measured at different depths [35]. Figure 6b shows our simulation results, from which it can be seen that when σ by increases gradually from 1 MPa to 20 MPa, σ bx increases gradually from 3 MPa to 60 MPa, and the tangential stress at point C (shown in Figure 7) increases from 10 to 190 MPa. Hence, different σ by and σ bx can be applied to the outer boundary of the numerical model to simulate the stress states at corresponding depths.
shown in Figure 5 is established in this study. One of the benefits of such a model is that only normal stress will exist at the outer boundary of the model and the tangential stress at outer boundary is zero.
Being limited by the number of mesoscopic elements and computing power for the rock damage process, the actual three-dimensional experiment was simplified to a two-dimensional model. Therefore, pillar stability modeling was performed in one plane of interest, 3.5 m below the tunnel floor. This means that only the major principal stress 1  and the minor principal stress 3  are considered. The size of the numerical model is 7.5 m × 7.5 m. The dotted red line in Figure 5a represents the tunnel trend, along which two 1.75 m diameter bore holes are excavated close to each other to create the 1 m wide pillar. The angle between the APSE tunnel trend and the minor principal stress σ3 is represented by  (equalling 6° for the stress trend in Table 2).

Numerical Model for Excavation Stage
As shown in Figure 5a, at the excavation stage, the vertical boundary stress  Figure 5b is the elastic modulus distribution with the meso-mechanical parameters listed in Table 1.
The excavation of the APSE tunnel causes the in situ stress to redistribute and concentrate within the surrounding rock. The stress state in the pillar decreases with pillar depth. As a result, the precise stress state around the experimental area is difficult to implement. Figure 6a shows the excavation-induced tangential stress measured at different depths. Figure 6b shows our simulation results, from which it can be seen that when by  increases gradually from 1 MPa to 20 MPa, bx  increases gradually from 3 MPa to 60 MPa, and the tangential stress at point C (shown in Figure 7) increases from 10 to 190 MPa. Hence, different by  and bx  can be applied to the outer boundary of the numerical model to simulate the stress states at corresponding depths.

Numerical Model for Heating Stage
After the excavation stage, four heating holes are drilled further to simulate the heating stage ( Figure 7). Heater powers used in the experiment and in the numerical simulation are listed in Table 3. The temperature at the outer boundary is set to the initial rock formation temperature 0 T of 15 °C. A natural air convection boundary is set on the surface of the deposit holes. The heater powers used in the simulation are based on the report and a related study [35,36]. In addition, an initial water pressure 0 p of 4.5 MPa is applied at the outer boundary and an atmospheric pressure of 0.1 MPa is applied around the holes. Outer boundary stresses by  of 12, 13, and 14 MPa are applied to simulate stress conditions at different depths (corresponding to 3.5, 2.5 and 2.0 m, respectively) in the pillar. Table 3. Heater powers used in the experiment and in the model [35,36].

Numerical Model for Heating Stage
After the excavation stage, four heating holes are drilled further to simulate the heating stage ( Figure 7). Heater powers used in the experiment and in the numerical simulation are listed in Table 3. The temperature at the outer boundary is set to the initial rock formation temperature T 0 of 15 • C. A natural air convection boundary is set on the surface of the deposit holes. The heater powers used in the simulation are based on the report and a related study [35,36]. In addition, an initial water pressure p 0 of 4.5 MPa is applied at the outer boundary and an atmospheric pressure of 0.1 MPa is applied around the holes. Outer boundary stresses σ by of 12, 13, and 14 MPa are applied to simulate stress conditions at different depths (corresponding to 3.5, 2.5 and 2.0 m, respectively) in the pillar. Table 3. Heater powers used in the experiment and in the model [35,36].  0  200  200  200  180  18  170  160  200  160  38  354  348  400  320  42  354  348  400  320  46  263  263  400  300  52  263  263  400  280  59  263  263  400  270  66  0  0  0  0 The calculated temperatures at points A, B, C, and D are shown together with the corresponding measurements in Figure 8. The figure shows that a similar temperature distribution is obtained, so that the pillar damage evolution at the heating stage can be simulated further.

Days from Heater
It should be noted that a precise description of the temperature field and a detailed discussion on thermal conduction and the thermal convection process based on energy conservation Equation (10) are not provided here, due to the complexity of the APSE site. For example, the monitoring system showed that there was water inflow in three of the cored boreholes and in both of the deposition holes. The water inflow rate in one of the heating holes was up to 5.36 L/min, leading to a dramatic cooling effect in the pillar. Hence, a packer was placed above the fracture and the water was led through the packer pipe out of the hole to weaken the cooling effect of heat convection induced by water flow [35]. All these options make it difficult to precisely recreate the temperature field in simulations. For this reason, obtaining the temperature variation close to the measurement (as shown in Figure 8) and then analysing the pillar damage evolution in such temperature fields are the main purposes of this study. deposition holes. The water inflow rate in one of the heating holes was up to 5.36 L/min, leading to a dramatic cooling effect in the pillar. Hence, a packer was placed above the fracture and the water was led through the packer pipe out of the hole to weaken the cooling effect of heat convection induced by water flow [35]. All these options make it difficult to precisely recreate the temperature field in simulations. For this reason, obtaining the temperature variation close to the measurement (as shown in Figure 8) and then analysing the pillar damage evolution in such temperature fields are the main purposes of this study.  Figure 9a shows the damage zone evolution at the excavation stage with β of −6°, in which the tensile damage is represented as negative numbers (−1 ≤ D < 0), while the shear damage is represented as positive numbers (0 < D ≤ 1). In the numerical model, the "damage zone" represents the mesoscopic element of which the absolute value of damage variable D is greater than 1. In physics, the "damage zone" represents the mineral grain or small rock block in which shear damage or tensile damage occurs. At Step 60, shear damage first emerges near points C and D of the shear stress concentration area. With the increase in the external load, this shear damage zone grows gradually (at Step 65) and propagates much deeper into the surrounding rock (at Step 70), resulting in a v-shaped notch in the sidewall (at Step 75). Meanwhile, tensile damage also occurs around the shear damage elements in the v-shaped notch area.

Damage Zone Evolution at Excavation Stage
It is noted that the v-shaped notch shows unsymmetrical distribution along the tunnel axis as a result of the deflection angle between the tunnel trend and the in situ major principal stress trend. The shear damage around the open hole is mostly located on the left side of the tunnel axis, while the shear damage around the confining hole is mostly located on the right side of the tunnel axis. Figure 10 shows the location of all the acoustic events recorded within a vertical distance of ±0.1 m from each section depth at both the excavation stage and the heating stage [34]. The acoustic events around the open hole are mostly located on the right side of the tunnel axis, which is obviously inconsistent with the simulation results shown in Figure 9a. It can be speculated that the  Figure 9a shows the damage zone evolution at the excavation stage with β of −6 • , in which the tensile damage is represented as negative numbers (−1 ≤ D < 0), while the shear damage is represented as positive numbers (0 < D ≤ 1). In the numerical model, the "damage zone" represents the mesoscopic element of which the absolute value of damage variable D is greater than 1. In physics, the "damage zone" represents the mineral grain or small rock block in which shear damage or tensile damage occurs. At Step 60, shear damage first emerges near points C and D of the shear stress concentration area. With the increase in the external load, this shear damage zone grows gradually (at Step 65) and propagates much deeper into the surrounding rock (at Step 70), resulting in a v-shaped notch in the sidewall (at Step 75). Meanwhile, tensile damage also occurs around the shear damage elements in the v-shaped notch area.

Damage Zone Evolution at Excavation Stage
It is noted that the v-shaped notch shows unsymmetrical distribution along the tunnel axis as a result of the deflection angle between the tunnel trend and the in situ major principal stress trend. The shear damage around the open hole is mostly located on the left side of the tunnel axis, while the shear damage around the confining hole is mostly located on the right side of the tunnel axis. Figure 10 shows the location of all the acoustic events recorded within a vertical distance of ±0.1 m from each section depth at both the excavation stage and the heating stage [34]. The acoustic events around the open hole are mostly located on the right side of the tunnel axis, which is obviously inconsistent with the simulation results shown in Figure 9a. It can be speculated that the trend of in situ stress (σ 1 , σ 3 ) set in Figure 9a, which is determined according to in situ measurement, is not consistent with the actual in situ stress trend around the deposit holes. The excavation of the drift tunnel has influenced the stress state around the holes, causing not only stress concentration and redistribution but also stress trend deflection.
Hence, further numerical simulations on damage zone evolution under adjusted stress conditions were carried out, in which the trends of in situ stress (σ 1 , σ 3 ) were 316, 226; 320, 230; and 326, 236, respectively (as shown in Figure 9b-d).
For the σ 1 trend of 316 • , the shear damage zone around the open hole is distributed uniformly on both sides of the tunnel axis. For the σ 1 trend of 320 • , the shear damage around the open hole is mostly located on the right side of the tunnel axis, which is consistent with the location of all the acoustic events. For the σ 1 trend of 326 • , the shear damage around the open hole is also located on the right side of the tunnel axis, but the damage area is relatively small and the v-shaped notch is not formed.
A comparison of Figures 9 and 10 indicates that the damage zone distribution with the σ 1 trend of 320 • shows relatively good agreement with the location of AE events. This conclusion is also consistent with the research reported by Andersson et al. [35]. Hence, the in situ stress (σ 1 , σ 3 ) trend of 320, 230 is used in the later simulation of damage evolution at the heating stage.
Materials 2016, 9,841 12 of 20 excavation of the drift tunnel has influenced the stress state around the holes, causing not only stress concentration and redistribution but also stress trend deflection. Hence, further numerical simulations on damage zone evolution under adjusted stress conditions were carried out, in which the trends of in situ stress ( 1  , 3  ) were 316, 226; 320, 230; and 326, 236, respectively (as shown in Figure 9b-d).
Step 60 Step 65 Step 70 Step 75 (a) Step 60 Step 65 Step 70 Step 75 (b) Step 60 Step 65 Step 70 Step 75 (c) Step 60 Step 65 Step 70 Step 75 (d)  located on the right side of the tunnel axis, but the damage area is relatively small and the v-shaped notch is not formed. A comparison of Figure 9 and Figure 10 indicates that the damage zone distribution with the 1  trend of 320° shows relatively good agreement with the location of AE events. This conclusion is also consistent with the research reported by Andersson et al. [34]. Hence, the in situ stress ( 1  , 3  ) trend of 320, 230 is used in the later simulation of damage evolution at the heating stage.

Damage Zone Evolution at Heating Stage
Based on the damage zone distribution at the excavation stage as shown in Figure 9c, the damage zone evolution at the heating stage was simulated. The trends of in situ stress 1  and 3  were 320° and 230°, respectively.
In Figure 11a, the initial in situ stress values of 1  and 3  are 36 MPa and 12 MPa, respectively, corresponding to Step 60 in Figure 9c. Only a few discrete damage points emerge in the area between deposit holes and heating holes, and the shear damage zone near points C and D

Damage Zone Evolution at Heating Stage
Based on the damage zone distribution at the excavation stage as shown in Figure 9c, the damage zone evolution at the heating stage was simulated. The trends of in situ stress σ 1 and σ 3 were 320 • and 230 • , respectively.
In Figure 11a, the initial in situ stress values of σ 1 and σ 3 are 36 MPa and 12 MPa, respectively, corresponding to Step 60 in Figure 9c. Only a few discrete damage points emerge in the area between deposit holes and heating holes, and the shear damage zone near points C and D remain unchanged. In Figure 11b, the initial in situ stress values of σ 1 and σ 3 are 39 MPa and 13 MPa, respectively, corresponding to Step 65 in Figure 9c. More discrete damage points, mainly shear damage, emerge in the area between deposit holes and heating holes, and the shear damage zone near points C and D has barely changed. In Figure 11c, the initial in situ stress values of σ 1 and σ 3 are 42 MPa and 14 MPa, respectively, corresponding to Step 70 in Figure 9c. Much more shear damage occurs in the area between deposit holes and heating holes, and the v-shaped notch also shows a little growth.

Damage Zone Evolution at Heating Stage
Based on the damage zone distribution at the excavation stage as shown in Figure 9c, the damage zone evolution at the heating stage was simulated. The trends of in situ stress 1  and 3  were 320° and 230°, respectively.
In Figure 11a, the initial in situ stress values of 1  and 3   are 36 MPa and 12 MPa, respectively, corresponding to Step 60 in Figure 9c. Only a few discrete damage points emerge in the area between deposit holes and heating holes, and the shear damage zone near points C and D remain unchanged. In Figure 11b, the initial in situ stress values of 1  and 3  are 39 MPa and 13 MPa, respectively, corresponding to Step 65 in Figure 9c. More discrete damage points, mainly shear damage, emerge in the area between deposit holes and heating holes, and the shear damage zone near points C and D has barely changed. In Figure 11c  respectively. Hence, the heating-induced damage zone areas are 0.12, 0.14, and 0.25 m 2 , which are larger than the excavation-induced damage zone areas, respectively. In addition, the heating-induced damage zone area increases as the initial in situ stress increases, corresponding to  Figure 12 shows the damage zone area variation with time at the heating stage. For initial in situ stresses σ 1 of 36, 39, and 42 MPa, the excavation-induced damage zone areas are 0.06, 0.12, and 0.18 m 2 , respectively. At day 60, the damage zone area increases to 0.18, 0.26, and 0.43 m 2 , respectively. Hence, the heating-induced damage zone areas are 0.12, 0.14, and 0.25 m 2 , which are larger than the excavation-induced damage zone areas, respectively. In addition, the heating-induced damage zone area increases as the initial in situ stress increases, corresponding to a shallower position in the pillar.  Figure 12 shows the damage zone area variation with time at the heating stage. For initial in situ stresses 1  of 36, 39, and 42 MPa, the excavation-induced damage zone areas are 0.06, 0.12, and 0.18 m 2 , respectively. At day 60, the damage zone area increases to 0.18, 0.26, and 0.43 m 2 , respectively. Hence, the heating-induced damage zone areas are 0.12, 0.14, and 0.25 m 2 , which are larger than the excavation-induced damage zone areas, respectively. In addition, the heating-induced damage zone area increases as the initial in situ stress increases, corresponding to a shallower position in the pillar.  Figure 13 shows the damage zone areas occurring each day under the three different in situ stress conditions (a)-(c), the summation of the daily damage zone area (d), and accumulated acoustic emission (AE) events each day [34]. The figure shows that the damage zone occurs mainly during days 1-20 and 40-60, which is consistent with the increasing temperature. The total damage zone area each day in the simulation is also consistent with the AE events.   Figure 13 shows the damage zone areas occurring each day under the three different in situ stress conditions (a)-(c), the summation of the daily damage zone area (d), and accumulated acoustic emission (AE) events each day [35]. The figure shows that the damage zone occurs mainly during days 1-20 and 40-60, which is consistent with the increasing temperature. The total damage zone area each day in the simulation is also consistent with the AE events.

Effect of Confining Pressure on Damage Zone Evolution
In underground repositories for radioactive waste, clay bentonite mixtures are usually used to backfill the deposit holes to seal off the canisters. This kind of bentonite mixture swells considerably when exposed to water, resulting in swelling pressure being applied to the rock walls. As shown in Figure 7, a confining pressure c P of 0.7 MPa is applied to the confining hole wall in the above simulations. In the KBS-3 concept the bentonite-soil mixture in the emplacement borehole could generate a swelling pressure of approximately 5 MPa [20]. Hence, numerical simulations on pillar damage evolution at the heating stage with confining pressures of 0, 2, 3, 4, and 5 MPa were carried out to study the impact of confining pressure on damage zone distribution ( Figure 14). The initial in

Effect of Confining Pressure on Damage Zone Evolution
In underground repositories for radioactive waste, clay bentonite mixtures are usually used to backfill the deposit holes to seal off the canisters. This kind of bentonite mixture swells considerably when exposed to water, resulting in swelling pressure being applied to the rock walls. As shown in Figure 7, a confining pressure P c of 0.7 MPa is applied to the confining hole wall in the above simulations. In the KBS-3 concept the bentonite-soil mixture in the emplacement borehole could generate a swelling pressure of approximately 5 MPa [20]. Hence, numerical simulations on pillar damage evolution at the heating stage with confining pressures of 0, 2, 3, 4, and 5 MPa were carried out to study the impact of confining pressure on damage zone distribution ( Figure 14). The initial in situ stress values of σ 1 and σ 3 were 42 MPa and 14 MPa, respectively, corresponding to Step 70 in Figure 9c. Figure 14 shows the damage zone distribution at day 60 with confining pressure P c of 0.7 and 0 MPa. When confining pressure P c decreases from 0.7 MPa to 0 MPa, more tensile damage occurs at the hole walls, leading to the formation of a deeper v-shaped notch. The shear damage zone in the area between deposit holes and heating holes also increases a little. (d) total damage zone area in simulation; (e) accumulated AE events each day [34].

Effect of Confining Pressure on Damage Zone Evolution
In underground repositories for radioactive waste, clay bentonite mixtures are usually used to backfill the deposit holes to seal off the canisters. This kind of bentonite mixture swells considerably when exposed to water, resulting in swelling pressure being applied to the rock walls. As shown in Figure 7, a confining pressure c P of 0.7 MPa is applied to the confining hole wall in the above simulations. In the KBS-3 concept the bentonite-soil mixture in the emplacement borehole could generate a swelling pressure of approximately 5 MPa [20]. Hence, numerical simulations on pillar damage evolution at the heating stage with confining pressures of 0, 2, 3, 4, and 5 MPa were carried out to study the impact of confining pressure on damage zone distribution ( Figure 14). The initial in situ stress values of 1  and 3  were 42 MPa and 14 MPa, respectively, corresponding to Step 70 in Figure 9c.  Figure 15a shows the variations of the damage zone area over time, in which the shear damage zone area, the tensile damage zone area, and the total damage zone area are given. For confining pressure P c of 0.7 MPa, the heating-induced total damage zone is 0.43 m 2 , including shear and tensile damage zone areas of 0.39 and 0.04 m 2 , respectively. For a confining pressure P c of 0 MPa, the heating-induced total damage zone is 0.55 m 2 , including shear and tensile damage zone areas of 0.47 and 0.08 m 2 , respectively. When confining pressure increases from 0 to 0.7 MPa, the total, shear, and tensile damage zone areas decrease by about 22%, 17%, and 50%, respectively. This indicates that the shear and tensile damage zones are both suppressed by the confining pressure. the heating-induced total damage zone is 0.55 m 2 , including shear and tensile damage zone areas of 0.47 and 0.08 m 2 , respectively. When confining pressure increases from 0 to 0.7 MPa, the total, shear, and tensile damage zone areas decrease by about 22%, 17%, and 50%, respectively. This indicates that the shear and tensile damage zones are both suppressed by the confining pressure. According to Andersson et al. [34], at the heating stage in APSE, the confining pressure began to decrease from day 61 on. During the release of the confining pressure, many AEs were recorded but they died off within two days of the final pressure drop. The numerical simulation results from the FLAC3D model also show that when the confining pressure increases from 0 to 0.8 MPa the crack initiation stress increases by about 4% and the damaged area also becomes smaller [22]. Hence, the conclusion shown in Figure 15 is consistent with the abovementioned measurement [34] and simulations [22]. According to Andersson et al. [34], at the heating stage in APSE, the confining pressure began to decrease from day 61 on. During the release of the confining pressure, many AEs were recorded but they died off within two days of the final pressure drop. The numerical simulation results from the FLAC3D model also show that when the confining pressure increases from 0 to 0.8 MPa the crack initiation stress increases by about 4% and the damaged area also becomes smaller [22]. Hence, the conclusion shown in Figure 15 is consistent with the abovementioned measurement [34] and simulations [22].

Effect of Biot's Coefficient on Damage Zone Evolution
As shown in Equation (4) and Figures 5 and 7, water flow is also considered in the proposed THMD model. The Biot's coefficient α was set to zero in all the above simulations to focus on the effects of excavation, heating, and confining pressure on damage evolution. Hence, a numerical simulation of pillar damage evolution at the heating stage with a Biot's coefficient α of 0.5 was carried out. The initial in situ stress values of σ 1 and σ 3 were 42 MPa and 14 MPa, respectively, corresponding to Step 70 in Figure 9c. Figure 16 is the damage zone distribution at day 60 with Biot's coefficients α of 0 and 0.5. When α increases from 0 to 0.5, since the effect of water pressure on pillar deformation is based on the effective stress principle, more tensile and shear damage appears at the hole walls, leading to a deeper v-shaped notch.

Discussion
In this paper, a damage-based model for rock under coupled thermal-hydrologicalmechanical conditions is proposed, in which the heterogeneity of rock is also characterized using Weibull distribution. The proposed model is then implemented in, and solved by, COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment.
The simulation results for the damage zone distribution at the excavation stage ( Figure 9) and the heating stage ( Figure 11) show that the location and depth of the v-shaped notch are mainly controlled by the trend and magnitude of in situ stress ( 1  , 3  ) respectively, while heating-induced stress is mainly responsible for the damage zone located in the area between the deposit holes and heating holes. The damage zone distribution with in situ stress ( 1  , 3  ) trends of 320, 230 shows relatively good agreement with the AE location. In addition, the v-shaped notch under in situ stress ( 1  , 3  ) of (36 MPa, 12 MPa), corresponding to the 3.5 m depth section in Figure 10, is still not formed at the end of the heating stage (Figure 11a, Step 60). The main reasons for the above discrepancies between our simulations and monitoring are likely: (1) simplification from a three-dimensional jointed rock mass to a two-dimensional heterogeneous rock; (2) simplifications of  Figure 17 shows the variations of damage zone area with time, in which the shear damage zone area, tensile damage zone area, and total damage zone area are given. For a Biot's coefficient α of 0, the heating-induced total damage zone is 0.43 m 2 , including shear and tensile damage zone areas of 0.39 and 0.04 m 2 , respectively. For a Biot's coefficient α of 0.5, the heating-induced total damage zone is 0.50 m 2 , including shear and tensile damage zone areas of 0.45 and 0.05 m 2 , respectively. When α increases from 0 to 0.5, the total, shear, and tensile damage zone areas increase by about 16%, 15%, and 25%, respectively. This indicates that the shear and tensile damage zones are both promoted by the water pressure.

Discussion
In this paper, a damage-based model for rock under coupled thermal-hydrologicalmechanical conditions is proposed, in which the heterogeneity of rock is also characterized using Weibull distribution. The proposed model is then implemented in, and solved by, COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment.

Discussion
In this paper, a damage-based model for rock under coupled thermal-hydrological-mechanical conditions is proposed, in which the heterogeneity of rock is also characterized using Weibull distribution. The proposed model is then implemented in, and solved by, COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment.
The simulation results for the damage zone distribution at the excavation stage ( Figure 9) and the heating stage ( Figure 11) show that the location and depth of the v-shaped notch are mainly controlled by the trend and magnitude of in situ stress (σ 1 , σ 3 ) respectively, while heating-induced stress is mainly responsible for the damage zone located in the area between the deposit holes and heating holes. The damage zone distribution with in situ stress (σ 1 , σ 3 ) trends of 320, 230 shows relatively good agreement with the AE location. In addition, the v-shaped notch under in situ stress (σ 1 , σ 3 ) of (36 MPa, 12 MPa), corresponding to the 3.5 m depth section in Figure 10, is still not formed at the end of the heating stage (Figure 11a, Step 60). The main reasons for the above discrepancies between our simulations and monitoring are likely: (1) simplification from a three-dimensional jointed rock mass to a two-dimensional heterogeneous rock; (2) simplifications of the horizontal in situ stress directions; (3) absence of the middle principal stress in the damage criterion; and (4) simplification of thermal and hydraulic boundary conditions.
Although the damage zone area that developed each day shows good agreement with the AE events, as shown in Figure 13, especially during days 1-20 and 40-60, the simulation result during days 21-39 is relatively different from the monitoring. After the temperature increasing stage from day 1 to 19, the temperature around the deposit holes reaches about 22 • C and basically remains stable during days 21-39. A relatively small number of AE events are recorded during this time. However, in our simulation no damage occurs at the same time. It is likely that the actual damage process of rock, especially of jointed and heterogeneous rock masses, takes a relatively long time to complete, including damage in local mesoscopic minerals, stress redistribution, permeability variation, water pressure, and rock temperature variation. This time-dependent effect of the rock mass, however, is not considered in our simulations. Figure 15b shows that when confining pressure P c increases from 0 to 0.7, 2, 3, 4, and 5 MPa, the total damage zone area decreases from 5.5 to 4.3, 3.8, 3.5, 3.4, and 3.2 m 2 , decreasing by about 22%, 30%, 37%, 39%, and 41%, respectively. This indicates that the heating-induced damage zone area decreases as the confining pressure increases, but the decreasing amplitude declines gradually as the confining pressure increases from 0 to 5 MPa. Considering the fact that the confining pressure applied on the deposit hole wall can counteract temperature increase-induced swelling, the finding in Figure 15b is understandable, because confining pressure's inhibition on damage development is related to the original increment of temperature and thermal stress. Figure 17 indicates that shear and tensile damage are both promoted by water pressure. The increments of total, shear, and tensile damage are 0.07, 0.06, and 0.01 m 2 , respectively. There is a relatively larger increment of shear damage because the effective stress has the effect of translating the Mohr circle to the left by the amount αp. The increment of tensile damage, however, is relatively small as a result of the compression induced by the in situ stress. Hence, the deformation and reopening of pre-existing fractures and the impact on water pressure distribution will have a great influence on damage evolution, which will be the subject of future research.
Notwithstanding the model limitations mentioned above, the results of the presented modeling work indicate that the application of numerical simulation tools provides a valuable insight regarding the evaluation and prediction of underground nuclear waste repositories under coupled THM conditions. To overcome the given limitations step-by-step is an essential part of future work [17].