Prediction of Fracture Evolution and Groundwater Inrush from Karst Collapse Pillars in Coal Seam Floors: A Micromechanics- Based Stress-Seepage-Damage Coupled Modeling Approach

State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China State Key Laboratory of Coal Resources and Safe Mining, School of Mines, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China School of Resources and Environment Engineering, Jiangxi University of Science and Technology, Ganzhou 341000, China


Introduction
A karst collapse pillar (KCP) is a karst cavity formed by the strong dissolution of groundwater in the underlying soluble limestone of coal measures. The overlying rock mass of a KCP is unstable and collapses under gravity and other effects, which then fills the dissolution space, leading to a cylindrical fractured rock mass pillar [1], as shown in Figure 1. A KCP is the result of a special geological phenomenon that occurs widely in the coalfields of North China. KCPs generally have a water-conducting section and a water-resisting section. Because of the existence of the water-resisting section, the vast majority (more than 80%) of KCPs exposed by mining are non-water-conducting [2].
However, under the combined action of mininginduced stress and floor-confined hydraulic pressure, the initial fractures in the water-resisting section of a KCP may propagate, connect, and coalesce, thereby forming a water-conducting channel, which converts the original non-water-conducting KCP into a water-conducting KCP, eventually causing groundwater inrush incidents from the coal seam floor [3][4][5]. Since the 1960s, a number of groundwater inrush accidents in China's coal mines have caused numerous casualties, property loss, and severe damage to the water resource environment in mining areas. One of the most typical cases causing the most severe loss involved groundwater inrush from a hidden KCP in the floor of the second-level working face in the Fangezhuang coal mine of the Kailuan mining area, China, on June 2, 1984. The maximum water inflow reached 2,053 m 3 /min, flooding a large-scale mechanized coal mine with an annual output of greater than 3 Mt in less than 20 h, causing the deaths of nine miners and resulting in a direct economic loss of several hundred million RMB. Therefore, a scientific and accurate understanding of the occurrence and mechanisms of groundwater inrush from KCPs in coal seam floors is of great significance to ensure the safe mining of coal resources in confined groundwater bodies [6][7][8][9].
In recent years, many researchers have conducted extensive studies on the mechanisms of groundwater inrush from KCPs. Wang and Yang [10] investigated the groundwater inrush processes of strong water-filling KCPs, non-water-filling KCPs, and hidden KCPs in the floor using FLAC 3D numerical simulations. Bai [11] established a "plug" mechanical model accommodating coupled fluid-solid interactions to study the mechanism of groundwater inrush from KCPs. Zhu and Wei [9] employed COMSOL software to simulate the process and mechanism of groundwater inrush from KCPs. Yao et al. [12] established a dynamic model of fluid-solid coupling with a variable mass for fractured rock masses and studied the permeability evolution caused by particle migration of a KCP. Ma et al. [13] studied groundwater inrush disasters caused by seepage instabilities in KCPs. Zuo et al. [14,15] studied the failure behavior of KCPs using a discrete element method. Ma et al. [16] investigated the particle erosion effect on the groundwater inrush mechanism of a KCP by using FLAC 3D software and a stress-seepage coupling model. These previous studies have provided a scientific basis for the prevention and control of groundwater inrush from KCPs in coal seam floors. However, existing studies on the mechanism of groundwater inrush from KCPs have mostly been based on the traditional elasticplastic mechanics theory or phenomenological damage mechanics theory. These studies were limited in their ability to reproduce the microstructural damage and resulting permeability evolution behavior of rocks. Therefore, a realistic process and the mechanisms of groundwater inrush from KCPs in coal seam floors are far from well understood.
Field observation data have shown that groundwater inrush from a KCP in a coal seam floor is caused by the propagation and coalescence of numerous initial microfractures and the induced permeability enhancement in the KCP under the combined action of mining-induced stress and confined hydraulic pressure [17,18]. Therefore, the process of groundwater inrush from a KCP essentially involves the multiscale coupling of rock stress, damage, fracture, and seepage, which cannot be analyzed in depth by existing rock mechanics theories or numerical calculation methods. The establishment of an effective theoretical model and a numerical method to accurately accommodate the fracture and permeability evolution and seepage characteristics of KCPs under hydromechanical coupling conditions is a key scientific problem involved in the current study of the mechanism of groundwater inrush.
In this study, a novel micromechanical model incorporating the stress-seepage-damage (SSD) coupling of rock is proposed to study the development of damage, fractures, permeability, and seepage in a coal seam floor with a KCP. Numerical results are verified by in situ high-precision microseismic monitoring observations. Furthermore, the key influencing factors, including the advancing distance of the working face, the confined hydraulic pressure, and the initial water-conducting height of the KCP, on the process and mechanism of groundwater inrush from a KCP are systematically discussed.  Figure 2.

