A Full Coupled Thermal–Hydraulic–Chemical Model for Heterogeneity Rock Damage and Its Application in Predicting Water Inrush

: A coupled thermal–hydraulic–chemical (THC) model was carried out in this paper to study the inﬂuence of rock heterogeneity and the coupling e ﬀ ect of temperature, groundwater, and hydrochemistry on rock damage. Firstly, the hydrochemical and hydraulic erosion equations were established. The equations of the coupled THC model were established by combining the hydrochemical and hydraulic erosion equations, the ﬂow equations, and the heat transfer equations. Weibull distribution was adopted to govern the heterogeneity of initial rock porosity distribution. Secondly, the inﬂuence of the hydrochemistry, the temperature and the initial porosity heterogeneity on porosity and ﬂuid velocity change was studied. Then the rock damage rule changed with time at di ﬀ erent pH values and temperature was studied. Finally, an actual deep coal mine model was established to apply the THC model to predict water inrush. Results indicate that: (1) The average porosity and average ﬂuid velocity approximately show linear growth and exponential growth with time, respectively, and their growth rates increase with decreasing pH value and increasing temperature in a certain acidity and temperature range. (2) The increase of initial porosity heterogeneity has little inﬂuence on porosity change, but it can increase the ﬂuid velocity growth rate. The porosity heterogeneity and ﬂuid velocity heterogeneity approximately show exponential growth with increasing time, and the rock heterogeneity growth contributes to form cracks. The increase of temperature and decrease of pH value have little inﬂuence on the porosity heterogeneity, but they can increase the growth rate of the ﬂuid velocity heterogeneity. (3) The rock damage shows linear growth with time, and its growth rate increases with decreasing pH value and increasing temperature in a certain acidity range and temperature range. (4) The increase of rock heterogeneity can increase the possibility of water inrush.


Introduction
As the shallow coal resource is increasingly exhausted in China, more and more coal resources will be mined in deep mines. The depth of deep mines is generally over 600 m. Mine water inrush means that groundwater in the aquifer flows into roadways or working face through the water flowing channels. Mine water inrush is a major threat to mine production safety. Predicting water inrush is an important measure to guarantee mine workers' security. Water inrush is more likely to occur under the interaction of high ground temperature, high ground water pressure, high ground stress, and a complex hydrochemical environment in deep mines. The impact of a single factor of high temperature [1], water pressure [2], or hydrochemistry [3] has great influence on rock properties, and their coupled processes will have a greater and more complex influence on rock properties. Therefore, it is necessary to study the coupled thermal-hydraulic-chemical (THC) effect on rock damage to predict water inrush. The framework of international projects like DEvelopment of COupled models and their VALidation against EXperiment (DECOVALEX) have dramatically promoted the study of the thermal-hydraulic-mechanical (THM) coupling process [4,5]. Based on the mass, energy and momentum conservation, and number of experiments, these researches have made great progress in the combined effect of fluid flow and heat transfer on rock or soil properties [6][7][8][9]. For instance, Xiao et al. [8] developed a THM coupling model to simulate the combined effect of fracture fluid flow and heat transfer on matrix stability and shearing dilation behaviors. Based on these models and experiments, the rock damage rule was applied to deep mining and predicting water inrush [10,11]. In recent years, considering the hydrochemical erosion, the research on coupled hydraulic-mechanical (HC) [12], THC [13], mechanical-hydraulic-mechanical (MHC) [14,15] and thermal-hydraulic-mechanical-chemical (THMC) [16][17][18][19] effect on rock damage has made great progress. For instance, Bond et al. [19] presented a study that included 2D and 3D high-resolution coupled thermo-hydro-mechanical-chemical models to study the impact of complex physical and chemical processes on the geological formation surrounding the nuclear waste disposal facility. However, the rock properties were generally homogeneous in previous multi-physics coupling works especially in coupled THC numerical simulation models. A series of research results have indicated that rock heterogeneity has a great influence on rock damage. In an inhomogeneous vesicle size, the microcracks grow from the largest vesicles [20]. The heterogeneity has a great influence on rock failure mode, and the complexity of crack propagation increases with increasing rock heterogeneity [21]. Therefore, it is necessary to establish a multi-physics coupling model for heterogeneity rock. Before this, a kind of method should be found to govern the rock heterogeneity. In recent years, Weibull distribution is usually adopted to govern the rock heterogeneity [22][23][24][25]: Zhu et al. [23] used Weibull distribution to describe rock heterogeneity and discerned that the homogeneity index was the most important factor to simulate material failure; Wang et al. [24] assumed that rock parameters follow Weibull distribution and proposed a time-dependent creep model to study time-dependent creep behavior of heterogeneous brittle rocks; Zheng et al. [25] adopted Weibull statistical distribution to assign rock heterogeneity to find an effective method to increase gas drainage efficiency. Therefore, in the process of establishing a multi-physics coupling model for heterogeneity rock, Weibull distribution can be adopted to govern the rock heterogeneity.
In summary, rock heterogeneity has a great influence on rock damage, but the rock properties were generally homogeneous in previous multi-physics coupling works especially in multi-physics coupling numerical simulation models. Meanwhile, the rock of aquifuge was also generally simulated by homogeneous material in previous models for predicting water inrush, which could make mining engineers underestimate the rock damage degree and the danger degree of water inrush. The aquifuge is a kind of rock matrix which neither transmits nor stores water, and is located between the coal seam and aquifer, which prevents groundwater from the aquifer from flowing into the coal seam. In view of this, a coupled thermal-hydraulic-chemical (THC) model is carried out in this paper to study the influence of rock heterogeneity and the coupling effect of temperature, groundwater, and hydrochemistry on rock damage. Weibull distribution is adopted to govern the heterogeneity of initial rock porosity distribution. Considering the chemical reaction between the solid skeleton and the fluid, the porous rock is viewed as a four-phase medium consisting of solid skeleton, fluid, the fluidized solids from hydraulic erosion and the fluidized solids from chemical erosion. The equations of the coupled THC model consist of the hydrochemical and hydraulic erosion equations, the flow equations, and the heat transfer equations. Two kinds of numerical models are established. One is a cylinder model, with a size of ϕ50 × 100 mm, the standard specimens size that rock uniaxial and triaxial compression tests usually adopt. The other one is an actual deep coal mine model. The cylinder model is established to study the basic behavior of coupled THC processes for heterogeneity rock damage. The actual deep coal mine model is established based on the results of the basic behavior to predict water inrush of coal seam 15 of the Huatai coal mine, a deep coal mine in China. These models show the change rule of rock heterogeneity, which indicates the influence of hydrochemistry, temperature, and the initial porosity heterogeneity on porosity and fluid velocity change, and reveals that heterogeneity has an important influence on rock damage. These findings could provide a new understanding of the impact of rock heterogeneity on rock damage for related researchers and make mining engineers heighten their vigilance to the dangers of water inrush induced by rock heterogeneity in an aquifuge.

