The water burst of roof on working face has been one of the significant geotechnical engineering problems that needs to be urgently resolved. The coupling effects of seepage and damage on the amount and intensity of water inrush from the roof are critically important. In this paper, the seepage-damage coupling mathematical model of the aquifer in the working face is studied, and the seepage-damage coupling mechanics model at different stages of the aquifer is established. Under the coupling of permeability and damage, the water-soil characteristics of the aquifer in the 101163 working face of Mindong were numerically simulated by establishing the constitutive relation between vertical stress and permeability coefficient. The numerical results show that the stress concentration factor of the mining stress field gradually increases with the coal seam mining. The water-flowing fractured zone of the overburden is close to the communication of the quaternary aquifer. When the coal seam is excavated 250–300 m. Three free surfaces appear in the groundwater pressure field, and a large falling funnel is formed to establish a deep flow S-well well flow model. The research on the mining stress field and seepage field is carried out in combination with the Jakob formula. It is found that two sectors with reduced permeability of the fan surface are formed in front of the work. The variation law of the apocalyptic permeability infiltration under different mining distances, different coal seam thicknesses, different water pressures, and different roof management modes is studied systematically. The research indicates that the seepage flow under the condition of seepage infiltration of the lower aquifer should be between 50% and 100% of the traditional calculation method. The research results can help to deepen the understanding of the process of water inrush under the coupling of stress and seepage.

1. Introduction

Water inrush from roof is an important factor influencing coal mine safety production. Study on characteristics of rock seepage-stress coupling and on characteristics of advanced seepage reduction of working face in mining engineering has important theoretical and application value. A large number of Chinese experts and scholars have carried out a lot of experiments and theoretical researches in this regard [1, 2], mainly on permeability coefficient and permeability law of rock under mining stress coupling and numerical simulation research for the establishment of rock permeability evolution model and seepage-stress coupling model. Many experts and scholars have also established different rock seepage coupling models from the perspective of mathematics and mechanics [36]. The classical viscoelastic-plastic constitutive theory is used to treat the seepage field on the assumption that no damage is caused to the pore and fissure structure in the rock, which is also called the classic seepage-stress coupling model. Based on the classic seepage-stress coupling model, the more complex coupling effect caused by damage changes in the pore structure system of the rock is also studied, which is called seepage-stress damage coupling model. The research on permeability changes in coal mining in countries other than China was carried out earlier, mainly on theoretical analysis, field tests, and numerical simulations. Baghbanan et al. [7, 8] established a relation model between permeability and stress of overlying rock strata during coal mining based on the relationship between stress changes and permeability coefficients in mining strata and they calculated changes of permeability coefficient of overlying rock before and after mining. Fang et al. [9] used EDEM to simulate the influence of stress on permeability and the fluid flow pattern in fractured rock mass is studied when the pore size of fracture is related to the trace length of fracture. Ghabezloo et al. [10] found that the water flow from fissure decreased significantly with the increase of the fracture normal stress through laboratory experiments and came up with empirical formulas for the permeability coefficient, stress, and fissure parameters. Heiland et al. [11, 12] studied the relationship between rock porosity, permeability, and effective stress of different rock types and different pore types. Li et al. [1316] revealed the relationship between the permeability coefficient and stress of fractured rock mass.

In recent years, scholars have made great efforts in the research on the path and mechanism of mine inrush. Various theoretical and mathematical models have played a good role in coal mine safety production. Scholars proposed complex rock seepage-stress coupling models including fracture mechanics model and damage mechanics model. However, these models failed to simulate the characteristics of complex composite rock masses and the entire process of fissure and water inrush induced by coal mining. The instability of rock mass and abrupt change of seepage flow are caused by the evolution of stress-seepage coupling system, so the stress-seepage coupling model of rock mass needs to be further optimized and different permeability should be adopted for rocks in different damage stages and in different regions. In particular, the confined hydraulic structure characteristics and coupling mechanism of aquifers in different locations are not clear and traditional numerical simulation constitutive models and damage criteria are not fully justified. In this paper, the geological conditions of the coal-bearing sandstone aquifer in Mindong Mining Area were taken as the research object to explore scientific issues such as the establishment of a coupled mechanical evolution model of the aquifer above the coal seam on the working face at 101163 and the establishment of a numerical model under the coupling of stress and seepage, to expound the evolution of seepage characteristics of confined aquifers and migration law of groundwater flow field and conduct a comprehensive study of the advanced seepage reduction characteristics of the working face under the coupling of stress and seepage.