Micromechanical
It is assumed that each microcrack is an ideal threedimensional (3D) coin-shaped or two-dimensional (2D) linear crack and that these microcracks are randomly distributed in the MREV. The microcracks in the MREV are divided into N groups according to their different inclinations. The initial geometric parameters of each group of microcracks primarily include the following: (i) the initial

Geofluids
statistical mean radius of the microcracks (a 0 ), which is the average of the initial lengths of all microcracks in the same group; (ii) the mean inclination angle of the microcracks (β), which is the average of the inclination angles of all microcracks in this group; and (iii) the number density of the microcracks (m), which is calculated by the ratio of the number of all microcracks in this group to the volume of the MREV. These initial geometric parameters are assumed to be constant for all the microcrack families so that the initial response of the MREV is isotropic. These initial microcrack parameters can be determined by optical measurement tests combined with digital image processing technology [19]. The MREV is subjected to the far-field total stress σ ij ði, j = 1, 2, 3Þ, and the hydraulic pressure in the microcracks and micropores is p. Under the combined action of σ ij and p, the microcracks in the MREV will undergo initiation, propagation, connection, and coalescence (i.e., microscopic damage), and this evolution is the root cause of the changes in the mechanical properties and seepage properties of the MREV. The governing equations of the SSD coupling for the MREV are derived in the following sections. The stress symbol follows the elastic mechanics rule (the tensile stress is positive, and the pore water pressure is positive).

Balance Equations.
For the deformation of the solid skeleton of the MREV, it is assumed that the total stress tensor is σ and the total strain tensor is ε. Considering that the microscopic damage (microcracking) in the rock does not affect the equilibrium equation and geometric equation of the solid, the following equations apply: where F is the body force tensor and u is the displacement tensor.
For the water seepage in the MREV, according to the mass balance principle, the difference in the mass of the fluid flowing into and out of the MREV is equal to the change in the fluid mass within the MREV as follows: where q is the seepage velocity in the rock and ζ is the change in the fluid content per unit volume of the rock. Assuming that the water seepage flow in the MREV follows Darcy's law, the seepage velocity can be expressed as follows: where υ is the dynamic viscosity coefficient, p is the water pore pressure, D denotes the damage tensor of the MREV (closely related to the microcracking behavior), and K represents the equivalent permeability tensor of the MREV (which is a function of the rock damage tensor variation). Furthermore, it is assumed that in any state of the rock damage process, the coupled stress-seepage response of the MREV satisfies Biot's theory of poroelasticity [20]. Based on this assumption, the poroelasticity coefficients of the MREV can be regarded as functions of the rock damage tensor as follows: where CðDÞ is the fourth-rank elastic stiffness tensor of the damaged MREV and exhibits the classic Voigt symmetries C ijkl = C jikl = C ijlk = C klij ði, j, k, l = 1, 2, 3Þ, MðDÞ denotes the scalar Biot's modulus of the damaged MREV, and αðDÞ is the symmetric (α ij = α ji ) second-rank Biot's effective stress coefficient tensor of the damaged MREV. Equation (4) shows that the damaged MREV is an anisotropic porous elastic medium, whose Biot's coefficients (i.e., α and M) can be expressed as functions of the anisotropic elastic stiffness coefficient (C) as follows [21][22][23]: where δ ij ði, j = 1, 2, 3Þ is the Kronecker symbol, K s is the bulk modulus of the solid constituent of the rock, K f is the bulk modulus of the fluid, which is approximately 3.3 GPa for water, ϕ is the porosity of the rock, and K * is the generalized drainage bulk modulus of the damaged rock, where K * ðDÞ = C iijj ðDÞ/9. Based on the Voigt symmetry of the elastic stiffness tensor, Biot's effective stress coefficient tensor α ij ðDÞ in Equations (1)-(4) are the general governing equations derived for the SSD coupling model of the MREV of the rock. Compared with a classic stress-seepage (SS) coupling model, this model is coupled with the influence of damage, thereby establishing a theoretical foundation for the in-depth study of the rock failure process and the corresponding flow evolution under hydromechanical coupling conditions. The specific relations among the damage tensor D, the elastic stiffness tensor CðDÞ, and the permeability tensor KðDÞ of the rock are further derived in the following sections based on microcracking damage mechanics.