The Equations of Coupled THC Model
Considering the chemical reaction between the solid skeleton and the fluid, the porous rock is viewed as a four-phase medium, and the hydrochemical and hydraulic erosion equations are established based on this four-phase medium. The equations of the coupled THC model consist of the hydrochemical and hydraulic erosion equations, the flow equations, and the heat transfer equations.

The Hydrochemical and Hydraulic Erosion Equation
According to the previous research [26], the porous rock was viewed as a three-phase medium consisting of solid skeleton, fluid, and fluidized solids. The fluidized solids came from the solid skeleton and generated due to the hydraulic erosion. However, the contribution of hydrochemical erosion to the fluidized solids was not considered and the chemical reaction between the solid skeleton and the fluid should be considered. Therefore, the fluidized solids are divided into two parts including the fluidized solids from hydraulic erosion and the fluidized solids from chemical erosion. Based on this four-phase medium, the establishment process of the hydrochemical and hydraulic erosion equation is as follows.

Basic Equations
The porosity and the volume fraction of the fluidized solids are defined as where φ is the porosity of the porous rock; V, V p , V f , V hs and V cs are the volume of the porous rock, pore, fluid, the fluidized solids from hydraulic erosion, and the fluidized solids from chemical erosion, respectively, m 3 , in which the volume of the fluidized solids from hydraulic erosion and chemical erosion are actually the reduced volume of the solid skeleton due to hydraulic erosion and chemical erosion; c h and c c are the volume fraction of the fluidized solids from hydraulic erosion and the fluidized solids from chemical erosion, respectively. Therefore, the equivalent densities of the four phases are defined as where ρ 1 , ρ 2 , ρ 3 and ρ 4 are the equivalent density of the solid skeleton, fluid, the fluidized solids from hydraulic erosion, and the fluidized solids from chemical erosion, respectively, kg/m 3 ; m s , m f , m hs and m cs are the mass of the solid skeleton, fluid, the fluidized solids from hydraulic erosion, and the fluidized solids from chemical erosion, respectively, kg; ρ s and ρ f are the solid density and fluid density, respectively, kg/m 3 ; V s is the volume of the solid skeleton, m 3 . The mass balance for the solid skeleton can be expressed as where v sk is the fluid velocity of the solid skeleton, m/s; I hr and R are the mass sources and sinks, kg/(m 3 ·s), I hr represents the removed mass due to the hydraulic erosion here, and R represents the removed mass due to the chemical erosion here. R also represents the chemical reaction rate. The mass balance for the fluidized solids from hydraulic erosion can be expressed as where v f s is the fluid velocity of the fluidized solids from hydraulic erosion, m/s; I hp is the mass sources and sinks, kg/(m 3 ·s), which represent the produced mass from the solid skeleton. Therefore, According to the previous research [26], c h 1 (13) where λ is the porosity diffusivity coefficient, m 2 /s.

