Fracture evolution of a thick soft protection layer and the water inrush mechanism in overburden under longwall mining


Abstract

In the Jurassic coal-bearing strata of western China, thick soft protection layers (TSPLs) are widespread and play a key role in controlling roof water inrush hazards. Understanding their time-dependent deformation and failure mechanisms is essential for hazard prevention. This study combines laboratory creep tests, mechanical modeling, and numerical simulations to investigate the fracture behavior of TSPL under mining. A modified triaxial testing system was developed to characterize both shear and tensile creep. A viscoelastic thick-plate mechanical model and a numerical model were established to capture progressive deformation and failure processes. Results show that tensile creep leads to failure in the decelerate creep stage with limited deformation, whereas shear creep is dominated by large-scale rheological deformation and failure in the accelerate creep stage. The rheological properties and the thickness of TSPL jointly control fracture patterns: thicker layers tend to develop inclined shear-induced fractures connecting the high-position separation zone to the working face, while thinner layers are prone to vertical tension-induced fractures linking to the goaf. Further, four typical types of mine roof hydraulic hazards are identified, including aquifer water inrush via the water-conducting fractured zone, water accumulation in separation layers without inrush, water inrush from separation layers into the goaf, and water inrush from separation layers at the working face. Mechanism-based prevention measurements tailored to each water inrush type are proposed to enhance mine safety and ensure sustainable coal extraction.

Highlights


  • Thick soft rock layers in Jurassic coal mines show time-dependent deformation that controls water inrush risks.

  • A viscoelastic mechanical model is established to predict progressive deformation and failure of a thick soft protection layer.

  • Thicker protection layers tend to experience shear failure near the working face, while thinner layers develop tensile cracks above the goaf.

  • Four distinct water inrush patterns were identified, each requiring different prevention methods.


1 INTRODUCTION

With the strategic westward shift of China's energy production base, the Jurassic coalfields in the Shanxi–Shaanxi–Inner Mongolia junction area have become critical energy hubs. Recent mining activities in the Ordos Basin, particularly in key energy bases such as Shendong, Huanglong, and Ningdong, have revealed water inrush from separation layers, with inflows ranging from 300 to 2000 m3/h. These water inrushes typically show sudden, intermittent, and cyclic characteristics, frequently forcing workface shutdowns and potentially triggering major safety accidents, even resulting in casualties (Qiao et al., 2021). With thicknesses ranging from 20 to 220 m, the Jurassic coal-bearing strata in western China, deposited in fluvial and shallow lacustrine environments, are primarily composed of soft mudstones and weakly cemented sandstones (Li et al., 2024; Zhang & Mi, 2021). These low-strength rocks (uniaxial compressive strengths, UCS < 25 MPa) form a thick soft protection layer (TSPL) between the conventional water-conducting fractured zone and the overlying Cretaceous aquifer. The TSPL acts as the key aquiclude controlling roof water inrush during fully mechanized top-coal caving in ultra-thick seams (10–30 m). The interaction between mining-induced stress and the time-dependent fracture evolution of TSPL governs the hydraulic connectivity between the Cretaceous aquifer and coal seams, exerting a decisive influence on water hazard development (Qiao et al., 2023a). Similar hydro-mechanical instability mechanisms have been observed in the Bowen Basin of Australia and the Paleozoic basin in the United States, underscoring the global relevance of this issue (Rajabi et al., 2024; Rios et al., 2024).

Since the early 20th century, strata control theory has developed into an integrated framework of theoretical modeling, experimental simulation, and field validation. Foundational models, such as cantilever beam theory, pressure arch theory, articulated block, and masonry beam models, have emerged through coupling equilibrium conditions, constitutive laws, and failure criteria (Jiang et al., 1993; Liu et al., 2017; Peng, 1984; Song, 1998). Among them, key stratum theory (Qian et al., 1996) has dominated Chinese mining practice for decades and has since evolved into multiscale coupling models, including the composite key stratum, short masonry beam, fracture–stress arch, and quasi-hyperbolic models (Ju et al., 2013; Li, Yang, et al., 2020; Miao et al., 2011; Song et al., 2021; Sun et al., 2021; Wang et al., 2019).

Traditional key strata are typically hard rock layers with high strength (UCS > 35 MPa, Elastics modulus, E > 30 GPa), thin thickness (<30 m), and a transition mechanism from a continuous to a discontinuous structure during failure. This framework has been proven to be effective in eastern and northern coalfields with burial depths <500 m and hard sandstone roofs (Qian et al., 2003, 2010). However, the Jurassic strata in western China significantly deviate from these assumptions. In the Huanglong coalfield, for instance, the Jurassic Anding Formation consists of thick mudstones (20–80 m) with low strength (UCS < 25 MPa, E < 5 GPa), undergoing pronounced deformation and creep under mining-induced tensile and shear stress paths (Lü et al., 2020; Zhang et al., 2024). Such conditions fall beyond the scope of classical beam or thin-plate models.

More importantly, unlike traditional key strata that passively transfer stress, the TSPL in Jurassic strata actively participate in stress redistribution through long-term shear slippage and rheological deformation. These layers show dual mechanical behavior, as both load-bearing and deformable media, challenging the core assumptions of classical key stratum theory (Qian et al., 2010). Their progressive degradation governs both strata pressure evolution and roof water inrush formation, especially by forming large-scale separation spaces that enhance hydraulic connectivity with the overlying Cretaceous aquifer (Fan et al., 2022; Qiao et al., 2021).

Consequently, water inrush mechanisms in western China's Jurassic coalfields must be analyzed from the perspective of the holistic deformation–fracture evolution within the TSPL, rather than attributing them solely to the instantaneous failure of traditional key strata. Although prior studies have explored the mechanical and hydraulic behavior of aquicludes, their focus has predominantly remained on instantaneous or elastoplastic responses, with limited consideration of rheological effects. For instance, Fan et al. (2019) established a Winkler foundation beam model to predict initial and periodic water inrush in separated layers. Sun et al. (2020) integrated mineralogical and mechanical data to model hydro-mechanical coupling in aquicludes, while Tao et al. (2024) investigated stress–permeability responses in mudstones under varying load paths. Guo et al. (2023) developed damage evolution models to evaluate aquiclude stability under infilling.

In practice, field observations of water inrush from separation layers show that water inrush frequently lags behind both peak microseismic activity and the onset of groundwater-level anomalies (Qiao et al., 2011). This temporal offset suggests that the TSPL does not fail instantaneously, but instead undergoes a progressive degradation process under long-term mining-induced stress. Specifically, the TSPL gradually enters a creep-dominated damage stage, where slow, sustained deformation leads to interlayer delamination and the formation of water-bearing cavities. Even without obvious macroscopic deformation, the continued accumulation of micro-damage eventually reduces the rock mass strength to a critical threshold, resulting in sudden hydraulic breakthrough. This process accounts for the typical “sustained water level drop, followed by delayed inrush” behavior observed in Jurassic coalfields. Therefore, the high-level water inrush should be reinterpreted as a rheologically controlled phenomenon, governed by time-dependent damage evolution rather than brittle failure alone.

Despite these advancements, the lack of a rheological mechanics model capable of describing their time-dependent deformation and failure under mining-induced disturbances has limited the accurate prediction of water inrush hazards. Moreover, the role of fracture evolution in TSPL in controlling roof water inrush types and disaster mechanisms remains insufficiently understood. Consequently, this study develops a refined theoretical model to characterize the time-dependent failure behavior of TSPL, providing a scientific basis for hazard assessment and prevention in Jurassic coal mining operations.

