## 3 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 \(\varphi 50\times 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  
\[\phi_{in}=\frac{k_{w}}{\lambda_{w}}\Bigg{(}\frac{\phi_{in}}{\lambda}\Bigg{)}^{ k_{w}-1}\exp\left(-\left(\frac{\phi_{in}}{\lambda_{w}}\right)^{k_{w}}\right) \tag{24}\]  
where \(\lambda_{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 \(\lambda_{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.  
Figure 1: The coupled thermal–hydraulic–chemical (THC) processes for heterogeneity rock damage.  
The permeability is generally a function of porosity [26], which is  
\[k=k_{c}\frac{{\phi_{in}}^{3}}{{(1-\phi_{in})}^{2}} \tag{25}\]  
where \(k_{c}\) is a permeability parameter independent of the porosity of the rock, m\({}^{2}\); \(\phi_{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\times 10^{-6}\) kg/m\({}^{3}\)/s when the temperature increased by 1 \({}^{\circ}\)C, and increased \(4.62\times 10^{-5}\) kg/m\({}^{3}\)/s when the pH value decreased by 1, and the chemical reaction rate was \(6.97\times 10^{-5}\) kg/m\({}^{3}\)/s when the pH value was \(6.6\) and the temperature was \(75\)\({}^{\circ}\)C.  
Figure 3: The initial porosity values at 692 points of the cylinder at three times: (**a**) \(k_{w}=2\); (**b**) \(k_{w}=5\); (**c**) \(k_{w}=8\).  
Figure 2: The initial porosity distribution at different shape parameters: (**a**) \(k_{w}=2\); (**b**) \(k_{w}=5\); (**c**) \(k_{w}=8\).