Constitutive Equations of Microcracking Damage
Evolution. The initial state of the MREV (before it is loaded) is the undamaged state. When the MREV is subjected to hydromechanical coupling loading, its internal microcracks will undergo initiation, propagation, connection, and coalescence (i.e., microcracking damage). To characterize the effect of microcracking damage, a second-order symmetric damage tensor defined by the change in the density of microcracks in the MREV [24] is used as follows: where N is the number of groups of microcracks in the MREV, m k and n k (k = 1 ⋯ N) are the number density and the normal unit vector of the kth group of microcracks, respectively, and d k is defined as the relative change in the density of the kth group of microcracks induced by microcrack growth, expressed as follows: where a 0 and a k are the initial average length of the kth group of microcracks and the length after propagation under external loading, respectively. Equations (7) and (8) directly link the damage tensor to the length of microcrack growth in each orientation, which is determined according to the dynamic conditions and laws of microcracking. In current studies, a classic model used to analyze macroscopic fracture propagation is the sliding wing crack model [25,26]. However, this model cannot realistically reproduce microcracking behavior. A large number of microscopic experiments in rocks have shown that microcracks in rocks rarely propagate in the wing crack pattern, but they frequently take on more complex forms [27]. In this study, to facilitate the establishment of the constitutive equations of microcracking damage evolution using continuum mechanics, the equivalent microcrack propagation model proposed by Costin [28] and Golshani et al. [29] is used. It is assumed that there is an equivalent local tensile stress field around a microcrack and that this field induces a self-similar mode of tensile propagation of the microcrack. The magnitude of the local tensile stress field is dependent on the farfield deviatoric stress and the confining pressure around the microcrack is as follows: whereσ t is the equivalent local tensile stress acting on the surface of the kth group of microcracks, σ n and S n are the far-field normal stress and deviatoric stress around the kth group of microcracks, respectively, δ is the Kronecker symbol, and f is a scalar proportional function related to the microcrack length a k as follows [24]: where ω is a dimensionless constant reflecting the characteristics of microcrack propagation and b is the critical length of microcrack propagation. These parameters can all be determined by conventional triaxial compression tests [24]. Under the combined action of the equivalent local tensile stress on the microcrack surface and the water pore pressure in the microcracks, the propagation criterion for the microcracks is obtained based on the linear elastic fracture mechanics theory as follows: where K I is the stress intensity factor of the microcrack and K IC is the mode I fracture toughness of the rock. Equation (11) is used to obtain the propagation length of each group of microcracks in the MREV of the rock in a given stress state, and then, the damage tensor of the MREV can be calculated using Equation (7). A modified Helmholtz free energy function [24,30] is used to represent the elastic potential function of the damaged MREV as follows: where ψ and ε are the Helmholtz free energy and strain tensor of the damaged MREV, respectively, λ and μ are the initial Lamé elastic constants of the rock, and A and B are the damage influence coefficients, which are introduced to describe the contribution of microcrack damage to the free energy of the MREV. By taking the first derivative of Equation (12) with respect to the strain tensor ε, the equivalent elastic stiffness tensor of the damaged MREV is explicitly expressed as follows: The anisotropic elastic stiffness coefficient matrix of the damaged MREV in Equation (13) contains 21 independent components, which can be expressed as functions of the damage tensor components as follows: Equations (7), (11), and (13) form the constitutive equations of microcrack damage evolution in the MREV, and these equations represent the damage state of anisotropic microcracks and the macroscopic constitutive response of rocks under hydromechanical coupling.