This study is organized as follows: Section 2 presents the engineering geological characteristics of TSPL. Section 3 establishes a mechanical model of TSPL to examine the influence of time-dependent behavior and the thickness on fracture evolution. Section 4 utilizes numerical simulations to investigate the failure mechanisms of TSPL and their role in roof water inrush. Section 5 classifies water inrush types associated with the failure of TSPL and proposes targeted engineering mitigation strategies. Finally, Section 6 summarizes the key findings of this study.

2 MINING AND GEOLOGICAL CONDITIONS

2.1 Engineering geological conditions

The Cuimu colliery, located in the Huanglong coal base of southern Shaanxi Province, China, serves as a representative case for investigating the mechanical behavior of the TSPL, as shown in Figure 1a. The coal seam thickness varies significantly from 5 to 25 m, with extraction conducted using the longwall top-coal caving method. The operational mining height is set at 5–7 m at the working face, supplemented by 5–8 m of top-coal caving, resulting in a total extraction height of 10–15 m. Panel dimensions range from 180 to 200 m in width and 700 to 1200 m in length, with an advancing rate of 3–6 m/day. The mine uses a section pillar layout system, where protective coal pillars of 25–35 m width are designed between adjacent panels.

Details are in the caption following the image
Location and rock strata parameters of the Cuimu colliery, China. (a) Geological position of the study area, (b) stratigraphic sequence of the Huanglong coal base, and (c) the cross-section of Cuimu colliery.

The stratigraphic sequence of the Huanglong coal base, as shown in Figure 1b, consists of the following key formations: (1) Jurassic Yan'an Formation (J2y): The principal coal-bearing unit containing mudstone and fine sandstone, hosting the No. 3 coal seam (5–25 m thick) at depths of 550–700 m. (2) Jurassic Zhiluo Formation (J2z): Composed of interbedded mudstone and sandstone. (3) Jurassic Anding Formation (J2a): A sequence of mudstone–sandstone with an average thickness of 110 m and a low uniaxial compressive strength (<15 MPa). It functions as the primary aquiclude, separating the coal seam from the overlying aquifers. (4) Cretaceous Formation: A massive sandstone and glutenite aquifer (200–500 m thick) with pore water pressures ranging from 3 to 5 MPa, directly influencing the water inrush hazards.

As depicted in the cross-section in Figure 1c, the vertical distance between the Cretaceous Luohe Formation (K1l) and the coal seam at the Cuimu colliery varies from 135 to 225 m. According to standard empirical formulas (2017), the height of the water-conducting fractured zone is estimated to range from 73 to 87 m, with a mining height of 10–15 m. Therefore, the thickness of the protection layer above the fractured zone ranges from 48 to 152 m, classifying it as a typical TSPL. The deformation and failure of this TSPL are crucial in controlling the evolution of roof water inrush hazards.

2.2 Hydrogeological hazards and water inrush events

Taking panel 22303 of the Cuimu colliery as an example, as shown in Figure 2, mining-induced hydrogeological responses showed a distinct multistage pattern:
  • 1.

    Initial response stage: When the panel advanced to 100 m, the Cretaceous aquifer water level began to decline at a rate of 0.8–1.2 m/day. Despite this initial drop, no immediate water inrush occurred. During this period, microseismic signals remained low, characterized by sporadic low-energy events and low-frequency activity. This indicates minor fractures and localized stress adjustments within the lower strata, without the formation of significant water-conducting channels.

  • 2.

    Delayed sudden inrush stage: As the panel advanced to 300 m, the water level continued to decline steadily. Approximately 36 days after the initial drop, a sudden water inrush occurred in working face, with inflow surging from 30 to 200 m³/h within 2 h. By this time, the water level had dropped by approximately 37 m. In situ monitoring of Xie et al. (2024) in Cuimu colliery demonstrates that mining-induced separation layers primarily develop at the Jurassic–Cretaceous boundary and within the Cretaceous aquifer. Consequently, water inrush events show a lagged response following water-level decline, often displaying a periodic pattern. Just before the water inrush, microseismic energy significantly increased, with a peak shortly before the surge, indicating intensified fracture propagation within the TSPL and the formation of a water-conducting channel.

  • 3.

    Periodic inrush stage: As mining progressed, water inrush events showed a clear periodic pattern. At 500 m of panel advance, inflow increased to 120 m3/h, followed by another surge to 220 m3/h at 680 m. A third inrush event of 220 m3/h occurred at 1050 m. This indicates the periodic failure of the TSPL, ultimately leading to recurring water inrush events.

  • 4.

    Water level–inflow coupling oscillation mechanism: After each inrush event, the aquifer water level rebounded at a rate of approximately 1 m/day. Before the subsequent inrush, the water level gradually declined again, at a rate of about 1.2 m/day, usually occurring 5–20 days before the next surge. This recurring pattern of “decline–inrush–pause–rebound” underscores the dynamic interaction between aquifer water levels and periodic water inrush events. Each water inrush was preceded by a phase of heightened microseismic activity, with both energy and event frequency showing a stepwise increase.

Details are in the caption following the image
Water level, inflow, and microseismic activity of panel 22303 during mining.

Traditionally, high-level water inrush is attributed to the gradual subsidence and elastic–plastic failure of overlying strata, forming a through-going fracture channel. However, this study reveals that water inrush events typically lag behind the peak of microseismic activity and the onset of abnormal groundwater-level decline—phenomena that cannot be fully explained by a “rigid fracture” mechanism. By integrating microseismic monitoring data with groundwater-level responses, it is evident that the aquiclude composed of soft rock at the interface between the Cretaceous and Jurassic formations gradually enters a time-dependent rheological degradation stage under mining-induced stress.

As shown in Figure 3a, the persistent creep deformation of the soft rock layer leads to progressive delamination of the overlying strata, thereby forming a water-bearing separation layer. As time progresses, the strength of the rock mass continues to degrade due to rheological damage. Even in the absence of visible deformation, a sudden through-going fracture may occur once the long-term strength threshold is reached, establishing a hydraulic connection (Figure 3b). This progressive damage process accounts for the typical “sustained water level drop, followed by delayed inrush” behavior, demonstrating that water inrush from separation layers is fundamentally controlled by the time-dependent rheological damage evolution of the aquiclude, rather than conventional brittle slip or instantaneous failure mechanisms. Building upon this finding, this study establishes a rheological mechanical model for the TSPL to systematically investigate the spatiotemporal evolution of fractures under various engineering geological conditions, thereby revealing the characteristic fracture development patterns of TSPLs and their associated water inrush mechanisms.

Details are in the caption following the image
Evolution mechanism of water inrush from separation layers. (a) Initial response stage and (b) delayed sudden inrush stage.

3 ANALYSIS OF THE FRACTURE EVOLUTION IN TSPL BASED ON A THICK-PLATE THEORY

3.1 Experimental study of shear–tensile creep of rocks

3.1.1 Experimental equipment

To investigate the shear and tensile creep behavior of TSPL under unloading conditions, a modified triaxial creep testing system was developed based on the Multi-Field Coupling Triaxial Testing Apparatus (MCTTA) produced by France Top Industry (Liu et al., 2024). Since overburden failure in collieries is predominantly characterized by tensile and shear fractures, two experimental schemes were designed: (1) a shear creep test using two pairs of shear metal plates to generate two shear planes along the axial direction of the specimen and (2) an unloading tensile creep test, conducted using a specially designed unloading-tensile auxiliary device, as shown in the diagram in Figure 4.

Details are in the caption following the image
Auxiliary device for the triaxial tensile and shear creep test. (a) Multifield coupling triaxial testing apparatus, (b) shear creep test using two pairs of shear metal plates, (c) design drawing of the unloading-tensile auxiliary device, (d) physical diagram of the tension–shear testing apparatus, (e) sealing permeation holes using hot melt adhesive, and (f) the specimen after creep tensile failure.