2. Analysis of Mechanism of Coupling between Mining Stress and Seepage

2.1. Analysis of Seepage Characteristics of Mining-Induced Fractured Rock Mass

With the development of fissures, working face can be divided into three areas from top to bottom, virgin rock area, fissure area, and caving area, as shown in Figure 1.

When aquifer I is in the original rock area, the mode of action of groundwater on rocks is disintegration and degradation. When aquifer II is in the mining-induced fissure area, the action mode of groundwater on rocks is shear action. When aquifer III is in the craving zone, the rock of the aquifer collapses as the mining face is exploited. As a result, the circulation channel of the aquifer is broken and the water in the aquifer in the craving area directly enters the working surface vertically to gushing water direction. The water in the aquifer on both sides enters the working surface in the form of lateral recharge and the advance bearing pressure leads to a change of permeability of aquifer and fissure shape.

2.2. Study on Influence of Mining Stress on Seepage Field

Figure 2(a) is a composite bearing layer composed of soft rock overburden and aquifer. Before the initial fracture, the composite bearing layer will be subject to the gravity of the overlying rock and the water pressure of the aquifer. The aquifer represented by the lower layer is composed of medium conglomerate. Therefore, the aquifer mainly contains solid conglomerate and internal cracks and pores. A microunit is selected at the interface between aquifer and aquifer. The fracture is affected by rock pressure of rock and water pressure above the rock. For virgin rock area, under the mining and water pressure, water inrush in coal seam is the development of rock stratum and the expansion of various structural panels at different angles from the floor. If the aquifer floor is not collapsed, we can still believe that the reason of water inrush lies in the expansion of various structural planes at different angles within the rock. It is most dangerous when joint or fissure is at 90° from the floor. Assuming that the cracks are elliptical, the adjacent cracks do not affect each other, and the local variation of material properties is ignored, the open cracks are subjected to internal water pressure and external stress, as shown in Figure 2.

Major semiaxis of the fissure is set as a and the minor semiaxis of the fissure is set as b. According to the mechanics of elasticity, the solution formula for the stress value at any point on the fissure is (1).where is the shear stress at point A, is the stress resultant on fissure, a is the eccentric angle, m is the number of strata from aquifer to surface, σx is the x direction normal stress, σy is direction normal stress, and τyx is shear stress.

The fissure unstable propagation meets the Griffith strength criterion, in order to simplify the calculation when considering that the bottom water-proof layer of the goaf is a uniform and horizontal cross-fixed plate. The entire water-proof layer is curved downward from the center of the circle, the roof is open fissure, and the bottom is fracturing fissure. The middle height of the plate can be considered as the contact surface where the fractures meet. The critical water pressure value is where is the critical water pressure, R is the tensile strength of the rock formation, h2 is the thickness of the effective water-proof layer, γ is the unit weight of the rock formation, Q is the pressure of the mine, Lx is the periodic weighting length, and L is the length of the working face.

During mining, fissure will appear in fissure area and then gradually develop. Groundwater in the confined aquifer enters fissures and produces seepage forces in the fractured rock mass, including hydrostatic pressure and hydrodynamic pressure. The hydrodynamic pressure is mainly resulted from the lateral recharge of the aquifer, causing groundwater to fill the fissure space and interact with the overlying rock layer continuously. As a result, the gravity of the overlying rock layer acts on the bending zone in the form of hydraulic pressure through the rock layer above the bending zone. The change stress of the rock mass will affect the structure of the rock mass, thus changing the permeability, which is manifested as the changes in the groundwater seepage field. The act of mining stress on the permeability coefficient can be obtained withwhere K is the permeability coefficient for the stress value, K0 is the permeability coefficient under zero stress, Δ is the change of this variable, σi is the consolidation pressure, E0 is the effective modulus of the particles, and σc is the confining pressure.

The relationship between the permeability coefficient and the internal particles can be obtained through the Hubbert relation expression as follows: where N is a dimensionless number related to the filling particles, ρ is the density, is the acceleration of gravity, µ is the coefficient of kinetic viscosity, and d is the average particle size.

The Bai. M’s method [17] is used to evaluate the influence of stress on the development of fissures. Assuming that the confining pressure and pore pressure result in effective stress, equation (5) can be expressed as the following formula with the cube model:where Rk is the particle radius. When the cubes are piled up, Nd2 grams are expressed as