Equations of the Permeability Enhancement Induced by
Microcracking Damage. Microcracking damage also significantly enhances the permeability of rocks. The MREV of the rock is regarded as a porous matrix medium with embedded microcracks. Therefore, the overall permeability of the MREV is the superposition of two parts as follows: where K is the overall permeability tensor of the MREV of the rock, K 0 is the permeability tensor of the porous matrix medium, and K c is the additional permeability tensor caused by the microcracking damage.
For the kth group of microcracks in the MREV, the additional permeability tensor caused by this group of microcracks (under a 2D condition) can be derived from a modified cubic law and Darcy's law as follows: where κ is the modification factor of the cubic law, V is the volume of the MREV of the rock, and E 0 is the elastic modulus of the porous matrix of the rock.
By superimposing Equation (16) for each group of microcracks in the MREV, the total permeability tensor of the MREV induced by microcracking damage is as follows: where R k ða k Þ is a connectivity coefficient reflecting the connectivity of the microcrack network and is represented as a function of the microcrack length as follows [31]: where t 1 and t 2 are material constants. Equations (15)-(17) establish a relation between the equivalent permeability tensor and the microcrack damage tensor. The equations show that the resulting permeability tensor is proportional to the (t 2 + 4)th power of the microcrack length, indicating that the propagation of microcracks will cause dramatic changes in the permeability of the MREV. Although an extensive number of previous studies have proposed different theoretical models to describe the permeability evolution [32], most of these models reflect only the effect of the stress state and cannot accommodate the influence of changes in the microstructures of the rock. The microcrack-based permeability evolution model proposed in this study compensates for these deficiencies, thus laying a foundation for further investigation of the process of groundwater inrush from KCPs in coal seam floors.

Numerical Implementation of the SSD Coupling Model.
All of the mechanical and hydraulic parameters (e.g., the elastic stiffness coefficients and the permeability coefficients) in the proposed SSD coupling model are dynamically adjusted and change with the development of the damage. This SSD coupling model is a highly nonlinear fluid-solid coupling problem. Therefore, an iterative algorithm is proposed to approximately solve this problem using the following basic procedure: (i) Using a finite element mesh, the model geometry is discretized into a series of MREVs, and the hydromechanical parameters and microcrack parameters for each of the MREVs are initialized (ii) The loading condition of the model is discretized into several small subload steps in the time domain, and the model is loaded sequentially by subload (iii) During each subload step, all hydromechanical parameters of the rock are constant, and a finite element method based on the full coupling analysis is adopted to calculate the average stress field, average pore pressure field, average strain field, and average seepage velocity field for each MREV of the model (iv) According to the stress state of each MREV of the rock obtained by finite element calculation, the 5 Geofluids stress intensity factor at the microcrack tip of each group in each MREV is calculated by Equation (11), and the propagation length of each microcrack is calculated (v) Equations (7), (14), and (15) are used to compute the microcrack damage tensor, the new induced elastic stiffness tensor, and the permeability tensor, respectively, for each MREV of the rock (vi) The updated hydromechanical parameters are reintroduced into the fully coupled finite element analysis model in Step (iii), and the equilibrium calculation is performed again to reobtain the stress and seepage fields in the rock after microcrack damage. Steps (iii) to (vi) are repeated iteratively until none of the microcracks in the model for this subload step propagate (vii) The next subload step is applied, and Steps (iii) to (vii) are repeated until the full load is completed The above implementation procedure is programmed using a combination of MATLAB and COMSOL software to replicate the dynamic processes of damage and fracture as well as the seepage evolution behavior of the rock under hydromechanical coupling. This program enables more realistic and visual numerical simulations of the process and mechanism of groundwater inrush from KCPs in coal seam floors.