The triaxial tensile test system relies on an auxiliary device that transforms axial pressure from the rock creep apparatus into axial tension while keeping confining pressure constant. The device consists of two-column and four-column tensile components, recessed connecting disks, and high-strength bolts. The two-column component has protruding guide rails, while the four-column component has matching grooves, allowing smooth vertical sliding with minimal lateral movement and friction (Huang et al., 2023).

The setup is mounted on the MCTTA's circular platform, with grooves ensuring stability. A cylindrical rock specimen is bonded between the tensile components, leaving a 2 mm gap to accommodate deformation. All parts are made of hardened steel, with SCM435 alloy bolts providing high tensile strength (≥80 MPa) and hardness (≥39 HRC). The device supports confining pressures of up to 120 MPa and fits MCTTA specifications (125 mm diameter, 103 mm height), accommodating specimens up to 25 mm in diameter and 50 mm in height. To prevent oil leakage, hot melt adhesive seals the permeation holes, maintaining stable confining pressure. The design ensures precise tensile loading under high-pressure conditions.

3.1.2 Experimental scheme

The creep experiments were performed on mudstone samples from the J2a layer of Cuimu colliery, adhering to the rock rheological test standards (Chinese Society for Rock Mechanics & Engineering, 2022). Preliminary tests established the peak tensile and shear strengths at 0.7 and 2.3 MPa, respectively, guiding the design of loading levels. A minimum of five stress levels were applied to ensure comprehensive creep stage coverage.

The tests used a stepwise axial stress increase under a constant 10 MPa confining pressure. The confining pressure was first applied at 0.1 MPa/s and held steady to simulate hydrostatic conditions. Axial stress was then incrementally increased at 0.5 MPa/s across multiple creep stages. Each creep stage maintained stress for 12 h to observe time-dependent deformation before proceeding to the next level. The process continued until the specimen entered accelerated creep, ultimately leading to failure.

3.1.3 Experimental results and analysis

The shear and tensile creep test results reveal significant differences (Figure 5). Under shear stress conditions, the specimens show typical three-stage creep characteristics: at τ = 1.0 MPa, instantaneous elastic strain and attenuation creep dominate; as stress increases to τ ≥ 3.3 MPa, strain accumulation becomes significantly enhanced; and at τ = 6.5 MPa, the specimen undergoes accelerated creep failure after approximately 1.5 h. Notably, due to the double-shear plane failure mechanism, the actual long-term shear strength is only half of the axial failure stress (approximately 3 MPa). In contrast, tensile creep demonstrates distinct behavior: even under high tensile stress (σt = −1.2 MPa), the specimen lacks a distinct acceleration stage but rather experiences sudden failure after 11 h.

Details are in the caption following the image
Creep curves under varying stress levels at a constant confining pressure of 10 MPa. (a) Shear creep curves and (b) tensile creep curves.

Existing studies indicate that rocks only show significant time-dependent deformation when the applied stress exceeds a specific creep threshold. Numerical modeling by Li (2019) under geological conditions similar to this study demonstrates that TSPL can experience vertical compressive stresses of up to 25 MPa and shear stresses fluctuating between 9 and 10 MPa under combined mining and water pressure effects. Meanwhile, in situ stress back-analysis by Ma (2020) reveals that static water pressure and separation space deformation can generate localized tensile stresses of approximately 5.5 MPa. Comparatively, the experimental results show that TSPL specimens enter the accelerated creep stage under tensile stresses as low as 1.2 MPa or shear stresses of 3 MPa, confirming that TSPL is indeed prone to creep rupture under mining-induced unloading conditions. These findings provide compelling evidence that roof water inrush involving TSPL should not be simply attributed to gradual movement or brittle collapse of overlying strata, but rather interpreted as a progressive rheological weakening process of soft rock under mining disturbances, ultimately leading to through-going failure.

To further analyze the viscoelastic response of the rocks during the shear and tensile creep test, the relaxation modulus was calculated, as shown in Figure 6. Based on the characteristic creep behaviors identified in both shear and tensile testing regimes, the Burger model was adopted for relaxation modulus fitting analysis (Sun, 1999).
(1)
where E( t) represents the relaxation modulus; E 1 corresponds to the instantaneous elastic modulus; E 2 denotes the delayed elastic modulus; η 1 is the viscous coefficient governing long-term creep; η 2 is the viscous coefficient controlling short-term creep; and t represents the creep time.
Details are in the caption following the image
Relaxation modulus of the rocks under the unloading test. (a) Shear creep test and (b) tensile creep test.

As shown in the fitting curves in Figure 6, this model effectively captures the key features of the experimental data (R² > 0.99), and the fitting parameters are listed in Table 1.

Table 1. Creep parameters fitted under unloading tension and shear at different stress levels.
Stress level E1 (MPa) E2 (MPa) η1 (MPa·s) η2 (MPa·s) R2
τ (MPa) 1.0 0.624 5.630 388.150 2.699 1.00
2.0 0.401 22.145 205.500 53.042 1.00
2.3 0.360 17.122 156.409 58.873 1.00
2.8 0.305 5.103 72.721 13.753 0.99
3.3 0.285 2.873 52.630 7.096 1.00
4.0 0.267 2.312 46.513 14.162 1.00
5.0 0.228 1.120 31.306 1.118 1.00
σt (MPa) −0.4 0.368 10.425 3.691 0.426 0.99
−0.8 0.646 87.263 2.801 39.485 1.00
−1.0 0.610 74.910 1.896 68.914 1.00
−1.2 0.586 83.992 0.388 39.838 0.99

3.2 Viscoelastic mechanical model of the TSPL

To investigate the deformation and failure behavior of TSPL during coal mining, a viscoelastic mechanical model based on the thick-plate theory is established, as illustrated in Figure 7. Given the significant creep effects of soft strata, their mechanical response must account for time-dependent deformation. The following fundamental assumptions are made to ensure model rationality and computational feasibility:
  • 1.

    The overlying thick soft rock layer is considered a homogeneous, isotropic continuous medium, with its deformation described using macroscopic mechanics, neglecting discrete fractures' influence on the initial stress–strain relationship.

  • 2.

    Although the soft rock layer may undergo large displacements over prolonged creep, the theoretical derivation assumes small deformations at the moment of disclosure.

  • 3.

    The time-dependent deformation of the soft rock layer follows a viscoelastic response, characterized by instantaneous elasticity, linear viscosity, and attenuated creep. The bending stiffness is treated as a time-dependent function to capture the creep-induced stiffness degradation.

  • 4.

    Creep primarily affects deflection, while bending moments and shear forces originate from the instantaneous stress distribution, consistent with the generalized Hooke's law in linear viscoelasticity.

