Abstract
Karst collapse pillars (KCPs) frequently cause severe groundwater inrush disasters in coal mining above a confined aquifer. An accurate understanding of the damage and fracture evolution, permeability enhancement, and seepage changes in KCPs under the combined action of mining-induced stress and confined hydraulic pressure is of great significance for the early prediction and prevention of groundwater inrush from KCPs in coal seam floors. In this study, a micromechanics-based coupled stress-seepage-damage (SSD) modeling approach, in which the macroscopic mechanical and hydraulic properties of the rock are explicitly related to the microcrack kinetics, is proposed to simulate the fracture evolution and the associated groundwater flow in KCPs. An in situ high-precision microseismic monitoring technology is used to verify the micromechanical modeling results, which indicate that the numerical model successfully reproduces the damage and fracture evolution in a coal seam floor with a KCP during the mining process. The presented model also provides a visual representation of the complex process of KCP activation and groundwater inrush channel formation. A numerical study shows that the damage and activation of a KCP start from the edge of the KCP, gradually develop toward the interior of the KCP, and eventually connect with the damage fracture zone of the floor, forming a primary water-conducting channel in the KCP, causing the confined groundwater to flow into the working face. Groundwater inrush from a KCP is a gradual process instead of a mutation process. A reduction in the distance between the working face and a KCP and increases in the confined hydraulic pressure and the initial water-conducting height of the KCP can significantly increase the risk of groundwater inrush from the KCP.
1. 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 mining-induced 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–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 m3/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–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 FLAC3D 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 FLAC3D 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 elastic-plastic 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.
2. Micromechanical Model of Stress-Seepage-Damage (SSD) Coupling of Rock
2.1. Microscopic Representative Element Volume (MREV) of Rock
In a rock mass, there are a large number of microscopic structures such as microcracks and micropores that fundamentally control the macroscopic mechanical and seepage properties of the rock. To facilitate the investigation of the overall hydromechanical response of rock containing numerous microscopic structures using continuum mechanics theory, an MREV, which mainly consists of a porous matrix medium (an equivalent, continuous, homogeneous, and isotropic elastic material) and embedded microcrack defects, is introduced at the microscopic scale, as shown in Figure 2.
It is assumed that each microcrack is an ideal three-dimensional (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 groups according to their different inclinations. The initial geometric parameters of each group of microcracks primarily include the following: (i) the initial statistical mean radius of the microcracks (), 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 (), 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 , and the hydraulic pressure in the microcracks and micropores is . Under the combined action of and , 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).
2.2. 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 is the body force tensor and 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 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, is the water pore pressure, denotes the damage tensor of the MREV (closely related to the microcracking behavior), and 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 is the fourth-rank elastic stiffness tensor of the damaged MREV and exhibits the classic Voigt symmetries , denotes the scalar Biot’s modulus of the damaged MREV, and is the symmetric () 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 ) can be expressed as functions of the anisotropic elastic stiffness coefficient () as follows [21–23]: where is the Kronecker symbol, is the bulk modulus of the solid constituent of the rock, is the bulk modulus of the fluid, which is approximately 3.3 GPa for water, is the porosity of the rock, and is the generalized drainage bulk modulus of the damaged rock, where .
Based on the Voigt symmetry of the elastic stiffness tensor, Biot’s effective stress coefficient tensor in Equation (5) and the generalized drainage bulk modulus are further expanded as follows:
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 , the elastic stiffness tensor , and the permeability tensor of the rock are further derived in the following sections based on microcracking damage mechanics.
2.3. 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 is the number of groups of microcracks in the MREV, and () are the number density and the normal unit vector of the th group of microcracks, respectively, and is defined as the relative change in the density of the th group of microcracks induced by microcrack growth, expressed as follows: where and are the initial average length of the th 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 far-field deviatoric stress and the confining pressure around the microcrack is as follows: where is the equivalent local tensile stress acting on the surface of the th group of microcracks, σn and Sn are the far-field normal stress and deviatoric stress around the th group of microcracks, respectively, is the Kronecker symbol, and is a scalar proportional function related to the microcrack length as follows [24]: where is a dimensionless constant reflecting the characteristics of microcrack propagation and 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 is the stress intensity factor of the microcrack and 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 and 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.
2.4. 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 is the overall permeability tensor of the MREV of the rock, is the permeability tensor of the porous matrix medium, and is the additional permeability tensor caused by the microcracking damage.
For the th 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, is the volume of the MREV of the rock, and 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 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 and 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 ()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.
2.5. 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 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.
3. 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 km2, destroying the integrality and water-resisting 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 m3/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.
3.2. 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 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.
(a)
(b)
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.
3.3. 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 , 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 (). All boundaries are zero-flux boundaries except for the confined aquifer, where a constant water pressure () 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 and , 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 denotes the characteristic mechanical parameters (e.g., elastic modulus) of the MREV in the KCP, is the average value of the characteristic parameters, and 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 () distribution in the floor obtained by numerical simulation. In the figure, the scalar fracture damage variable () is calculated by the relative change in the fracture length as follows: where denotes the maximum length of all microcracks in one MREV and represents the degree of fracture damage to the th 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.
4. Numerical Simulations of Groundwater Inrush from the KCP in the Coal Seam Floor
4.1. Simulation Setup
The 2D numerical model developed above is applied here to further investigate the effects of the distance between the working face and the KCP (), the hydraulic pressure of the confined aquifer (), and the initial water-conducting height of the KCP () on the fracture damage, seepage evolution, and groundwater inrush process of the floor with the KCP. Three schemes of numerical simulations are performed as follows: (i)Influence of : ; ; and (ii)Influence of : ; ; and(iii)Influence of : ; ; and
4.2. Influence of the Distance between the Working Face and the KCP
Figure 10 shows the distributions of the fracture damage of the floor when the working face is at different distances () from the KCP. The data in this figure show the following: (i)The advancement of the working face results in a remarkable fracture damage zone in the floor. However, damage regions first appear in the KCP only when the working face advances to less than 15 m from the KCP due to the influence of the advanced abutment pressure caused by mining. Then, as the working face advances further, the damage and activation of the KCP rapidly occur from the combined action of the mining-induced stress and the hydraulic pressure(ii)The damage and activation of the KCP initiate from the edge of the KCP (near the side of the working face), gradually develop toward the interior of the KCP, and eventually connect with the fracture damage zone of the 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 is defined as the average of the fracture damage values of all the MREVs constituting the KCP as follows: where represents the fracture damage value of the th MREV in the KCP and is the total number of MREVs in the KCP. Figure 11 plots the resulting variation in the average degree of activation () of the KCP with the distance () between the working face and the KCP. The data show that the average degree of activation () of the KCP increases as a power function as decreases. When decreases from 6 m to 1 m, 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 . Correspondingly, Figure 13 shows the distributions of the seepage vector field in the floor when is 5, 3, 2, and 1 m. The data in the figures show the following: (i)The KCP is the main channel for the seepage and conduction of confined water in the floor to the working face. The mining disturbance leads to damage and activation of the KCP, greatly improving the permeability of the KCP, such that a large amount of confined water gradually flows and gathers toward the region of the KCP, causing a significant increase in the pore pressure and the seepage velocity(ii)The variation in the seepage field in the KCP depends on the distance between the working face and the KCP. When the working face is far from the KCP, the confined water has difficulty flowing 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
(a)
(b)
(c)
(d)
Furthermore, Figure 14 shows that the average pore pressure () in the KCP varies with the depth of the floor () at distances of . represents the average pore water pressure at each cross section of the KCP. The data show that at a floor depth of , the pore pressure in the KCP is always maintained at a constant value of (equal to the initial confined hydraulic pressure) because this depth is located in the water-conducting section of the KCP. When the floor depth () is less than 20 m, the average pore pressure in the KCP decreases as 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 (), as the distance () decreases, in the KCP increases. As the distance decreases from 7 m to 1 m, at the top of the KCP () 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.
4.3. 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 ( and 7.5 MPa) during the working face advancement process (m). Correspondingly, Figure 16 further plots the variation in the average degree of activation () of the KCP with the confined hydraulic pressure () at different distances from the working face to the KCP. The data in the figures show the following: (i)At the same distance () from the working face to the KCP, the greater the hydraulic pressure in the confined aquifer is, the greater the degree of activation () of the KCP. In the case of , as increases from 2.5 MPa to 7.5 MPa, increases from 0.042 to 0.37, a nearly 10-fold increase. For and , a macroscopic fracture channel (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 and (see the left panel in Figure 15(d)).(ii)The greater the confined hydraulic pressure () is, the faster the development of the activation of the KCP with the advancing working face. When , as decreases from 8 m to 2 m, increases by only approximately 0.05. In comparison, when , increases by approximately 0.33 under the same conditions, a nearly sixfold increase
Furthermore, Figure 17 shows the distributions of the normalized pore pressure () 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 and , 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 and , 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 () with the distance () under different confined hydraulic pressures () is shown in Figure 18, where is defined as the average pore water pressure of all the MREVs that constitute the KCP. It can be seen that gradually increases as decreases. Moreover, when decreases to a certain extent (see points A, B, and C in Figure 18), 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 , a sudden increase in the average pore pressure in the KCP occurs at distances , 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 ( and 8 m) of the KCP during the working face advancement process (). The data show that when , the mining-induced fracture damage region in the KCP is sparse, and no water-conducting channel is formed in the KCP (see the left panels in Figure 19). However, when , the distribution of the fracture damage in the KCP is intensive, and a macroscopic groundwater inrush channel is formed in the KCP at a distance of (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 () of the KCP with the initial water-conducting height () 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 () increases with . Moreover, the larger the value of , the faster the degree of activation of the KCP () develops with the advancement of the working face (). When , as decreases from 8 m to 2 m, increases by only 0.046. In comparison, when , under the same conditions, 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 water-conducting heights ( and 8 m) of the KCP. When , the confined water in the floor eventually fails to conduct effectively in the KCP (see the left panels in Figure 21). However, when , at a distance of , 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 () in the KCP with the distance () from the working face to the KCP under different initial water-conducting heights () of the KCP are shown in Figure 22. The data show that as the initial water-conducting height of the KCP increases, the pore pressure of the KCP significantly increases. Moreover, in the cases of and 8 m, as decreases to 8 m and 6 m, respectively (see points A and B in Figure 22), shows a sudden stepwise increase, while in the case of , there is no similar sudden increase in the average pore pressure in the KCP. These results indicate that a KCP with a larger initial water-conducting height is more likely to be damaged and activated to induce groundwater inrush.
5. 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(ii)The developed model and calculation program are used to perform micromechanics-based numerical simulations of the patterns of coal seam floor fracture and KCP activation. Additionally, a high-precision coal mine explosion-proof microseismic detector is used to perform continuous and dynamic monitoring and recording of microseismic signals at the fracture of a coal seam floor in a mine with a KCP. The failure mode and failure depth of the coal seam floor along the strike of the working face obtained from the simulations are consistent with the microseismic monitoring results. Moreover, the numerical simulation results enable the creation of physically realistic visualizations of the fracture and damage evolution, as well as the seepage evolution pattern of the KCP in the floor during the mining process. The results also provide a true representation of the complex process of the formation of water-conducting channels and the evolution of groundwater inrush disasters in coal seam floors with KCPs(iii)The numerical simulation results of groundwater inrush from a KCP show that under the combined action of mining-induced stress and confined hydraulic pressure, the damage and activation of the KCP in the floor start from the edge of the KCP, gradually develop toward the interior of the KCP, and eventually connect with the fracture damage zone of the floor. The activation region of the KCP tends to be inclined toward the side of the working face during the mining process. Additionally, the mining-induced fractures greatly improve the permeability of the KCP, resulting in a significant increase in the pore pressure and seepage velocity of the KCP, and finally form a main water-conducting channel in the KCP that causes the confined groundwater to flow into the working face. Moreover, groundwater inrush from a KCP is a gradual process instead of a mutation process and commonly occurs at a certain distance away from the KCP. The progress of the groundwater is closely related to the distance between the working face and the KCP (), the confined hydraulic pressure (), and the initial water-conducting height of the KCP (). The smaller the value of and the larger value of the and are, the greater the degree of activation of the developing KCP with the advancement of the working face, the larger the pore pressure in the KCP, and the higher the risk of groundwater inrush from the KCP. Therefore, it is of great importance to strengthen the dynamic monitoring of precursory information, such as the activation of a KCP and the conduction of the confined water in the KCP, for the early forecasting and prevention of groundwater inrush from the KCP
In summary, the presented micromechanics-based SSD coupling model provides an alternative virtual experimental tool to perform realistic predictions of fracture evolution and associated groundwater flow from karst collapse pillars in coal seam floors above a confined aquifer. Moreover, it allows for replications of key coupling processes at the microscale in fluid-infiltrated rocks. The proposed numerical approach has broad application prospects in solving a variety of complicated hydromechanical problems in practical rock engineering where fracture-fluid interaction effects are important.
Data Availability
The raw/processed data/code required to reproduce the findings of this study cannot be shared at this time as the data also forms part of an ongoing study.
Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
We gratefully acknowledge the support provided by the National Natural Science Foundation of China (NSFC, Grant No. 51874288), the Jiangsu Natural Science Foundation (Grant No. BK20181356), and the Xuzhou Key Research and Development Projects (Grant No. KC18224).