The Flow Equation
The continuity equation for flow in porous medium is The Darcy's law is where ρ f is the fluid density, kg/m 3 ; v is the Darcy fluid velocity, m/s; k is the permeability, m 2 ; µ is the fluid dynamic viscosity, N·s/m 2 ; p is the fluid pressure, Pa; z is a unit vector in the vertical direction; g is the gravitational acceleration, m/s 2 ; and Q s is the volumetric flow rate per unit volume of reservoir for a fluid source, 1/s.

The Heat Transfer Equation
According to the previous research [27], the thermal conductivity and thermal conductivity of porous medium are defined by an averaging model to account for both solid matrix and fluid properties, which are where k e f f is the effective thermal conductivity, W/(m·K); k s and k f are the thermal conductivity of solid matrix and fluid, respectively, W/(m·K); k d is the dispersion thermal conductivity, W/(m·K); (ρ f c f ) e f f is the effective volumetric heat capacity, J/(m 3 ·K); c s is the specific heat capacity of solid material in porous media, J/(kg·K); c f is the specific heat capacity of fluid, J/(kg·K). If the dispersion thermal conductivity is not considered, Equation (19) becomes The governing equations of heat transfer for porous media are where T is the absolute temperature, K; q is the conductive heat flux, J/(s·m 2 ), Q is the heat source or sink, W/K 3 . In summary, the governing equation of hydrochemical and hydraulic erosion is Equation (16), the governing equations of fluid flow are Equations (17) and (18), and the governing equation of heat transfer is Equation (22).
Compared with the full coupled THMC processes reported by Othman [28], some of the coupled THMC processes and some new coupled processes are studied using the equations of the coupled THC model, as shown in Figure 1. Appl. Sci. 2019, 9,

Numerical Model Setup
The coupled THC numerical model is established by using COMSOL Multiphysics 5.1 (COMSOL Co., Ltd., Shanghai, China). Porous rock is assumed to be a continuous medium. The fluid properties are assumed to be constant. The dispersion thermal conductivity is not considered. The porosity heterogeneity is adopted to describe the rock heterogeneity, and assumed to follow the Weibull distribution. The geometric model is a cylinder, with a size of φ50 × 100 mm, the standard specimens size that rock uniaxial and triaxial compression tests usually adopt. The bottom of the cylinder is applied fixed constraint, the side surface and top of the cylinder are applied boundary pressure to simulate the confining pressure and axial pressure, respectively. The bottom and top of the cylinder are applied to the inlet and outlet water pressure, respectively.
The probability density function of the Weibull distribution for initial porosity distribution is where w λ is the scale parameter, which governs the average porosity; w k is the shape parameter, which governs the porosity heterogeneity. The random numbers for Weibull distribution are generated by the external functions of COMSOL, then they are used to define the porosity in the material setting function of COMSOL. When w λ is 0.12, Figure 2 shows the initial porosity distribution at different shape parameters. Apparently, the porosity heterogeneity decreases with increasing the shape parameter. In addition, the computer time is the random number seeds, therefore, the calculation results of the random numbers for the Weibull distribution are different each time. Nevertheless, the number of points in a porosity range is approximately same at different times, as shown in Figure 3. More and more points concentrate around 0.11 with an increasing shape parameter, and all the average initial porosity values at different shape parameters and times are 0.11, which indicates that the scale parameter has no influence on the average porosity. Therefore, when the scale parameter is constant, setting different shape parameters can simulate the rocks that have the same average porosity but different heterogeneity.