Details are in the caption following the image
Thick soft protection layers (TSPL) model of four sides' fixed support.
Based on these assumptions, the viscoelastic thick-plate bending model was created based on Reissner's thick-plate theory (He & Shen, 1993; Reissner, 1944). Under a uniform load q, the dynamic governing equation is expressed as
(2)
where w is the deflection, D is the bending stiffness, , E is the elastic modulus, h is the thickness of TSPL, h =  H-H w, H is the distance from the coal seam roof to the Cretaceous basement, H w is the cumulative height of the fractured zone and the caving zone, q is the uniform load subjected to TSPL, and v is Poisson's ratio.
Since the goaf filling provides long-term support to the high-level strata, the uniform load is simplified to streamline calculations and facilitate parameter determination:
(3)
where k s represents the support coefficient of the gangue filling body, ranging from 0 to 1. For simplification, it is assumed to be a constant (Pan et al., 2022). q 0 accounts for the dead weight of the beam and the overlying strata load.
The laboratory results in Section 3.1 demonstrated that the time-dependent relaxation response under unloading conditions could be accurately captured by the Burgers model, thus supporting its applicability in modeling large-scale rheological behavior of TSPL. Accordingly, in the construction of the viscoelastic thick-plate model, the same constitutive form was adopted for the relaxation modulus, enabling a more faithful representation of the progressive deformation of TSPL under mining-induced stress disturbance. Therefore, the time-dependent bending stiffness D ( t) is given as
(4)
Incorporating this into the governing equation leads to the viscoelastic thick-plate deflection equation:
(5)
The deflection could be expressed as
(6)
Before the initial breaking, the main roof had been directly fixed in the surrounding rock strata. Therefore, the main roof boundary conditions prior to the initial breaking were four-sided fixed supports and a uniform pressure. Essentially, the boundary conditions could be written as follows:
where and are the rotation angles. a and b correspond to the plate length in the x-direction and the plate width in the y-direction, respectively. In engineering terms, a represents the advance distance of the working face in the strike direction, while b denotes the panel width in the dip direction.
Considering the boundary conditions, a double-sine series is chosen as the trial solution, naturally satisfying the fixed boundary displacement conditions:
(7)
where m and n are the spatial mode numbers of the deflection function; W mn ( t) represents the amplitude variation of the corresponding mode over time.
Substituting into the governing equation yields
(8)
Solving this equation provides the final deflection expression as follows:
(9)
To further simplify the model, only the primary terms are retained by setting m =  n = 1:
(10)
From Reissner's thick-plate theory (He & Shen, 1993), the shear force ( Q x, Q y) and bending moment ( M x, M y, M xy) are given by
(11)
(12)

Considering that creep effects primarily influence deformation rather than the distribution of bending moments and shear forces, the initial elastic stiffness D0 is retained for these calculations. While the effective modulus E (t) governs long-term deformation, the instantaneous stress state remains controlled by the initial modulus E0, ensuring that the stress–strain relationship at any given moment follows generalized Hooke's law. As a result, despite the time-dependent evolution of deflection, bending moments and shear forces are still computed based on D0, a common assumption in long-term structural mechanics analysis that maintains both computational rationality and physical consistency.

By solving the governing equation, the final expressions for deflection, bending moment, and shear force are obtained, which form the basis for further analysis of failure mechanisms in TSPL. By solving Equations ( 11) and ( 12), the following results are obtained:
(13)
(14)
The shear stress ( τ xz, τ yz, τ xy) and the normal stress of TSPL ( σ x, σ y) were calculated using the following:
(15)
(16)

For the normal stress, positive values represent compressive stress and negative values represent tensile stress.

3.3 Model validation and fracture evolution analysis

3.3.1 Mechanical parameter calibration based on in situ monitoring

To validate the proposed viscoplastic mechanical model, in situ monitoring data of strata displacement from Xie et al. (2024) were used to calibrate model parameters and assess its predictive accuracy. The monitoring borehole is positioned within the panel, at approximately 938 m from the cutting face along the panel strike (advance direction) and 7.8 m from the ventilation roadway along the panel dip direction. Given these conditions, the monitoring point was located at coordinates (938.0, 7.8). This thickness was determined based on a mining height of 15 m, with a calculated Hw = 87.46 m. Given H = 213.7 m, h = 213.7-87.56 = 126.24 m. The strata were continuously monitored starting from July 24, 2022. By October 30, 2022, the working face had advanced 260 m beyond the borehole, leading to an accumulated monitoring period of 100 days. Over this duration, the monitored deformation reached 0.045 m, corresponding to a total observation time of 3600 × 24 × 100 s. The working face spanned 200 m in width, and the full exposure length of the rock strata reached 1212 m, serving as the basis for model parameter calibration.

Due to the scale effect, rheological parameters obtained directly from laboratory-scale samples cannot be directly applied to field-scale rock mass models without correction. Therefore, a regression-based inversion approach was used, combining the mechanical model with in situ monitoring data to determine representative field-scale parameters. The accuracy and convergence of this inversion process are highly sensitive to the selection of initial values—large deviations between the initial estimates and the true values can lead to poor fitting results or even non-convergence. To address this issue, the average values of the Burgers model parameters fitted under different stress levels in triaxial shear creep tests were adopted as initial inputs for the inversion. This ensured a physically meaningful starting point and improved the stability of the optimization. As a result, the final calibrated parameters used in the mechanical model are as follows: E1 = 5.0 × 1012 Pa, E2 = 5.0 × 108 Pa, η1 = 6.5 × 1014 Pa·s, η2 = 5.0 × 1014 Pa·s, v = 0.3, ks = 0.85, q0 = 3.15 × 106 Pa, h = 126.14 m, a = 1212 m, and b = 200 m. A comparison between the theoretical predictions and monitored deformation results (as illustrated in Figure 8) demonstrates a high degree of agreement.

Details are in the caption following the image
Comparison of strata displacement from a mechanical model versus field monitoring data.

3.3.2 Stress distribution and potential fracture positions of the TSPL with time

Based on the model, the deflection and stress distributions of TSPL with dimensions a = 600 m and b = 200 m over a 1-month timescale (5, 10, 20, 30, 40 days) are shown in Figures 9-12. Due to the symmetry of the rock plate, only half of the plate in the strike direction (y = 0 to b/2 m) is displayed for better illustration of the stress distribution. As the neutral plane serves as the boundary where the stress transitions from compression above to tension below, the maximum tensile stress occurs at the bottom interface of the strata (z = −h/2). Thus, the distribution of σx and σy at the bottom interface (z = −h/2) is computed and plotted in Figure 10. On the other hand, the maximum shear stress occurs at the neutral plane (z = 0). Therefore, the distribution of shear stresses at the neutral plane (z = 0) is also calculated and shown in Figures 11 and 12.

Details are in the caption following the image
Deflection distribution and evolution of thick soft protection layers with time.
Details are in the caption following the image
Tensile stress distribution and evolution of thick soft protection layers with time. (a) σ x and (b) σ y.
Details are in the caption following the image
Horizontal shear stress τ xy distribution and evolution of thick soft protection layers with time.
Details are in the caption following the image
Vertical shear stress distribution and evolution of thick soft protection layers with time. (a) τ xz and (b) τ yz.

Over time, both deflection and stress initially increase rapidly, followed by a gradual increase after approximately 30 days, indicating that the deformation rate of the TSPL slows down and the degree of deformation tends to stabilize. This suggests that after 30 days, the development of separation spaces reaches a stable state, which reduces the rate of water accumulation within the separation layers. This phenomenon explains why, at Cuimu colliery, the water level in panel 22303 decreased rapidly over 36 days (indicating water filling the separation layers) and subsequently stabilized and even began to rise, as the deformation of the TSPL stabilized, leading to a slower or halted rate of water recharge.

The tensile stress distribution reveals that high-tensile-stress zones are primarily concentrated in the central goaf region (center of the plate, x = a/2, y = b/2). In contrast, σy is compressive in this condition. τxy along the strata interfaces shows significant concentrations near the open-off cutting (x = 0), the working face (x = a), and the strike roadway regions (y = 0 and y = b).

As proposed by Qian and Xu (1998), as the working face advances, early-formed separation fractures in the central goaf region are gradually compacted, while a continuous zone of bedding separation fractures develops around the periphery of the goaf, forming an “O-shaped” separation zone. The mechanical model presented in this study indicates that in TSPL, lateral shear stress concentrations along the rock edges also contribute to the formation of O-shaped interlayer separation fractures. However, whether these shear-induced O-shaped fractures evolve into an actual O-shaped separation space depends on the coordinated deformation behavior of the TSPL.

