###### Abstract
[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]: [PERSON] 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; [PERSON] 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; [PERSON] 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 \(\varphi 50\times 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.