Numerical Model Setup
The coupled THC numerical model is established by using COMSOL Multiphysics 5.1 (COMSOL Co., Ltd., Shanghai, China). Porous rock is assumed to be a continuous medium. The fluid properties are assumed to be constant. The dispersion thermal conductivity is not considered. The porosity heterogeneity is adopted to describe the rock heterogeneity, and assumed to follow the Weibull distribution. The geometric model is a cylinder, with a size of ϕ50 × 100 mm, the standard specimens size that rock uniaxial and triaxial compression tests usually adopt. The bottom of the cylinder is applied fixed constraint, the side surface and top of the cylinder are applied boundary pressure to simulate the confining pressure and axial pressure, respectively. The bottom and top of the cylinder are applied to the inlet and outlet water pressure, respectively.
The probability density function of the Weibull distribution for initial porosity distribution is where λ w is the scale parameter, which governs the average porosity; k w is the shape parameter, which governs the porosity heterogeneity. The random numbers for Weibull distribution are generated by the external functions of COMSOL, then they are used to define the porosity in the material setting function of COMSOL. When λ w is 0.12, Figure 2 shows the initial porosity distribution at different shape parameters. Apparently, the porosity heterogeneity decreases with increasing the shape parameter. In addition, the computer time is the random number seeds, therefore, the calculation results of the random numbers for the Weibull distribution are different each time. Nevertheless, the number of points in a porosity range is approximately same at different times, as shown in Figure 3.
More and more points concentrate around 0.11 with an increasing shape parameter, and all the average initial porosity values at different shape parameters and times are 0.11, which indicates that the scale parameter has no influence on the average porosity. Therefore, when the scale parameter is constant, setting different shape parameters can simulate the rocks that have the same average porosity but different heterogeneity.
(a) (b) (c) Figure 2. The initial porosity distribution at different shape parameters: The permeability is generally a function of porosity [26], which is where c k is a permeability parameter independent of the porosity of the rock, m 2 ; in φ is the initial porosity.
The model parameters are listed in Table 1. According to the previous research, the pH value of fluid and temperature have influence on the chemical reaction rate between the fluid and rock. However, there are various qualitative and quantitative analyses and different conclusions about the chemical reaction rate. In general, the chemical reaction rate increases with the decrease of the pH value of groundwater [29], and with the increase of the temperature, the chemical reaction rate increases at first and then decreases [3,30]. In this paper, according to [29,31], the chemical reaction rate increased by 3.06 × 10 −6 kg/m 3 /s when the temperature increased by 1 °C, and increased 4.62 × 10 −5 kg/m 3 /s when the pH value decreased by 1, and the chemical reaction rate was 6.97 × 10 −5 kg/m 3 /s when the pH value was 6.6 and the temperature was 75 °C.   The permeability is generally a function of porosity [26], which is where k c is a permeability parameter independent of the porosity of the rock, m 2 ; φ in is the initial porosity.
The model parameters are listed in Table 1. According to the previous research, the pH value of fluid and temperature have influence on the chemical reaction rate between the fluid and rock. However, there are various qualitative and quantitative analyses and different conclusions about the chemical reaction rate. In general, the chemical reaction rate increases with the decrease of the pH value of groundwater [29], and with the increase of the temperature, the chemical reaction rate increases at first and then decreases [3,30]. In this paper, according to [29,31], the chemical reaction rate increased by 3.06 × 10 −6 kg/m 3 /s when the temperature increased by 1 • C, and increased 4.62 × 10 −5 kg/m 3 /s when the pH value decreased by 1, and the chemical reaction rate was 6.97 × 10 −5 kg/m 3 /s when the pH value was 6.6 and the temperature was 75 • C.

Results
Based on the above analysis, setting different temperatures, pH values, and shape parameters are used to study the impact of temperature, hydrochemistry, and rock heterogeneity on rock porosity, fluid velocity, rock cracks, and rock damage. Then the above model is computed using COMSOL, and the following results are obtained.

The Porosity
The influence of the hydrochemistry, the temperature and the initial porosity heterogeneity on porosity change is obtained in this section.

Results
Based on the above analysis, setting different temperatures, pH values, and shape parameters are used to study the impact of temperature, hydrochemistry, and rock heterogeneity on rock porosity, fluid velocity, rock cracks, and rock damage. Then the above model is computed using COMSOL, and the following results are obtained.