In the xz plane, τxz concentrations predominantly occur near the open-off cutting (x = 0) and the working face (x = a), with secondary peaks at the mid-dip position (y = b/2). In the yz plane, τyz concentrations primarily occur near the strike-oriented roadways (y = 0 and y = b) and the goaf, highlighting areas prone to shear-induced deformation and potential structural instability.

3.3.3 Water inrush prediction and model validation

Further validation of the model was conducted using field data from panel 22303, where a significant water inrush event was recorded. As the working face advanced to approximately 300 m, water inflow was observed. Based on rock mechanics tests, the long-term tensile strength and long-term shear strength of the rocks from TSPL under an initial confining pressure of 10 MPa were determined to be approximately 1.2 and 3.0 MPa, respectively.

Using the calibrated mechanical parameters in Section 3.3.1, the viscoelastic mechanical model was applied to calculate the failure length of the TSPL. Given that H = 185 m and Hw = 73.25 m, with a mining height of 10 m, the TSPL thickness was determined to be h = 185−73.25 = 111.75 m. The maximum stress locations were identified based on the analysis in Section 3.3.2. Subsequently, the stress evolution within the TSPL of panel 22303 was computed as the working face advanced, as shown in Figure 13.

Details are in the caption following the image
Stress evolution and prediction of the water inrush location of panel 22303.

The results indicate that τyz exceeded the shear strength threshold at approximately 290 m, leading to initial failure within the TSPL. In contrast, σx reached the tensile strength threshold later, at approximately 381 m. Notably, in the xz plane, shear stress τxz initially increased and then decreased, but remained below the shear strength threshold. Similarly, σy showed an increasing and then declining trend without exceeding the tensile strength limit.

Thus, the first water inrush event, observed at approximately 300 m, aligns well with the shear failure threshold in the yz direction, confirming that the initial inrush was triggered by shear fracturing of the thick soft rock layer rather than tensile failure. This consistency further demonstrates the reliability of the proposed mechanical model in accurately capturing the progressive failure mechanisms of TSPL under mining-induced stress.

Based on the evolution of tensile and shear stress with face advancement (Figure 13), the shear fracture mechanism of TSPL can be divided into two distinct stages:

Stage I (a < b): When the working face advancement distance is less than the panel length, τxz dominates over τyz. If shear failure occurs in this stage, it primarily develops along the strike direction, following the ventilation and transportation roadways of the panel. However, if the τxz does not exceed the shear strength before reaching the critical advancement distance (i.e., the “first square” position), xz-plane shear failure will not occur afterward, regardless of further face advancement.

Stage II (a > b): Once the advancement exceeds the panel length, shear failure transitions to being dominated by τyz. Beyond the critical position, shear fractures primarily develop in the dip direction, aligning with the cutting eye and working face. This finding aligns well with the O–X fracture pattern (Qian et al., 2003).

3.3.4 Influence of thickness on fracture evolution of the TSPL

The evolution of normal stress in rock strata with increasing advance distance, considering different thicknesses (40–140 m), is illustrated in Figure 14. From the σx graph, it is evident that in the early mining stage, σx remains positive (compressive) within a certain range. However, as the advance distance increases, σx transitions from compression to tension, marking a shift in the stress state that influences tensile fractures' initiation. The critical advance distance at which this transition occurs increases with rock thickness, indicating that in a thicker protection layer, a greater mining advancement distance is required before tensile failure occurs along the panel strike direction. This suggests that in thicker strata, the risk of tensile failure along the strike direction is lower at shorter advance distances, requiring a larger span before tensile stress becomes dominant.

Details are in the caption following the image
Evolution of normal stress with advance distance for different thick soft protection layers thicknesses. (a) σ x and (b) σ y.

In contrast, Figure 14b shows that σy remains tensile throughout the entire advance distance range, indicating that tensile stress along the dip direction governs fracture initiation in the early mining stage. Moreover, in the thinner protection layer, the magnitude of σy is greater than σx, reinforcing the idea that tensile failure along the dip direction is the primary mechanism in TSPL. As the thickness increases, σy shows a more pronounced decay, and for a protection layer with a thickness of 120 m, it transitions into compressive stress within the given mining range. This shift suggests that in an ultra-thick protection layer, the influence of tensile failure along the dip direction diminishes, while tensile stress along the strike direction continues to dominate even at a larger advance distance. Notably, for a 160 m thick protection layer, σx stabilizes at a negative value, indicating that in a very thick protection layer, tensile failure is primarily governed by stress along the strike direction.

Figure 15 illustrates the evolution of interlayer shear stress (τxy) with increasing advance distance for different protection layer thicknesses (h = 40–160 m). For thin to moderately thick strata (h < 80 m), shear stress decreases significantly with increasing thickness, indicating that thicker strata in this range are less prone to interlayer shear failure and tend to fail in a more cohesive, structural manner. In contrast, for a thicker to an ultra-thick protection layer (h > 100 m), shear stress increases with thickness, suggesting that greater thickness enhances the likelihood of interlayer shear failure, leading to separation within the TSPL.

Details are in the caption following the image
Evolution of τ xy with advance distance for different thick soft protection layers thicknesses.

Figure 16 illustrates the evolution of τxz and τyz with advance distance for different TSPL thicknesses (h = 40–160 m). As the thickness increases, both τxz and τyz shear stresses show a monotonic decrease, indicating that shear stress magnitudes diminish with increasing thickness. However, their evolution with advance distance differs: τxz first increases when the advance distance is shorter than the panel width, and it decreases once the advance distance exceeds the panel width. On the other hand, τyz continuously increases with advance distance. This suggests a shift in the control of shear stress during mining—from τxz dominance in the early stages to τyz dominance in the later stages. As mining progresses, shear failure becomes more likely along the coal wall of the working face, increasing the risk of large-scale shear collapse.

Details are in the caption following the image
Evolution of shear stress with advance distance for different thick soft protection layers thicknesses. (a) τ xz and (b) τ yz.

From a stress perspective, thin-to-moderately thick protection layers (h < 80 m) are more prone to failure due to tensile stress, with tension concentrated in the goaf, leading to tensile-induced fracturing. In contrast, for thick to ultra-thick protection layers (h > 100 m), tensile stress approaches zero or transitions into compressive stress, while shear stress decreases at a slower rate. This shift in stress distribution makes shear failure the dominant mode, with stress concentrations near the working face and roadways. Consequently, mining-induced failure transitions from tensile-dominated fracturing in the goaf for a thinner protection layer to shear-dominated collapse near fixed boundaries for a thicker protection layer, highlighting the increasing risk of shear failure and water inrush occurrence at the working face.

The mechanical model further reveals the deformation and failure mechanisms of the TSPL under varying thickness conditions. Under mining-induced stress, the rock mass is typically subjected to both tensile and shear stresses, and the failure mode depends on the interaction and evolution of these stress components. If tensile failure occurs before significant shear deformation develops, through-going fractures may rapidly form, leading to traditional water inrush centered in the goaf. In contrast, if the rock mass first undergoes shear-dominated large deformation, resulting in pronounced delamination and water accumulation, a subsequent tensile or shear failure may trigger the water inrush from separation layers. When tensile effects dominate, the inrush tends to initiate near the goaf area; when shear dominates, it is more likely to occur at the working face. Notably, if the TSPL is sufficiently thick and stiff, stress redistribution can dissipate energy progressively, preventing the critical evolution of tensile or shear failure. In such cases, although a delamination space may form and water-level changes occur, no through-going fractures develop and water inrush does not take place.

4 NUMERICAL SIMULATIONS AND ANALYSIS OF THE FRACTURE EVOLUTION OF THE TSPL