Equation (5) can be expressed as equation (7) in case of insufficient compaction.

Under uniform pressure, equation (8) can be applied to obtain particle radius.where ν is Poisson’s ratio of the drainage of the medium. Assuming that K0 can be obtained through experiment, the particle radius can be expressed by

The expression of the influence of the mining stress on the seepage field can be obtained through simultaneous formulas (7)–(9) as shown in formula equation (10).

For the caving zone, the permeability of the overlying rock is related to its porosity, and there is a direct relationship between porosity and stress which has been extensively studied [1824]. Their relationship is where Ф is the porosity of the rock under stress, Ф0 is the porosity under the zero stress, Фr is the limit value of the porosity under the high-pressure stress, αФ is the stress influence coefficient, which can be set as 5.0 × 10−8 Pa−1, and συ is the average effective stress.where σ1, σ2, and σ3 are three main stresses, p is the water pressure in the aquifer, which is set as 3.5 MPa, and the relationship between the permeability and the porosity of the rock mass is shown in equations (13) and (14):where ξs is the mutation coefficient of the permeability in the plastic zone, which is set as 4.5–5.5, Ks is the permeability on plastic zone, K0 is the permeability under the zero stress, m2, is the coefficient of kinetic viscosity which can be set as 10−3 Pa·s, and is the unit weight of water, Pa/cm, which can be set as 100.

3. Numerical Simulation of Water Inrush in Working Face under the Coupling of Mining Stress Field and Seepage Field

3.1. Establishment of Constitutive Model
3.1.1. Governing Equation for Rock Mass

The effective stress in the saturated zone iswhere σ is the total stress and σ′ is the stress acting on the solid matrix.

Equation of solid displacement field is as follows:where G is the shear modulus, ▽ is gradient, αx is the Biot coefficient, p is the underground water pressure, and G can be derived from the following formula:where E is elasticity modulus and the Biot coefficient αx is mainly related to the compression performance of the material, which can be derived from the following formula:where is Poisson’s ratio of drainage of the solid component and B is the Skempton coefficient, which can be set as 0.4.

The equation of groundwater seepage field is as follows:where ▽p represents the influence of the interstitial pressure produced by the fluid flow on matrix stress field and ρ is the density of fluid.

3.1.2. Stress Correlation of Permeability

The traditional rock constitutive model under uniaxial stress [2528] is shown in Figure 3.

As shown in Figure 3, rock experiences elastic change from 0 to fc0 at a uniaxial stress state, but the rock is actually anisotropic material with internal friction and can yield and break under high strength pressure. After undergoing the transition from elasticity to plasticity, its permeability will also increase sharply. The relation expression of the impact of plastic failure on permeability corresponds to the constitutive relationship in the figure. The formula is as equation (14).

On the basis of previous results, a study was carried out on the relationship between horizontal permeability of sandstone and stress state. Considering that the pressure loading path has been drawn up according to the trend of confining pressure of coal seam in No. III aquifer in Mindong Mining Area before the experiment, in the numerical simulation study, only the relationship between permeability and vertical pressure needs to be considered, without considering the relationship between permeability and confining pressure. Therefore, the numerical simulation can be carried out according to the following formula.

In the first stage, when γH <  > 1.7 γH,

In the second stage, when  > 1.7 γH,where is vertical stress, MPa.

3.2. Damage Criteria

The criteria for the plastic damage of rock mass are the modified C-M criterion and the traditional tensile failure criterion.where I1 is the first invariant of effective stress, J2 is the second invariant of effective stress deviation, α1 and K1 are the experimental constants related only to the viscous force c of the internal friction angle ѱ, and ƒt is the tensile strength of the rock mass. Compared with the original criterion, the modified C-M criterion puts intermediate principal stresses (σ2) and pore pressure p into account in three-dimensional space; as a result of this, it overcomes the main shortcoming of the C-M criterion and has been widely used in the numerical analysis of rock mechanics and engineering at home and abroad. Under any stress conditions, the tensile damage criterion is a prioritized judging factor.

3.3. The Establishment of Numerical Model

There are too many rock layers with different thickness in the Mindong Mining Area. Therefore, comprehensive histogram is adopted to replace the rock data for a certain borehole in this numerical simulation. The overlying rock that does not affect the simulation results is combined. The parameters of rock stratum and rock mechanical properties are shown in Table 1.