The Porosity
The influence of the hydrochemistry, the temperature and the initial porosity heterogeneity on porosity change is obtained in this section. In order to observe and analyze the porosity distribution conveniently, the porosity values at 692 points at each time step are sorted according to the ascending order in Figure  5. In these figures, the porosity values at most points increase with increasing time and pH value. Meanwhile, the difference between the minimum and maximum values of porosity at the 692 points increases with increasing time, which indicates that the heterogeneity of porosity increases with increasing time.   The average and variance of porosity and fluid velocity at the 692 points are shown in Figure 6. Figure 6a-c show the average porosity and porosity variance at different shape parameters. The variance change can reflect the porosity heterogeneity change. In this figure, the curves between average porosity and time are approximately straight lines, and their slopes increase with the increase of the pH value, which indicates that the increase of the pH value can increase the porosity growth rate. At different shape parameters, the curves have little change, which indicates that the heterogeneity of the initial porosity has little influence on porosity change. The porosity variance firstly experiences a slight then rapid increase, and approximately shows exponential growth with increasing time. Therefore, the rock heterogeneity increases with increasing time. In addition, by increasing the shape parameter, the increased range of the variance decreases. Specifically, after 5 × 10 6 s, when the pH value was 6.6, and w k was 2, the variance increased from 2.09 × 10 −4 to 3.49 × 10 −3 (increased by 3.28 × 10 −3 ); when w k was 5, the variance increased from 4.12 × 10 −5 to 1.35 × 10 −3 (increased by 1.31 × 10 −3 ); when w k was 8, the variance increased from 1.77 × 10 −5 to 4.37 × 10 −4 (increased by 4.20 × 10 −4 ). The porosity variance values at different pH values are listed in Table 2. As illustrated in this table, the pH value has little influence on the porosity variance change, which indicates that the pH value has little influence on rock heterogeneity change. The average and variance of porosity and fluid velocity at the 692 points are shown in Figure 6. Figure 6a-c show the average porosity and porosity variance at different shape parameters. The variance change can reflect the porosity heterogeneity change. In this figure, the curves between average porosity and time are approximately straight lines, and their slopes increase with the increase of the pH value, which indicates that the increase of the pH value can increase the porosity growth rate. At different shape parameters, the curves have little change, which indicates that the heterogeneity of the initial porosity has little influence on porosity change. The porosity variance firstly experiences a slight then rapid increase, and approximately shows exponential growth with increasing time. Therefore, the rock heterogeneity increases with increasing time. In addition, by increasing the shape parameter, the increased range of the variance decreases. Specifically, after 5 × 10 6 s, when the pH value was 6.6, and k w was 2, the variance increased from 2.09 × 10 −4 to 3.49 × 10 −3 (increased by 3.28 × 10 −3 ); when k w was 5, the variance increased from 4.12 × 10 −5 to 1.35 × 10 −3 (increased by 1.31 × 10 −3 ); when k w was 8, the variance increased from 1.77 × 10 −5 to 4.37 × 10 −4 (increased by 4.20 × 10 −4 ). The porosity variance values at different pH values are listed in Table 2. As illustrated in this table, the pH value has little influence on the porosity variance change, which indicates that the pH value has little influence on rock heterogeneity change.      The porosity heterogeneity growth indicates that the porosity has an evolution of the polarity between 0 and 1. In order to better observe the cracks, it is assumed that the points with over 0.5 of porosity value can form the cracks. Figure 7 shows the distribution of porosity over 0.5 after 5 × 10 6 s when k w is 2, 5, and 8. In this figure, the area formed by these points with porosity over 0.5 can be considered as the cracks. Regarding the porosity evolution of the polarity between 0 and 1, it can be considered that some pore sizes decrease and some increase. The points that trend to 1 of porosity value contribute to form cracks, which is corresponding to the results reported by Heap et al. [20]: in an inhomogeneous vesicle size, the microcracks grow from the largest vesicles. The cracks increase with a decreasing shape parameter, which indicates that the rock heterogeneity growth contributes to rock fracture and increases the complexity of cracks, and is corresponding to the results reported by Tang and Liu [21]: the heterogeneity has a great influence on rock failure mode, and the complexity of crack propagation increases with increasing rock heterogeneity. The porosity heterogeneity growth indicates that the porosity has an evolution of the polarity between 0 and 1. In order to better observe the cracks, it is assumed that the points with over 0.5 of porosity value can form the cracks. Figure 7 shows the distribution of porosity over 0.5 after 5 × 10 6 s when w k is 2, 5, and 8. In this figure, the area formed by these points with porosity over 0.5 can be considered as the cracks. Regarding the porosity evolution of the polarity between 0 and 1, it can be considered that some pore sizes decrease and some increase. The points that trend to 1 of porosity value contribute to form cracks, which is corresponding to the results reported by Heap et al. [20]: in an inhomogeneous vesicle size, the microcracks grow from the largest vesicles. The cracks increase with a decreasing shape parameter, which indicates that the rock heterogeneity growth contributes to rock fracture and increases the complexity of cracks, and is corresponding to the results reported by Tang and Liu [21]: the heterogeneity has a great influence on rock failure mode, and the complexity of crack propagation increases with increasing rock heterogeneity.

The Influence of Temperature
When k w is 2, and the pH value is 6, the distribution of porosity at different a temperature and time is shown in Figure 8. This figure indicates that the porosity of most parts increases with increasing temperature and time in a certain temperature range. When k w is 2, the average porosity and porosity variance of the 692 points at different pH values and temperature are shown in Figure 7d

The Fluid
The influence of the hydrochemistry, the temperature, and the initial porosity heterogeneity on the fluid velocity is obtained in this section. Affected by the porosity change, the fluid velocity change rule is similar to the porosity change rule in some aspects, but there are still some differences between them in many aspects.

The Fluid
The influence of the hydrochemistry, the temperature, and the initial porosity heterogeneity on the fluid velocity is obtained in this section. Affected by the porosity change, the fluid velocity change rule is similar to the porosity change rule in some aspects, but there are still some differences between them in many aspects.

The Influence of Hydrochemistry
When k w is 2 and T is 75 • C, the distribution of porosity at different pH values and time is shown in Figures 9 and 10