4.1 Model construction and parameter settings

A numerical model was developed using PFC to investigate the fracture evolution of the TSPL during mining. The Cretaceous strata in the study area primarily consist of thick coarse sandstone and conglomerate, with an average thickness of 260 m. As shown in Figure 17, the model incorporated a portion of the Cretaceous formations with a thickness of 105.40 m. The remaining overlying strata and the aquifer water pressure of 3 MPa were simplified as an equivalent uniform vertical load acting on the top surface of the Anding Formation, thereby simulating both in situ stress and water pressure effects. J2a has an average thickness of 126.1 m. Below this, J2z and J2y, with a combined thickness of 78.50 m, are characterized by thin sand–mud interbeds less than 5 m thick. The mining thickness in the model is set at 15 m.

Details are in the caption following the image
Schematic diagram of the numerical model. (a) Scheme 1 and (b) Scheme 2.

The model dimensions were 600 m (length) × 300 m (height), and excavation proceeded in 50 m increments. Each excavation step was followed by 10 000 computational cycles to simulate stress redistribution and fracturing processes. Particle diameters ranged from 0.8 to 1.5 mm, ensuring sufficient resolution to track microcrack initiation and propagation. The upper boundary of the model was left unconstrained (free to displace vertically) to simulate the natural subsidence of the overlying strata, with a uniform pressure of 5 MPa applied to simulate stress conditions. Zero-displacement constraints were applied to the left, right, and lower boundaries. Excavation was simulated incrementally, with 50 m protective coal pillars on both sides of the excavation zone to account for mining-induced stress evolution.

Creep deformation was realized by integrating the linear parallel bond model with the Burger model (PBM-BM) to simulate both instantaneous and time-dependent behavior of sedimentary rocks. The particle-scale constitutive assignment followed geological reality: particles with diameters <1.2 mm (representing clay minerals like illite and montmorillonite) were assigned pure Burger model parameters to capture viscoelastic behavior, while larger particles (>1.2 mm, representing quartz/feldspar grains) were assigned linear parallel bond parameters. As shown in Figure 18a, at lower stress levels, both PBM and Burger contacts jointly contributed to representing elastic and viscoelastic deformation. As stress increased and PBM contacts fractured, deformation transitioned into the viscoplastic stage governed by the Burger elements. When stress exceeded the failure threshold, contact bonds failed and residual strength was provided by interparticle slip resistance. Although the input viscosity coefficients varied linearly with time, the emergent macroscopic viscosity behavior was nonlinear, thereby capturing the creep characteristics of weak rocks effectively.

Details are in the caption following the image
Constitutive model and calibration of mechanical parameters for mumerical simulation. (a) Isochronous curves of the constitutive model, (b) calibration process of mechanical parameters, (c) calibration based on strata deformation, and (d) calibration based on fractured zone.

To ensure the reliability of simulation results, mechanical parameters were calibrated using a dual-validation approach addressing both deformation patterns and strata failure characteristics, as illustrated in Figure 18b. Specifically, subsidence curves of the J2a strata during mining were obtained from on-site monitoring by Xie et al. (2024) at the same coal mine. Corresponding displacement monitoring points were set in the model at the same elevations, and simulation parameters were iteratively adjusted until the simulated subsidence curves matched the field-measured results, achieving calibration of deformation behavior. In terms of fracture evolution, the simulated fracture zone height was compared to the measured fracturing height ratio (Hw/M) of 19.45 reported in the study of Qiao et al. (2023b), with the model result being 18.87, indicating good agreement and validating the damage calibration. The Burger model parameters were defined with stiffness moduli (Kelvin section: Normal stiffness modulus Kk = 5 GPa, shear stiffness modulus Gk = 5 GPa; normal viscosity coefficient ηkn = 0.5 GPa·s, shear viscosity coefficient ηks = 0.5 GPa·s. Maxwell section: Normal stiffness modulus Km = 7.9 GPa, shear stiffness modulus Gm = 7.5 GPa; normal viscosity coefficient ηmn = 3 GPa·s, and shear viscosity coefficient ηms = 3 GPa·s.) Other mechanical properties are provided in Table 1.

4.2 Simulation scheme

Since the TSPL is primarily located within J 2a, this formation serves as the key focus of the numerical simulation. Based on the occurrence characteristics of J 2a in the study area and the theoretical findings from the mechanical model, two simulation schemes (Figure 17) were designed to investigate the influence of TSPL on roof water inrush:
  • Scheme 1: The upper section of the Anding Formation consists of a continuous 80 m thick mudstone layer.

  • Scheme 2: To ensure that the TSPL experiences the same mining-induced stress conditions as in Scheme 1, the total thickness of the Anding Formation and its distance from the coal seam remain unchanged. However, in this scheme, the upper section contains only a 30 m thick mudstone layer, while the remaining portion consists of interbedded sand–mud layers.

The primary difference between the two schemes lies in the internal lithological structure of the TSPL, specifically the thickness of individual strata within the formation. By varying the thickness of these layers while maintaining the overall TSPL thickness, this simulation aims to validate the theoretical conclusions drawn from the mechanical model and elucidate the mechanisms of roof water inrush under different TSPL geological conditions.

4.3 Analysis of the simulation results

The simulated fractures, vertical displacement (dv), and contact forces for both schemes at their respective advance distances are illustrated in Figure 19. At 200 m of face advance, fractures develop within the near-field J2z and J2y above the goaf, forming a fracture network that periodically propagates with mining progression (Figure 19a). As the working face advances to 300 m, a pronounced displacement difference emerges between J2a and the overlying Cretaceous strata, accompanied by the formation of inclined shear fractures within the thick mudstone layer of TSPL (Figure 19b). At 350 m of face advance, inclined shear fractures within the thick mudstone layer of TSPL progressively extend and eventually coalesce, forming a water-conducting channel with an inclination of approximately 60°, which connects the separation water to the working face (Figure 19c). As the working face advances to 500 m, a new, higher-positioned separation layer forms within the Cretaceous strata, while periodic fracturing of the TSPL results in cyclic connectivity between the high-position separation water and the goaf, leading to periodic separation-layer water inrush events (Figure 19d).

Details are in the caption following the image
Fractures, vertical displacement, and contact force of the overburden during mining. (a) Advance distance of 200 m of Scheme 1, (b) advance distance of 300 m of Scheme 1, (c) advance distance of 350 m of Scheme 1, (d) advance distance of 500 m of Scheme 1, (e) advance distance of 300 m of Scheme 2, (f) advance distance of 400 m of Scheme 2, and (g) advance distance of 500 m of Scheme 2.
Details are in the caption following the image
Figure 19 (continued)
Fractures, vertical displacement, and contact force of the overburden during mining. (a) Advance distance of 200 m of Scheme 1, (b) advance distance of 300 m of Scheme 1, (c) advance distance of 350 m of Scheme 1, (d) advance distance of 500 m of Scheme 1, (e) advance distance of 300 m of Scheme 2, (f) advance distance of 400 m of Scheme 2, and (g) advance distance of 500 m of Scheme 2.

For Scheme 2, before the working face advances 300 m, the fracture patterns in near-field J2y and J2z remain largely consistent with those in Scheme 1. As mining progresses, a noticeable displacement difference develops between the TSPL and the overlying Cretaceous strata, leading to the formation of a separation space that accumulates water. However, after 350 m of advancement, shown in Figure 19e, as J2a undergoes significant deformation, the fracture evolution within the TSPL begins to diverge. As shown in Figure 19f, at 400 m of face advance, unlike Scheme 1, where inclined shear fractures develop within the TSPL, Scheme 2 shows no such shear fractures. Instead, failure is primarily governed by tensile fracturing, forming pathways that connect the high-position separation space to the goaf.