Model Validation with In Situ Microseismic
Monitoring Data 3.
1. An Engineering Geological Survey of a Working Site above a Confined Aquifer. In situ microseismic monitoring was carried out at the 10-108 working face of the 10# coal seam in the Tuanbai coal mine of the Huozhou mining area, Shanxi, China. This working face has a strike length of 900 m and a width of 180 m. The groundwater source for the 10# coal seam is mainly the Ordovician limestone aquifer in the floor. The mine is located in a complete hydrogeological unit consisting of the Fenhe fault basin and the surrounding Luliangshan and Huoshan Mountains. The Ordovician limestone is extensively developed in the exposed area of the Luliangshan Mountains in the west, and a large area of Ordovician limestone is exposed in the Fenhe River Valley from the north to the south. Quaternary sand gravel pore aquifers and Ordovician limestone aquifers are interconnected and replenish a large amount of seepage to the Ordovician limestone. Therefore, the Ordovician limestone aquifer under the 10# coal seam has abundant supply sources, which poses a threat to the mining of the 10# coal seam. The thickness of the aquiclude between the 10# coal seam floor and the underlying Ordovician limestone aquifer is approximately 27~45 m. The lithology of the aquiclude is primarily composed of aluminum mudstone and tight sandstone. Under normal conditions (without the influence of the tectonic fracture connection), the aquiclude has poor permeability and good water-resisting properties, thus playing an important role in the safe mining of the 10# coal seam above the confined aquifer. However, from the influence of multiple tectonic movements, the structures in the mine area are very developed, with an average of 91 faulted structures and 60 KCPs per km 2 , destroying the integrality and waterresisting properties of the aquiclude.
Previous 3D geological exploration observations showed that an X-187-08 KCP is located near the 10-80(2) roadway in the 10-108 working face of the mine, as shown in Figure 3. This KCP has an oblate cross section with an average radius of 20 m and a development height of 115 m, directly connecting the 10# coal seam and the underlying Ordovician limestone aquifer. Initially, the KCP does not conduct groundwater, but it may be activated by the mining effect of the 10-108 working face to allow groundwater inflow, posing a serious threat to safe mining at the working face. A coal mine adjacent to this mine experienced a groundwater inrush accident caused by a KCP in the floor of working face 2-1101 in March 2007, with a maximum water inflow of 1200 m 3 /h. The source of the groundwater inrush was the Ordovician limestone aquifer in the floor.
Therefore, to prevent groundwater inrush incidents in the floor of the 10-108 working face, it is very important to understand the mining-induced fracturing in the coal seam floor and the activation of the KCP during the process of advancing the working face. An in situ high-precision microseismic monitoring test was performed in the 10-108 working face to capture and locate the fracture characteristics of the coal seam floor and the KCP during mining, thereby providing a critical field basis for the prediction and prevention of groundwater inrush from the KCP in the coal seam floor.

In Situ Microseismic Monitoring Observations.
A Comise microseismic monitoring system, which is explosion-proof in underground mines, was used in this study for the in situ microseismic monitoring test, as shown in Figure 4. This system was supplied by Shandong University of Science and Technology in China. The microseismic monitor is equipped with several three-component microseismic detectors that convert the received microseismic signal (elastic wave signal) into an impulse signal and store it as a voltage. The arrival times of P waves through the microseismic detectors at different locations were collected and used to invert the temporal and spatial distribution of the microseismic source locations according to the Geiger algorithm [33].
The monitoring subjects are the floor of the 10-108 working face and the X-187-08 KCP. Six sets of microseismic monitoring boreholes (with diameters of 100 mm), labeled A, B, C, D, E, and F, were arranged in the 10-108(2) roadway, as shown in Figure 5. Two detectors, each with a diameter of 55 mm, were placed in each of the monitoring boreholes. Among them, borehole C collapsed from the construction in the KCP, such that the detector could not be placed. After the installation of the detectors, the boreholes were immediately grouted to fix the detectors and seal the boreholes. The five sets of microseismic monitoring boreholes (10 detectors in total) arranged in the 10-108(2) roadway formed an effective internal and external field monitoring space in the floor of the 10-108 working face. The microseismic monitoring 6 Geofluids host was placed 70 m in front of monitoring borehole F. By continuously detecting and recording the microseismic event signals from the floor, the fracture characteristic parameters of the floor and of the X-187-08 KCP during the mining process of working face 10-108 were identified. In Figure 6, plots of the profile projection diagram of the distribution of microseismic events in the floor along the strike of the working face in the mining process of the 10-108 working face are shown. The data show that the mining-induced microseismic events are densely distributed in the shallow region of the floor and become sparser in deeper regions. In the absence of KCPs, the fracture depth of the floor reaches a maximum value of 15.5 m during the periodic weighting of the roof (see points A and B in Figure 6). When the working face advances to the KCP region, the number of microseismic events in the floor increases significantly, and the fracture depth of the floor is greater than that of the normal intact floor. The deepest fracture zone in the floor is located at the edge of the KCP, with a maximum fracture depth of approximately 19.1 m (see points C and D in Figure 6). This result indicates that the aquiclude remains intact with an average minimum thickness of 16.9 m between the mining-induced fracture zone of the floor and the Ordovician limestone aquifer. Therefore, it is inferred that the mining of the 10-108 working face is relatively safe.