Based on the histogram and the length of the working face, a geometric model with a length of 800 m and a width of 358 m was established. In a vertical direction, a model was established hierarchically based on the different properties of the rock, simulating the excavation of the coal seam of 300 m.

3.4. Analysis of Numerical Simulation Results
3.4.1. Analysis of Changes of Mining Stress Field

Figure 4 is a stress-displacement diagram at different working face lengths, from which, we can see that excavation of coal seams results in stress concentration around working face, mainly distributed on coal pillars 0–40 m away from the coal wall, where coal seam under maximum stress > roof > floor. As the length of the working face increases, the maximum stress near the coal pillar gradually increases. When the lengths of working face are 120 m, 160 m, and 200 m, the maximum stresses are 18.7 MPa, 23.2 MPa, and 25.1 MPa, respectively. The maximum value of the roof caving displacement also gradually increases, as the length of working face prolongs to 0.29 m, 0.37 m, and 0.48 m, respectively.

3.4.2. Analysis of Change Law of Groundwater Flow Field

Considering that the purpose of the study is to research the change of the flow rate when rock is damaged, a simulation chart of inrush flow of the working surface under different conditions is made and shown in Figure 5.

It can be found that, with the working face advanced, inrush water volume reduced. The changing pattern of the inrush water is similar when the width of working face is at 120 m, 160 m, and 200 m. As the working face width increases, the inrush water flow at the working face decreases and the margin of decreases reduced. It can be considered that, with the width of the working face increasing, the caving increases, but the stress peak value is greater, and the flow rate is reduced, indicating that the water inrush path of the working face has a cascaded structure. As the local fissure closes, the amount of water inflow decreases, as the thickness of the coal seam increases, the amount of inrush water at the working surface increases. Although the increase in the thickness of the working face results in an increase in the stress concentration factor, the roof development height increases and so does the cross section of the inrush channel, resulting in the fact that an increase in the thickness of the coal seam and water inrush on the working surface increases too. This is related to the geological conditions and distribution characteristics of the aquifer. With the increase in water pressure, there are complex changes in the inrush water on the working surface and the water inrush is greatest when the water pressure is 1.50 MPa. When the water pressure is 3.50 MPa, 2.50 MPa, and 0.50 MPa, the changing curves of inrush water will cross at a later stage. All of the above indicates that, under the coupling of water pressure and rock mechanics, different effects on different inrush water volumes occur at different stages, which is consistent with the results of our previous theoretical analysis. Different roof management methods also indicate that the filling mining not only ensures that inrush water volume decreases rapidly but also ensures that the volume of the inrush water is small. The initial inrush is mainly caused by the static water on the roof. With filling mining, we can ensure that inrush water volume is always at a low value and finally reduces to 0.

3.4.3. Study on the Change Law of Permeability in Aquifer

As the study focuses on the change of permeability of damaged rock, a research on the relationship between different mining distances and permeability for coal seam with thickness of 8 m is carried out, as shown in Figure 6.

As shown in the figure, as the coal seam is excavated to 100 m, the permeability of the aquifer above the working face increases to about 1.6 times to 2 times of the original. At this time, the aquifer rock above the working face collapsed and the original aquifer low channel was destroyed. The water in the aquifer above the working face will all flow into the working face. In the later stage, the confined water mainly enters the working face through lateral recharge and acts on it; when the excavation reaches 200 m, the area above the working face is still the red permeability increasing zone, but two small permeability mutation areas are formed in the front and back of the working face. As the working surface advances to 300 m, we found that there is a sector permeability decreasing zone about 0–50 m in the front and back of the working face, respectively. There is an arc permeation increasing zone outside the sector permeability decreasing zone which is due to the formation of the rock structure. Leading bearing pressures are formed in the front and back of the working face, and lateral bearing pressures are formed on both sides of the work face. Leading bearing pressure and lateral bearing pressure cause vertical pressure to the aquifer around the working surface. As the horizontal distance from the coal wall of the working surface varies, the bearing pressure is different too, and so does the lateral confining pressure. Therefore, lateral permeability coefficient reduces in the aquifer within the range of the leading bearing stress area. When the working face is at 400 m in depth, it can be found that the sector permeability reduction area has expanded to the entire aquifer section. As a result, the volume of dynamic recharge water in the aquifer is lower than that in traditional calculations. Considering that the permeability is proportional to the flow rate, the permeability is reduced by 40%–50% at a maximum if the change in head pressure is not taken into account, so the seepage flow should be between 50% and 100% of that in traditional calculation.