The fracture patterns in Scheme 1 and Scheme 2 arise from differing force transfer mechanisms in the thick, soft, and permeable layers (TSPLs). In Scheme 1, where the TSPL comprises thick sublayers (80 m), force chain analysis reveals a well-developed pressure arch within the mudstone layer (Figure 19b). This arch redistributes mining-induced stress, concentrating shear forces along its inner boundary. The resulting shear stress triggers inclined fractures that propagate toward the working face, eventually connecting with overburden fractures and enabling water inrush.

In Scheme 2, where the TSPL consists of thinner sublayers (30 m), the pressure arch is less stable (Figure 19e). The interbedded structure disperses stress, preventing localized shear concentration. Instead, vertical tensile stress dominates the TSPL's mid-span, promoting continuous tensile fractures. These fractures bridge the separation space above with the goaf below, forming a direct water inrush pathway. The contrasting failure modes highlight how sublayer thickness influences fracture development and water intrusion mechanisms.

In summary, these findings indicate that the TSPL in Scheme 2, composed of thinner sublayers (30 m), primarily develops continuous tensile fractures at the mid-span of the rock layers, facilitating water inrush into the goaf. In contrast, the TSPL in Scheme 1, with thicker sublayers (80 m), predominantly forms inclined shear fractures above the working face, leading to water inrush directly at the face. Furthermore, the shear-induced collapse of overlying strata above the working face increases the risk of severe mining pressure hazards, such as roof falls and support crushing. These conclusions align with the stress evolution and distribution patterns predicted by the thick-plate mechanical model, further validating its applicability in assessing overburden failure mechanisms.

5 DISCUSSION

5.1 Classification criteria and identification of water inrush patterns

Based on the viscoelastic mechanical model established in this study, the deformation w (t) and mining-induced normal stress σ (t) and shear stress τ (t) of TSPL show significant time effects and thickness dependence. Over time and across different locations, the failure mode transitions from tension-dominated to shear-dominated, leading to multi-stage, multi-mechanism coupling in roof water inrush types during coal mining. To systematically characterize the temporal evolution mechanisms of water inrush under different conditions, we define the following critical time thresholds, as shown in Figure 20: tw is the mining time at which TSPL deformation w (t) leads to the formation of a hazardous water-bearing separation layer wc; tsh is the mining time at which shear stress in the xy plane τxy (t) reaches the interface shear strength [τxy]; tt is the mining time at which tensile stress σx/y (t) reaches tensile strength [σt]; and tsv is the mining time at which shear stress in the yz/xz planes τyz/xz (t) reaches shear strength [τxz/yz].

Details are in the caption following the image
Critical time thresholds of thick soft protection layers with mining time.
Based on the viscoelastic mechanical model established in this study, a process for identifying roof water hazard types is proposed, as shown in Figure 21. First, the range of the protection layer can be determined based on the height of the fractured zone. According to the definition of thick-plate geometric conditions, if the ratio of the protection layer thickness to the working face width is greater than 1/10, it is classified as a TSPL. Furthermore, the sequence in which these characteristic stages occur during the advancement of the working face can be determined, allowing for the identification of water inrush types.
  • 1.

    Aquifer water inrush via a water-conducting fractured zone

    In this scenario, the TSPL undergoes rapid tensile or shear failure with minimal overall deformation, preventing the formation of large-scale enclosed separation spaces. Consequently, vertical mining-induced fractures directly propagate into the overlying aquifer, resulting in water inflow. Essentially, the TSPL functions as part of the water-conducting fractured zone, establishing a direct hydraulic connection to the aquifer (Hebblewhite, 2020). This process leads to a gradual and continuous water inflow, being similar to traditional inrushes with low-level aquifer and thin aquicludes (He et al., 2021; Li, Yang, et al., 2020). A typical example is the aquifer water inrush observed in the Cretaceous strata of the Gaojiabao Coal Mine in the Binchang mining area (He et al., 2025). The mechanical conditions are as follows:

    (17)


    The severity of this type of water inrush primarily depends on the aquifer's water abundance and recharge capacity. Effective prevention strategies include controlling mining height and optimizing roof support techniques to reduce the formation of large-scale water-conducting fractures. In areas with high-permeability zones, additional grouting reinforcement may be required to enhance roof stability and minimize hydraulic connectivity, thereby mitigating water inrush risks.

  • 2.

    Water accumulation in separation layers without Inrush

    Throughout the entire mining process, the TSPL experiences large-scale bending deformation, leading to the expansion of separation spaces and the accumulation of water. However, as the TSPL does not undergo tensile or shear failure, no water-conducting pathway is formed, resulting in water remaining confined within the separation layer. Under stable mining conditions with a consistent advance rate, inrush events do not occur. However, during periods of low advance rates or production halts, the long-term rheological deterioration of the TSPL may eventually generate water-conducting pathways, leading to delayed inrush incidents. A notable case is the Zhaoxian Coal Mine, where abnormal declines in the Cretaceous aquifer water level were observed in panels 1307, 1305, and 1303 without subsequent inrush events. However, panel 1304, after a 1-month production halt, experienced a sudden water inrush upon mining resumption (Qiao et al., 2021). The mechanical conditions are as follows:

    (18)


    Since the separation space has already accumulated water but remains structurally intact during normal mining operations, direct drainage can be achieved through separation drainage boreholes to prevent excessive water accumulation. Additionally, it is critical to avoid mining delays or prolonged stoppages, ensuring a steady and rapid advancement of the working face.

  • 3.

    Water inrush from separation layers in the goaf

    In this type, the TSPL undergoes significant bending deformation before failure, resulting in large-scale separation spaces and extensive water accumulation. As mining advances, tensile failure occurs at the mid-span of the TSPL, leading to the development of vertical fractures above the goaf. These fractures directly connect the accumulated water in the separation layer to the goaf, triggering instantaneous and high-intensity water inrush. Such hazards are prevalent in areas where the protective layer is relatively thin or the TSPL contains multiple weak interlayers. In these cases, τxy induces interlayer slip, reducing the thickness of separated sub-layers within the TSPL and promoting tensile failure. A typical example is Jining No. 2 Coal Mine, where large-scale water inflow occurred behind the goaf, exerting minimal impact on face mining activities (Qiao et al., 2011). The mechanical conditions are as follows:

    (19)


    As the water mainly seeps into the goaf behind the working face, it generally has a limited direct impact on ongoing mining activities. However, in cases of down-dip mining, accumulated water in the goaf may migrate toward the working face, posing a risk of flooding. Under such conditions, flow control measures should be implemented to prevent uncontrolled water migration, and in some cases, drainage boreholes in the goaf may be required to redirect water away from active mining areas.

  • 4.

    Water Inrush from separation layers in the working face

    Similar to the previous type, this scenario also involves large-scale bending deformation of the TSPL, leading to significant water accumulation in the separation layer. However, as mining progresses, shear failure occurs directly above the working face, forming inclined fractures that establish a direct hydraulic connection between the accumulated water and the working face. This results in sudden and intense water inrush at the working face. These hazards are more prevalent in regions with relatively thick and lithologically uniform protective layers. In such cases, τxy does not induce significant interlayer slip, and the TSPL maintains a continuous rock mass failure mode dominated by shear failure. Typical cases include the Cuimu and Guojiahe collieries panel 1309, where inrush events occurred directly at the working face, accompanied by strong mining-induced stress, significantly impacting mining operations (Jin et al., 2023; Qiao et al. 2023a). The mechanical conditions are as follows:

    (20)


    This type of water inrush represents the most complex and hazardous water inrush mechanism observed in typical coal mines of the Huanglong coal base. When the roof consists of weakly cemented soft rock, this process can lead to catastrophic water–sand inrush events, directly threatening mining operations. Therefore, a surface drainage system should be implemented, including pre-extraction boreholes to lower the water level and reduce recharge to the separation space. Furthermore, the fracture evolution and failure positions of TSPL should be predicted using mechanical modeling to improve hazard prediction and response capabilities.