Validation of the Numerical Model.
Based on the engineering geological characteristics of the 10-108 working face and the X-187-08 KCP above, an idealized numerical model along the strike of the working face was developed to investigate the mining-induced fracture behavior in the floor with the KCP, as shown in Figure 7. The geometry of the model is 200 m × 100 m, with four idealized rock strata, including (from top to bottom) the overlying strata, the coal seam, the aquiclude, and the confined aquifer. The mechanical and hydraulic parameters for each of the strata are listed in Table 1. The boundary conditions correspond to a roller along the left, right, and lower sides of the calculation model and traction applied at the top boundary (σ y ). All boundaries are zero-flux boundaries except for the confined aquifer, where a constant water pressure (p 0 ) is applied.
In the numerical model, the KCP develops from the top of the confined aquifer to the bottom of the coal seam and includes a water-resisting section and a water-conducting section, with heights of h 1 and h 2 , respectively. The initial hydraulic pressure in the water-conducting section of the KCP is consistent with that of the confined aquifer. Since the mechanical properties of the filling materials vary randomly from point to point in an actual KCP (i.e., heterogeneity), the Weibull distribution function is selected to replicate the heterogeneity of the mechanical properties of the KCP as follows: where ε m denotes the characteristic mechanical parameters (e.g., elastic modulus) of the MREV in the KCP, ξ m is the average value of the characteristic parameters, and λ m denotes the coefficient of homogeneity. Relative to the KCP, the surrounding rock of the KCP can be assumed to be a homogeneous material. Table 2 lists the mechanical parameters of the KCP used in the present numerical simulation. The proposed calculation program for the SSD coupling model is employed to solve the above numerical model. Figure 8 shows the mining-induced fracture damage (D) distribution in the floor obtained by numerical simulation. In the figure, the scalar fracture damage variable (D) is calculated by the relative change in the fracture length as follows: where max ða k Þ denotes the maximum length of all microcracks in one MREV and D i represents the degree of fracture damage to the ith MREV. The data in Figure 8 show that for the intact floor, the maximum depth of the fracture damage zones reaches 16 m during the periodic weighting of the roof. When the working face advances to the KCP, the fracture damage of the floor is further aggravated. Because of the heterogeneity of the KCP, the development of fracture damage in the KCP is characterized by a random and disordered distribution. The depth of the fracture damage zones developed at the edge of the KCP is higher than that inside the KCP, with a maximum depth of approximately 21 m. Furthermore, Figure 9 compares the numerical observations and the microseismic monitoring results of the fracture depth of the floor, and these results mostly agree. This finding indicates that the proposed micromechanical modeling approach and the employed simplified 2D numerical model for the coal seam floor with the KCP are effective and reliable. Thus, they can be further used for in-depth study of the process and mechanism of groundwater inrush from a KCP in a coal seam floor. However, more realistic 3D numerical models based on high-performance computers that can represent the actual engineering situation need to be included in future works.