In order to further understand the change of the aquifer’s leading permeability reduction with different coal seam thicknesses, hydraulic pressures, and roof management methods, it has simulated the change of permeability of coal seam at the thicknesses of 8 m, 16 m, 24 m, and 32 m at the depth of 300 m. At the mining depth of 300 m, sector permeability reduction zones of different size are formed in aquifers in coal seam with different thicknesses. As the coal seam thickness increases, the area of sector permeability reduction zones decreases, which is different from the traditional view that the thicker the coal seam, the greater the leading bearing stress. From the simulation results of permeability changes under different water pressures, at the mining depth of 300 m, sector permeability reduction zones of different size are also formed in aquifers under different water pressures; we can see that when the water pressure increases, the area of sector permeability reduction zone gradually enlarges.

As shown in Figure 7, when the roof is managed with the direct caving method, sector permeability reduction zone is formed in the front and back of the working face. When the filling mining method is used, the gob is filled with a lag of 30 m, and the sector permeability reduction zone is not formed in the front and back of the working face. This fully verifies that the aquifer’s advanced antiseepage results from the increase in the leading bearing stress.

When the confined aquifer is within the falling zone, the aquifer above the working face collapses and the circulation channel of the aquifer is damaged. Moreover, all the water in the aquifer above the working face will flow into the working face. In the later stage, the confined water mainly enters and acts on the working face by lateral replenishment and with the working face continues to push forward. The aquifer behind the working under the ahead support pressure will form a fan permeability reduction zone, which corresponds to a curve surface infiltration increasing area outside the fan permeability reduction area. From the above analysis, it can be found that, as the mining distance, the thickness of the coal seam, the hydraulic pressure, and the roof management method vary, the sizes of the sector permeability reduction zone and the arc permeability increasing zone change. When the size of sector permeability reduction zone is larger than the height of the aquifer (since the lateral recharge of the aquifer into the working face is similar to tandem), the amount of lateral recharge water in the aquifer entering the working face reduces. Similarly, the lateral bearing pressure will also have the same effect on the coal seam side, so advanced antiseepage can be seen in the aquifer around the working face.

4. Conclusions

The results obtained made it possible to draw the following conclusions:(1)In this work, the seepage-damage coupling mathematical model of the working face roof aquifer is studied, and the seepage-damage coupling mechanical model of the aquifer in different stages is established. The variation law of the aquifer seepage leading to seepage reduction under different mining distance, different coal seam thickness, water pressure, and roof management mode is systematically studied.(2)In the upper part of the mining fracture area, the vertical tensile stress produces transverse cracks and fissures, forming a certain subsidence space. There are many vertical cracks in the middle of the mining fracture area, which is due to the bending of the rock layer, the vertical cracks are formed by the compression at the upper end of the middle, and the vertical cracks are formed by the tension at the lower end of the middle. Finally, many transverse and longitudinal fractures are formed in the lower part of the mining fracture area. By establishing the mechanical structure model of bending bearing stratum under hydrostatic pressure, the mechanism of hydrostatic pressure of confined water on aquifuge of bending subsidence strata is studied, the relationship between overburden fracture and water pressure and mine pressure is quantified, the shear force expression of aquifuge in bending subsidence zone is obtained, and the coupling relationship between mining stress and seepage field in mining fracture area is studied. The influence expression of mining stress on seepage field is given.(3)The numerical simulation study of the change of the mining stress field and seepage field was carried out. During the mining of coal seams, the stress concentration factor of the mining stress field gradually increased, and as a result, floor became larger than roof in stress concentration area. In this paper, leading antiseepage of the aquifer is discussed and the change of leading antiseepage of the aquifer with different mining distances, coal seam thickness, water pressure, and roof management methods is studied. The study finds that the seepage flow of the leading antiseepage of aquifer accounts for 50%–100% of that of the traditional flow.

Data Availability

This is an open-access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The authors acknowledge the financial support for this work provided by the Natural Science Research Program of Anhui Universities (KJ2020A0319), Open Foundation of State Key Laboratory of Mining Response and Disaster Prevention and Control in Deep Coal Mine (SKLMRDPC19KF10), National Key Research and Development Program (2019YFC1904304), and Shanxi Science and Technology Major Special List Project (20191101016).