Details are in the caption following the image
Flow diagram for determination of water inrush types.

5.2 Model assumptions and limitations

In the established viscoelastic mechanical model, the goaf support coefficient ks is treated as a constant, representing the effective stiffness or support provided by the collapsed strata within the goaf to the overlying TSPL. This simplification allows for a more tractable analytical framework; however, it does not fully reflect the complex spatial-temporal variations of ks under real mining conditions.

In practice, ks is not a fixed material constant, but rather a lumped parameter that encapsulates multiple factors, including the vertical distance between the aquiclude and coal seam, mining height, strength, and lithology of near-field strata, and their structural integrity. As such, ks inherently reflects the evolving geomechanical and hydrogeological environment of the stope, and may vary significantly both along the working face and over time as mining progresses. Particularly, the time-dependent evolution of the goaf compaction state and fractured zone closure implies that ks may show rheological or time-delay characteristics, introducing a spatiotemporal effect into the mechanical behavior of the protective layer.

Field observations support this interpretation. For example, as documented in Section 2.2, the 22303 working face showed periodic water inrush events during mining, which likely correspond to fluctuations in the effective goaf support and stress redistribution. These periodic variations in hydraulic response strongly suggest that ks is not constant, but dynamically evolves with the advancement of the face and the cumulative deformation of surrounding rock masses. Moreover, the spatial heterogeneity of ks may partially explain the coexistence of different water inrush types, such as separation-layer inrush near the coal wall and aquifer-derived inrush in the goaf, in different sections of the same working face.

Therefore, while the current model adopts a constant ks for the purpose of theoretical clarity and numerical stability, this assumption limits the model's capacity to capture the full spectrum of dynamic water inrush phenomena observed in the field. Future work should aim to develop a variable ks(x,t) formulation that incorporates both spatial geological variability and temporal mechanical evolution. This could be achieved through high-resolution field monitoring, coupled hydro-mechanical simulations, and inversion-based calibration techniques. Such refinement would enable the model to more accurately simulate and predict the periodic or hybrid water inrush behaviors often encountered in thick soft aquiclude mining scenarios.

6 CONCLUSIONS

In this investigation, laboratory experiments, mechanical modeling, and numerical simulations were used to examine the fracturing characteristics of thick hard rock strata in detail. Using a specially developed auxiliary device, tensile and shear creep tests of TSPL were conducted and a mechanical model of TSPL fracture was established based on the thick-plate theory. The effects of bending stress and shear stress on fracture types were analyzed, and the influence of different TSPL fracture modes on water inrush types was discussed. The main conclusions are as follows:

Field water inrush monitoring in the study area indicates that initially, the water levels in the Cretaceous aquifer decline without immediate inrush due to significant deformation but delayed fracture of the TSPL. As mining advances, water-conducting pathways form, triggering sudden inflows. The coupled oscillatory mechanism, characterized by a cycle of water-level decline, abrupt inrush at the working face, cessation of inrush, and subsequent water-level rebound, was clearly observed.

The Burgers model successfully describes the viscoelastic behavior of these rocks under constant confining pressure and axial shear/tensile stress conditions. Stress distribution analysis based on the viscoelastic mechanical model reveals that tensile stress concentrates in the central goaf region, while shear stress accumulates near the open-off cutting, the working face, and along the strike roadways. Model calculations indicate that shear failure in the TSPL triggers the initial water inrush at approximately 300 m of panel advance in panel 22303. The model further divides the shear fracture process into two stages, highlighting a transition from failure along the strike direction to failure along the dip direction. Finally, the influence of rock thickness on fracture evolution is critical, with thinner strata predominantly showing tensile failure and thicker strata showing a higher prevalence of shear failure near fixed boundaries.

A PFC numerical model was developed for the TSPL of Cuimu Colliery to analyze the impact of different TSPL structures on water inrush mechanisms. The results indicate that in TSPL with 80 m thick sublayers, a stable pressure arch forms, leading to shear fractures that propagate from the separation layer toward the working face. These fractures ultimately connect with mining-induced fractures, triggering water inrush directly at the face. In contrast, TSPL with 30 m thin sublayers lacks a stable pressure arch, resulting in continuous tensile fractures at the mid-span. These fractures establish a direct connection between the high-position separation space and the goaf, facilitating water inrush into the goaf. The numerical simulation results align closely with theoretical calculations and in situ monitoring data, further validating the findings.

This study categorizes water inrush hazards into four distinct types based on the time-dependent failure mechanisms of TSPL: (1) aquifer water inrush via the water-conducting fractured zone, (2) water accumulation in separation layers without inrush, (3) water inrush from separation layers into the goaf, and (4) water inrush from separation layers at the working face. Effective mitigation requires a mechanism-based prevention strategy tailored to each water inrush type, integrating drainage, monitoring, early warning, and engineering reinforcement. This targeted approach enhances risk control, ensuring safer and more sustainable mining operations. Future research will focus on refining the determination of the ks coefficient and incorporating fully coupled hydro-mechanical modeling to better capture fluid–solid interactions in TSPL failure and water inrush processes (Table 2).

Table 2. Physical and mechanical parameters of the different lithologies.
Lithology Elastic modulus (GPa) Poisson's ratio Tensile strength (MPa) Cohesion (MPa) Friction
Coarse sandstone 27.17 0.30 5.46 30.12 0.70
Conglomerate rock 17.17 0.20 6.47 26.78 0.70
Mudstone 15.00 0.50 3.50 2.50 0.40
Fine sandstone 25.00 0.25 7.50 4.70 0.48
Siltstone 14.35 0.50 2.50 3.40 0.43
Sand–mudstone interlayer 21.43 0.35 4.75 3.76 0.50
Coal 25.00 0.25 7.50 4.70 0.48
Interlayers 1.00 0.50 0.10 0.10 0.10

ACKNOWLEDGMENTS

This work was financially supported by the National Natural Science Foundation of China (Grant no. 42472334), the DeepEarth Probe and Mineral Resources Exploration-National Science and Technology Major Project (Grant no. 2024ZD1004208), and the China Postdoctoral Science Foundation (Grant no. 2025M771774).

    CONFLICT OF INTEREST STATEMENT

    The authors declare no conflict of interest.

    Biography

    • image

      Wei Qiao, PhD, professor and doctoral supervisor, currently serves as vice dean of the School of Resources and Geosciences, China University of Mining and Technology (CUMT). He obtained his PhD in Geological Engineering from CUMT in 2011, completed postdoctoral research at the Shandong Energy Group (2012–2015), and was a visiting scholar at the University of Toronto (2017–2018). He is a recipient of the “Gu Dezhen Youth Award” of the Geological Society of China, a selected academic leader of the Jiangsu Qinglan Project, and serves as an evaluation expert for the Ministry of Education. He is also a member of several international societies, including IAEG, IMWA, ISRM, and IAH. He has led over 30 national and enterprise-funded projects, participated as a key researcher in national key programs, and published more than 40 SCI/EI-indexed papers. He has co-authored one monograph, contributed to the Handbook of Coal Mine Water Prevention, and holds over 20 invention patents. He has received multiple provincial- and ministerial-level science and technology awards, including two First Prizes and two Second Prizes.




    附件【Deep Underground Science and Engineering - 2026 - Liu - Fracture evolution of a thick soft protection layer and the water.pdf