Numerical Simulations of Groundwater
Inrush from the KCP in the Coal Seam Floor    In addition, the damage region in the KCP is not directly linked vertically to the floor of the goaf, and the region shape is very irregular, which is mainly caused by the heterogeneity of the KCP and the propagation characteristics of the mining-induced stress (i.e., the mining-induced stress on the left side of the KCP is higher than that on the right side) (iii) The damage development and activation characteristics of the KCP show that groundwater inrush from the activated KCP in the floor does not necessarily occur when the KCP is exposed directly to the working face (i.e., the working face is linked to the KCP), but this process is very likely to occur at some distance away from the KCP To quantitatively characterize the degree of activation of the KCP, a variable D K is defined as the average of the fracture damage values of all the MREVs constituting the KCP as follows: where D i represents the fracture damage value of the ith MREV in the KCP and N is the total number of MREVs in the KCP. Figure 11 plots the resulting variation in the average degree of activation ( D K ) of the KCP with the distance (L) between the working face and the KCP. The data show that the average degree of activation ( D K ) of the KCP increases as a power function as L decreases. When L decreases from 6 m to 1 m, D K increases sharply from 0.016 to 0.47, an approximate 30-fold increase. This result indicates that as the working face approaches the KCP, the risk of groundwater inrush from the KCP increases significantly.
The mining-induced fracture further leads to changes in the seepage field of the KCP in the floor. Figure 12 shows the distributions of the pore pressure field in the floor at different distances of L = 20, 15, 11, 7, 5, 3, 2, and 1 m. Correspondingly, Figure 13 shows the distributions of the seepage vector field in the floor when L is 5, 3, 2, and 1 m. The data in the figures show the following:   through the water-resisting section of the KCP, which acts as a "rock plug" that seals the KCP and effectively blocks the hydraulic connection between the confined aquifer and the coal seam. However, as the working face gradually approaches the KCP, the KCP is continuously damaged and activated, and an increasing number of water-resisting regions in the KCP are converted into water-conducting regions, which gradually conduct the confined water upward along the KCP (iii) Groundwater inrush from a KCP is a gradual process rather than a mutation process. The progress of this groundwater inrush depends on a key inducing factor, the distance between the working face and the KCP, that reflects the mining influence of the working face (see Figure 13). Consequently, in actual projects, it is necessary to predict and prevent groundwater inrush from the KCP and to strengthen the dynamic monitoring of precursory information, such as the activation of the KCP and the conduction of confined water in the KCP Furthermore, Figure 14 shows that the average pore pressure (p K ) in the KCP varies with the depth of the floor (h) at distances of L = 7, 5, 3, and 1 m. p K represents the average pore water pressure at each cross section of the KCP. The data show that at a floor depth of 20 m < h < 24 m, the pore pressure in the KCP is always maintained at a constant value of p K = 5 MPa (equal to the initial confined hydraulic pressure) because this depth is located in the water-conducting section of the KCP. When the floor depth (h) is less than 20 m, the average pore pressure in the KCP decreases as h decreases (i.e., closer to the coal seam). This finding reflects the characteristic of gradual upward permeation of the confined water along the KCP. Moreover, at the same floor depth (h), as the distance (L) decreases, p K in the KCP increases. As the distance L decreases from 7 m to 1 m, p K at the top of the KCP (h = 0 m) increases sharply from 0.0079 MPa to 4.41 MPa, a nearly 500-fold increase. This finding also confirms that the distance between the working face and the KCP largely controls the process of groundwater inrush from the KCP.

Influence of the Confined Hydraulic
Pressure. Similar to the hydraulic fracturing mechanism of rock [18,34], the hydraulic pressure of the confined aquifer exerts an important influence on the damage and fracture of the coal seam floor. Figure 15 shows the fracture damage distribution of the floor under different confined hydraulic pressures (p 0 = 2:5 MPa and 7.5 MPa) during the working face advancement process (L = 8, 6, 4, and 2 m). Correspondingly, Figure 16 further plots the variation in the average degree of activation ( D K ) of the KCP with the confined hydraulic pressure (p 0 ) at different distances from the working face to the KCP. The data in the figures show the following:  12 Geofluids (i.e., the groundwater inrush channel) is formed in the KCP (see the right panel in Figure 15(b)). However, no obvious water-conducting channel appears in the KCP when p 0 = 2:5 MPa and L = 1 m (see the left panel in Figure 15(d)).  Figure 17 shows the distributions of the normalized pore pressure (p K /p 0 ) in the floor for different confined hydraulic pressures. Comparing the left and right groups of figures, the higher the confined hydraulic pressure is, the higher the water-conducting height of the confined water and the pore water pressure in the KCP, which further accelerate the damage and activation of the KCP. Thus, the risk of groundwater inrush from the KCP increases accordingly. For p 0 = 7:5 MPa and L = 6 m, the confined groundwater breaks through the blockage of the water-resisting section in the KCP and flows into the entire KCP (see the right panel in Figure 17(b)). However, for p 0 = 2:5 MPa and L = 2 m, the water-conducting height of the confined groundwater in the KCP is very small (see the left panel in Figure 17(d)).
In accordance with Figure 17, the variation in the normalized average pore pressure in the KCP ( p K /p 0 ) with the distance (L) under different confined hydraulic pressures (p 0 = 2:5, 5:0, and 7:5 MPa) is shown in Figure 18, where p K is defined as the average pore water pressure of all the MREVs that constitute the KCP. It can be seen that p K /p 0 gradually increases as L decreases. Moreover, when L  Figure   14 Geofluids decreases to a certain extent (see points A, B, and C in Figure 18), p K /p 0 exhibits a sudden and significant step increase. However, this observation does not imply that groundwater inrush occurs immediately at points A, B, and C because at this point, the water-conducting channel has not yet been connected to the damage and fracture zone of the floor (see Figure 17). In addition, for p 0 = 2:5, 5:0, and 7:5 MPa, a sudden increase in the average pore pressure in the KCP occurs at distances L = 3, 6, and 9 m, respectively. These results indicate that under a high confined hydraulic pressure, groundwater inrush may occur when the KCP is still far from the working face.