The Fluid
The influence of the hydrochemistry, the temperature, and the initial porosity heterogeneity on the fluid velocity is obtained in this section. Affected by the porosity change, the fluid velocity change rule is similar to the porosity change rule in some aspects, but there are still some differences between them in many aspects.  Compared with the average and variance of the porosity, the average and variance of the fluid velocity at the 692 points have different change rules. As shown in Figure 7g-i, the average fluid velocity approximately shows exponential growth with increasing time, and the growth rate increases with a decreasing pH value from 6.6 to 4. When w k increases from 2 to 8, after 5 × 10 6 s, the average fluid velocity at pH 4 increases from 4 × 10 −3 m/s to 0.93 m/s, 0.84 m/s, and 0.83 m/s, respectively, which indicates that the average fluid velocity growth rate has a little increase with the increase of the heterogeneity of the initial porosity. As shown in Figure 7d-f, the fluid velocity variance also shows exponential growth with increasing time, and the growth range is larger than that of porosity variance in the same condition. For instance, when w k is 2, the fluid velocity variance at pH 6 increases from 3.13 × 10 −6 m/s to 2.26 × 10 −2 m/s (increased by about 7.21 × 10 3 times) from 0 s to 5 × 10 6 s, but in the same condition, as shown in Figure 7a- when the pH values are 6.6, 6, 5, and 4, respectively. Therefore, although the pH value has little influence on the growth rate of the porosity variance (as shown in Table 2), the increase of the pH value can increase the growth rate of the fluid velocity variance. Compared with the average and variance of the porosity, the average and variance of the fluid velocity at the 692 points have different change rules. As shown in Figure 7g-i, the average fluid velocity approximately shows exponential growth with increasing time, and the growth rate increases with a decreasing pH value from 6.6 to 4. When k w increases from 2 to 8, after 5 × 10 6 s, the average fluid velocity at pH 4 increases from 4 × 10 −3 m/s to 0.93 m/s, 0.84 m/s, and 0.83 m/s, respectively, which indicates that the average fluid velocity growth rate has a little increase with the increase of the heterogeneity of the initial porosity. As shown in Figure 7d-f, the fluid velocity variance also shows exponential growth with increasing time, and the growth range is larger than that of porosity variance in the same condition. For instance, when k w is 2, the fluid velocity variance at pH 6 increases from 3.13 × 10 −6 m/s to 2.26 × 10 −2 m/s (increased by about 7.21 × 10 3 times) from 0 s to 5 × 10 6 s, but in the same condition, as shown in Figure 7a-c, the porosity variance increases from 2.09 × 10 −4 to 3.40 × 10 −3 (increased by about 16.24 times). In addition, after 5 × 10 5 s, the fluid velocity variance increases from 3.13 × 10 −6 m/s to 8.31 × 10 −3 m/s, 2.26 × 10 −2 m/s, 1.29 × 10 −1 m/s, and 7.89 × 10 −1 m/s (increased by 2.65 × 10 3 times, 7.21 × 10 3 times, 4.13 × 10 3 m/s times, and 2.52 × 10 5 times) when the pH values are 6.6, 6, 5, and 4, respectively. Therefore, although the pH value has little influence on the growth rate of the porosity variance (as shown in Table 2), the increase of the pH value can increase the growth rate of the fluid velocity variance.

The Influence of Temperature
When k w is 2, and the pH value is 6.6, the distribution of fluid velocity at a different temperature and time is shown in Figure 11. This figure indicates that the fluid velocity at most parts increases with increasing temperature and time in a certain temperature range. When k w is 2, the average fluid velocity and fluid velocity variance at the 692 points at different pH values and temperature are shown in Figure 7j-l. The curves between average fluid velocity and time shows exponential growth with increasing time, and the growth rate increases with increasing temperature, which indicates that the increase of temperature can increase the fluid velocity growth rate in a certain temperature range. The fluid velocity variance approximately also shows exponential growth with increasing time, and the variance has a wide varying change (e.g., after 5 × 10 5 s, increases from 3.13 × 10 −6 m/s to 4.49 × 10 −2 m/s, 1.28 × 10 −1 m/s, and 7.72 × 10 −1 m/s (increased by 1.43 × 10 4 times, 4.09 × 10 4 times, and 2.47 × 10 5 times) at 90 • C when the pH values are 6.6, 5, and 4, respectively). Therefore, the decrease of the pH value can increase the fluid velocity variance. Meanwhile, the fluid velocity variance growth rate increases with increasing temperature.

The Influence of Temperature
When w k is 2, and the pH value is 6.6, the distribution of fluid velocity at a different temperature and time is shown in Figure 11. This figure indicates that the fluid velocity at most parts increases with increasing temperature and time in a certain temperature range. When w k is 2, the average fluid velocity and fluid velocity variance at the 692 points at different pH values and temperature are shown in Figure 7j-l. The curves between average fluid velocity and time shows exponential growth with increasing time, and the growth rate increases with increasing temperature, which indicates that the increase of temperature can increase the fluid velocity growth rate in a certain temperature range. The fluid velocity variance approximately also shows exponential growth with increasing time, and the variance has a wide varying change (e.g., after 5 × 10 5 s, increases from 3.13 × 10 −6 m/s to 4.49 × 10 −2 m/s, 1.28 × 10 −1 m/s, and 7.72 × 10 −1 m/s (increased by 1.43 × 10 4 times, 4.09 × 10 4 times, and 2.47 × 10 5 times) at 90 °C when the pH values are 6.6, 5, and 4, respectively). Therefore, the decrease of the pH value can increase the fluid velocity variance. Meanwhile, the fluid velocity variance growth rate increases with increasing temperature.

The Rock Damage
The porosity change has been widely used to define the rock damage [32]. In this study, the rock damage is defined as where D is the rock damage; φ is the average porosity; in φ is the initial average porosity.
When w k is 2, the curves of rock damage change with time and their fitting polynomials at different pH values and temperature are shown in Figure 12. In this figure, the rock damage increases with time, which can be expressed as where f D is the fitting rock damage; t is time, 10 6 s; d k is the rock damage coefficient, which reflects the rock damage growth rate. Figure 12 shows that the rock damage coefficient increases with a decreasing pH value from 6.6 to 4 and increasing temperature from 75 °C to 100 °C, which indicates that in a certain acidity range and temperature range, the rock damage growth rate increases with decreasing pH value and increasing temperature.

The Rock Damage
The porosity change has been widely used to define the rock damage [32]. In this study, the rock damage is defined as where D is the rock damage; φ is the average porosity; φ in is the initial average porosity. When k w is 2, the curves of rock damage change with time and their fitting polynomials at different pH values and temperature are shown in Figure 12. In this figure, the rock damage increases with time, which can be expressed as where D f is the fitting rock damage; t is time, 10 6 s; k d is the rock damage coefficient, which reflects the rock damage growth rate. Figure 12 shows that the rock damage coefficient increases with a decreasing pH value from 6.6 to 4 and increasing temperature from 75 • C to 100 • C, which indicates that in a certain acidity range and temperature range, the rock damage growth rate increases with decreasing pH value and increasing temperature.

Application in Predicting Water Inrush
In the previous work of predicting water inrush in deep coal mines, the rock properties were generally homogeneous. In this section, based on the actual engineering background of coal seam 15 in the Huatai coal mine of China, applying the above results, the influence of rock heterogeneity on water inrush is studied. As shown in Figure 13, the depth of coal seam 15 is 1180 m, the length of the working face is 140 m, and the distance between the aquifer and the coal seam is about 85 m. The angle of the coal seam is 15°. The fractured area is induced by the coal mining. Groundwater can flow into the coal seam through the fractured area. Based on the statistical analysis of long-term experimental and engineering measured data, the researchers and engineers found that the thickness of the fractured area can be calculated by some functions of the length of the working face, the depth of the coal seam, and the angle of the coal seam, and therefore obtained the empirical equations. The empirical formulas are [33]:

Application in Predicting Water Inrush
In the previous work of predicting water inrush in deep coal mines, the rock properties were generally homogeneous. In this section, based on the actual engineering background of coal seam 15 in the Huatai coal mine of China, applying the above results, the influence of rock heterogeneity on water inrush is studied. As shown in Figure 13, the depth of coal seam 15 is 1180 m, the length of the working face is 140 m, and the distance between the aquifer and the coal seam is about 85 m. The angle of the coal seam is 15 • . The fractured area is induced by the coal mining. Groundwater can flow into the coal seam through the fractured area. Based on the statistical analysis of long-term experimental and engineering measured data, the researchers and engineers found that the thickness of the fractured area can be calculated by some functions of the length of the working face, the depth of the coal seam, and the angle of the coal seam, and therefore obtained the empirical equations. The empirical formulas are [33]:

Application in Predicting Water Inrush
In the previous work of predicting water inrush in deep coal mines, the rock properties were generally homogeneous. In this section, based on the actual engineering background of coal seam 15 in the Huatai coal mine of China, applying the above results, the influence of rock heterogeneity on water inrush is studied. As shown in Figure 13, the depth of coal seam 15 is 1180 m, the length of the working face is 140 m, and the distance between the aquifer and the coal seam is about 85 m. The angle of the coal seam is 15°. The fractured area is induced by the coal mining. Groundwater can flow into the coal seam through the fractured area. Based on the statistical analysis of long-term experimental and engineering measured data, the researchers and engineers found that the thickness of the fractured area can be calculated by some functions of the length of the working face, the depth of the coal seam, and the angle of the coal seam, and therefore obtained the empirical equations. The empirical formulas are [33]:  The conception of the water inrush coefficient is proposed based on the statistical analysis of long-term water inrush data and defined in the Regulation for Coal Mine Water Prevention and Control [34], where K is the water inrush coefficient, MPa/m; P is the water pressure, MPa; M is the thickness of the aquifuge, m. If the water inrush coefficient is over 0.1 MPa/m, the water inrush will occur. In the 1970s, considering the thickness of the disturbed area and the height of water intrusion into the aquifuge, the water inrush coefficient was revised as where M e f f is the thickness of the effective aquifuge, m; h 2 is the height of water intrusion into aquifuge, m. The height of water intrusion into the aquifuge is 0 m in this case. The water pressure is 6.17 MPa calculated by Equations (32) and (33). The bottom boundary of the aquifer is applied 6.17 MPa, and no flow is applied to the left, right, and top boundaries of the geometric model. From the results of previous analysis in Section 3, when the scale parameter is constant, setting different shape parameters can simulate the rocks that have the same average porosity but different heterogeneity, and the larger the shape parameter is, the more homogeneous the rock. Therefore, the scale parameter is set as 0.12, the shape parameters are set as 2 and 8, respectively, to simulate the different initial rock heterogeneity. The shape parameter should be as large as possible to simulate homogeneity rock, which is set as 50. The model parameters are listed in Table 3. These parameters are from the rock compression test. The rock specimens are from the roof and floor rock strata of coal seam 15. The testing machine is MTS815.03. In addition, the effective thermal conductivity and the effective volumetric heat capacity are calculated by Equations (21) and (20), and the permeability is calculated by Equation (25), and these three parameters are governed by the porosity generated by the Weibull distribution random numbers. The thermal conductivity and heat capacity of the rock and fluid properties parameters are listed in Table 2. Because the scale parameter is constant, when the shape parameters are 2, 8, and 50, the average porosity at 619 points in an effective aquifuge are all 0.2, and the average fluid velocity at these points are all 6 × 10 −4 m/s. Figure 14 shows the fluid velocity distribution of homogeneity rock and different heterogeneity rock, in which, Figure 14b,d,f show the fluid velocity over average value of 6 × 10 −4 m/s. In this figure, the fluid velocity heterogeneity and the maximum fluid velocity increases with decreasing the shape parameters. Meanwhile, the distribution area of fluid velocity over the average value gradually extends to the disturbed area by decreasing the shape parameter, which will increase the possibility of water inrush. In the results section, the rock heterogeneity increases with increasing temperature. Therefore, the danger of water inrush will increase with the increase of rock heterogeneity due to the combined effect of high temperature and high water pressure in deep mines.