4.4.
Influence of the Initial Water-Conducting Height of the KCP. Figure 19 plots the distributions of the fracture damage of the floor under different initial water-conducting heights          damage in the KCP is intensive, and a macroscopic groundwater inrush channel is formed in the KCP at a distance of L = 2 m (see the right panels in Figure 19). Therefore, the greater the initial water-conducting height of the KCP, the more easily the KCP is damaged and activated. This is because an increase in the initial water-conducting height is equivalent to the thinning of the aquiclude and is more prone to induce damage and fracture of the floor. Furthermore, the variations in the average degree of activation ( D K ) of the KCP with the initial water-conducting height (h 2 ) at different distances from the working face to the KCP are shown in Figure 20. The data show that the degree of activation of the KCP ( D K ) increases with h 2 . Moreover, the larger the value of h 2 , the faster the degree of activation of the KCP ( D K ) develops with the advancement of the working face (L). When h 2 = 0 m, as L decreases from 8 m to 2 m, D K increases by only 0.046. In comparison, when h 2 = 8 m, under the same conditions, D K increases by approximately 0.14, which is approximately twice as much as the former increase. Figure 21 further shows the distributions of the pore pressure field in the floor under different initial waterconducting heights (h 2 = 0 m and 8 m) of the KCP. When h 2 = 0 m, the confined water in the floor eventually fails to conduct effectively in the KCP (see the left panels in Figure 21). However, when h 2 = 8 m, at a distance of L = 4 m, the confined water breaks through the blockage of the water-resisting section in the KCP and flows into the entire KCP (see the left panels in Figure 21). This finding indicates that the larger the initial water-conducting height of the KCP is, the easier it is for the KCP in the floor to evolve into a water-conducting channel and the higher the risk of water inrush from the KCP in the floor.
The variations in the normalized average pore pressure ( p K /p 0 ) in the KCP with the distance (L) from the working face to the KCP under different initial water-conducting heights (h 2 = 0, 4, and 8 m) of the KCP are shown in Figure 22. The data show that as the initial waterconducting height of the KCP increases, the pore pressure of the KCP significantly increases. Moreover, in the cases of h 2 = 4 m and 8 m, as L decreases to 8 m and 6 m, respectively (see points A and B in Figure 22), p K /p 0 shows a sudden stepwise increase, while in the case of h 2 = 0 m, there is no similar sudden increase in the average pore pressure in the KCP. These results indicate that a KCP with a larger initial waterconducting height is more likely to be damaged and activated to induce groundwater inrush.

Conclusions
In this study, a micromechanics-based SSD coupling model is developed and verified by in situ high-precision microseismic monitoring data. This micromechanical modeling approach is used to perform a series of numerical simulations of groundwater inrush from a KCP in a coal seam floor. The main conclusions are as follows: (i) Microcracking evolution in rock is the root cause of the damage, fracture, permeability enhancement, and seepage changes in the rock. Based on realistic physical mechanisms of microcracking in rock, a microcrack-based damage tensor has been introduced in the classic Biot's theory of poroelasticity to establish a novel SSD coupling model in which the macroscopic mechanical and hydraulic properties of the rock are directly linked to the microcrack kinetics. MATLAB programming is used to embed the proposed SSD coupling model in COMSOL software to create a dynamic simulation of rock damage, fracture, and seepage evolution under hydromechanical coupling conditions. This model overcomes the limitation of traditional models, i.e., that they are unable to quantitatively describe the influence of rock damage on the SS coupling process. The model also lays a foundation to further investigate the process and mechanism of groundwater inrush from KCPs in coal seam floors