Discussions
This coupled THC model is mainly used to study rock heterogeneity and the coupling effect of temperature, groundwater, and hydrochemistry on rock damage, and the results indicate that the rock heterogeneity growth contributes to the formation of cracks and increases the possibility of water inrush. Our findings show that the points with a larger porosity value contribute to form cracks, which is corresponding to the finding that the microcracks grow from the largest vesicles in an

Discussions
This coupled THC model is mainly used to study rock heterogeneity and the coupling effect of temperature, groundwater, and hydrochemistry on rock damage, and the results indicate that the rock heterogeneity growth contributes to the formation of cracks and increases the possibility of water inrush. Our findings show that the points with a larger porosity value contribute to form cracks, which is corresponding to the finding that the microcracks grow from the largest vesicles in an inhomogeneous vesicle size reported by Heap et al. [20]. Meanwhile, rock heterogeneity growth contributes to rock fracture and increases the complexity of cracks, which is corresponding to the finding that the heterogeneity has a great influence on the rock failure mode, and the complexity of crack propagation increases with increasing rock heterogeneity reported by Tang and Liu [21]. Certainly, compared with previous heterogeneity rock models that adopted the Weibull distribution to govern the heterogeneity [22][23][24][25], some new findings of the coupling effect of temperature, groundwater, and hydrochemistry on rock heterogeneity, rock damage, and fluid velocity are obtained when THC coupling processes are considered. These findings could provide a new understanding of the impact of rock heterogeneity on rock damage for related researchers and heighten the vigilance of mining engineers to the dangers of water inrush induced by rock heterogeneity in aquifuge.
However, there are some limitations to this study which are listed as follows: 1.
In alkaline conditions and high temperature that make the chemical reaction rate decrease, the rock damage change rule is not studied. In addition, the impact of thermal expansion on rock damage is not considered.

2.
The relationship between the external stress and porosity change is not established. Therefore, the change in porosity due to mechanical compression or expansion is not studied. The rock damage is only defined by the porosity change. The strength parameters (e.g., the elastic modulus), are not considered to define the damage.
Future research should consider the mechanical-induced thermal, hydraulic, and chemical processes, as well as the alkaline condition and higher temperature to define and study rock damage. In addition, it is better to associate the actual heterogeneity obtained by imaging techniques with the artificial distribution. Future research should use some imaging techniques, such as Computed Tomography (CT) or Nuclear Magnetic Resonance (NMR) imaging techniques, to show the actual rock particle and pore size distribution, then based on the imaging results, adopt some statistical approaches to represent the actual rock homogeneity or heterogeneity in numerical modeling.

Conclusions
A coupled THC model was carried out in this paper to study the rock heterogeneity and the coupling effect of temperature, groundwater, and hydrochemistry on rock damage. The main conclusions are obtained:

1.
With increasing time, the average porosity approximately shows linear growth, and the average fluid velocity approximately shows exponential growth, and their growth rates increase with decreasing pH value and increasing temperature in a certain acidity and temperature range.

2.
The Weibull distribution is adopted to govern the heterogeneity of the initial rock porosity distribution, and the variance is used to express the heterogeneity. The increase of initial porosity heterogeneity has little influence on porosity change, but it can increase the fluid velocity growth rate. The porosity heterogeneity and fluid velocity heterogeneity approximately show exponential growth with increasing time. The increase of temperature and decrease of pH value have little influence on the porosity heterogeneity, but they can increase the growth rate of the fluid velocity heterogeneity. The rock heterogeneity growth contributes to rock fracture, which is corresponding to the results reported by the previous research.

3.
The porosity change is used to define the rock damage. The rock damage growth rate increases with a decreasing pH value and increasing temperature in a certain acidity range and temperature range.

4.
The water inrush is more likely to occur through the rock with greater heterogeneity.