From Stress Redistribution to Energy Accumulation: Lateral Pressure-Driven Chain Evolution in Gas-Bearing Coal

1 Introduction

As shallow coal resources are progressively depleted, underground coal mining worldwide is increasingly extending to greater depths [1]. Under deep mining conditions, the combined effects of high in situ stress, high gas pressure, and mining-induced disturbance substantially increase the risk of dynamic disasters such as coal and gas outbursts and rock bursts [2]. The multifield coupled response mechanism of gas-bearing coal seams under mining-induced disturbances has thus become a critical challenge hindering safe deep mining. Underground deep coal mining is a complex multifield coupling problem involving multiple physical processes, such as mechanical stress, fluid flow, heat transfer, acoustics, and chemistry [3,4]. Previous studies, which are based primarily on elastic or elastoplastic theories, focused on the two-way coupling between stress and seepage fields [5,6,7]. However, during actual mining operations, coal‒rock masses inevitably undergo damage degradation under stress, significantly altering their mechanical properties and permeability. Consequently, models that treat damage as an independent physical field—exploring its stress-driven evolution and its subsequent feedback on the seepage field—are becoming increasingly mainstream approaches for investigating disaster mechanisms in deep coal-rock masses [8,9]. Li and Zhou [10] derived a smoothed particle hydrodynamics formulation of the seepage equation and a two-dimensional stress–seepage–damage coupled constitutive model. Zhao et al. [11] investigated the coupled effects of seepage and damage in fractured rock masses using a dual-media model based on equivalent continuous and discrete fracture media. Collectively, these studies highlight the critical role of damage evolution in driving dynamic changes in other physical fields and the initiation of disasters.

The occurrence of dynamic hazards in underground coal mining is affected by multiple factors, with geological structural features constituting a significant category. These features include compressive structures, tensile structures, fault structures, fold structures, and heterogeneous lithological distributions. The lateral pressure coefficient ( ξ ) can be used to characterize the complexity and heterogeneity of the underground geological structural distribution, which has a decisive effect on the stress distribution within the roadway surrounding rock and on roadway stability [12,13,14]. Zhao et al. [15] established a coupled strata–coal seam excavation model, simulating plastic zone expansion at different lateral pressure coefficients. They reported that the greater is the deviation in the lateral pressure coefficient from 1, the more significant the asymmetric failure of the surrounding rock. Wang et al. [16] integrated meso-damage mechanics, statistical strength theory, continuum mechanics, and digital image correlation to determine the lateral pressure effect in columnar jointed basalt. Zheng et al. [17] investigated the influence of the lateral pressure coefficient on the damage state and the spatiotemporal evolution characteristics of the gas extraction flow field.

Previous studies have shown that the lateral pressure coefficient significantly affects the stress distribution of roadway-surrounding rock, plastic zone morphology, fracture propagation, and gas migration behavior [18,19,20,21]. However, most studies have focused on a single physical field or on local coupling relationships, whereas the continuous transmission processes among stress redistribution, damage initiation and propagation, permeability evolution, gas pressure reconfiguration, and energy accumulation have not been systematically characterized [22,23,24,25,26]. In particular, how the lateral pressure coefficient controls the dominant direction, temporal sequence, and synergistic intensity of this process remains ambiguous. In this study, the term “chain evolution path” refers to the progressive coupled process in which stress redistribution induced by excavation promotes damage evolution, damage further alters seepage pathways, and seepage evolution regulates gas pressure and energy conditions. Elucidating this process is important for establishing differentiated hazard control strategies under different geological stress regimes.

Accordingly, this study develops an isothermal seepage–stress–damage coupling framework for gas-bearing coal and applies it to roadway excavation at different lateral pressure coefficients. The main purpose of the present model is to elucidate the relative evolution laws and controlling mechanisms of the coupled fields at different lateral pressure coefficients ( ξ values of 0.5, 0.8, 1.0, 1.2, and 1.5) rather than to provide a strict one-to-one prediction for a specific engineering site. By combining the characteristic response observed in triaxial tests with roadway-scale numerical simulations, the evolution of stress, deformation, damage, permeability, gas pressure, and elastic energy density is analyzed, and the controlling effect of the lateral pressure coefficient on the evolution path of the multifield chain is discussed.

2 Seepage-Stress-Damage Coupling Model of Gas-Bearing Coal

The present model is established within an isothermal continuum-damage framework to describe the coupled interaction among deformation, damage, and gas seepage in gas-bearing coal while maintaining computational tractability. In this framework, the damage variable is used as an internal state parameter to characterize material degradation and to control stiffness reduction and permeability enhancement, whereas energy analysis is used mainly to characterize the relative variation in local energy-storage conditions during excavation disturbance. A strictly thermodynamically consistent constitutive derivation based on an explicit free-energy potential is beyond the scope of the present study and will be further explored in future work.

Coal and gas outbursts represent a typical multi-field coupling process. When thermal effects are neglected, this physical phenomenon involves gas desorption/adsorption in coal, deformation of the coal matrix itself, damage evolution, and permeability changes. The stress field primarily alters seepage field parameters, including porosity ( ϕ ) and permeability ( κ ). In return, the seepage field influences the stress field through pore pressure ( p ) from free-state gas and coal matrix deformation induced by gas adsorption/desorption ( ε s ). Strictly speaking, the damage field is not an independent physical field but rather an extension of the stress field. However, for clarity of discussion, it is treated as a separate physical field. The stress field determines the damage variable ( D ) through damage criteria and assigns it to the damage field. The damage field, in turn, causes the degradation of physical parameters, which subsequently affects both seepage field parameters ( κ ) and stress field parameters ( E ).

2.1 Governing Equation for Stress Field

The fluid-solid coupling damage process in gas-bearing coal is described based on the model developed by Zheng et al. [27]. The governing equations for the seepage field and solid mechanical field are presented below. E21+νui,jj+E21+ν1−2νuj,ji−αp,i−E31−2νεs,i+fi=0(1) where E and ν are the elastic modulus and Poisson’s ratio of the coal, respectively; ε s represents the volumetric strain induced by gas adsorption in the coal matrix; p denotes the gas pressure; α is the Biot’s coefficient, α = 1 − K / K s ; K and K s are the bulk modulus of the coal and coal grain, respectively; f i represents the external force acting on the system.

2.2 Governing Equation for Seepage Field

The seepage field incorporates both free-state gas migration and the adsorption/desorption processes of adsorbed gas. The gas mass balance equation can be described as follows: ∂m∂t+∇ρg⋅μ⇀=Qs(2) where, m is the gas content in coal; ρ g represents the density of gas; μ ⇀ is the fluid velocity; Q s denotes the internal gas mass source/sink term. Since no internal gas injection or active extraction source was imposed in the present simulations, and gas migration was controlled only by the boundary conditions together with adsorption/desorption-related storage effects, Q s was set to 0.

The gas content consists of free gas and adsorbed gas. The gas content in coal can be expressed as: m=ρgϕ+ρgaρcVsg(3) where, ρ g a is the density of the gas under atmosphere pressure; ρ c is the density of the coal; V s g is the content of the adsorbed gas.

The change of coal permeability with porosity follows cubic law, the evolution model for the permeability and porosity of gas-bearing coal is given by Zhang et al. [5]: ϕ=11+S1+S0ϕ0+αS−S0(4) κ=κ0ϕϕ03(5) where, S = ε v + p / K s − ε s , S 0 = p 0 / K s − ε s 0 . The subscript, 0, denotes the initial value of the variable.

The permeability evolution law adopted in this study is mainly used as an effective phenomenological relationship to capture the first-order influence of porosity variation and damage-induced crack opening on seepage capacity. The focus of the present study is on the comparative seepage response under different lateral pressure coefficients, rather than on a strict microscale representation of the matrix–cleat dual-porosity structure of coal. More rigorous parameter calibration and physical refinement based on independent triaxial gas-permeability tests remain important tasks for future work.

Based on the gas mass balance equation and adsorption-induced expansion strain, the governing equation for fluid flow in gas-bearing coal is derived as follows: ϕcρ−α−ϕεLPL1+Sp+PL2+paρcVLPLpp+PL2+α−ϕKs1+S∂p∂t−∇⋅κμ∇p−cρκμ∇p2=Qs−α−ϕp1+S∂εv∂t(6) where, ε L , V L and P L are the Langmuir constant of strain, volume and pressure, respectively; ε v are the volumetric strain of the coal; ϕ and ϕ 0 are the porosity and initial porosity of the coal; κ and κ 0 are the permeability and initial permeability of the coal; c ρ represents the compressibility coefficient of the gas.

The primary difference between various fluids lies in their compressibility coefficients. The isothermal compression coefficient of gas c ρ can be described as follows: cρ=1p−1C∂C∂p(7) where, C is the compressibility factor of the gas.

2.3 Governing Equation for Damage Field

Research indicates that damage in gas-bearing coal leads to the deterioration of its mechanical properties and alters its internal pore-fracture structure, including permeability and porosity, which in turn affects gas seepage behavior [28]. Under complex stress states, tensile and shear damage are evaluated using the maximum tensile stress criterion and the Mohr-Coulomb criterion, respectively.

When the stress state of a meso-element satisfies the maximum tensile stress criterion, tensile damage occurs; when it satisfies the Mohr-Coulomb criterion, shear damage occurs. Tensile damage is prioritized—if a meso-element satisfies both criteria simultaneously, tensile damage is assigned. The damage variable ( D ) under complex stress states is expressed as [29]: D=0F1<0,F2<01−εt0εt2F1=0,dF1>01−εc0εc2F2=0,dF2>0(9) where, ε t 0 and ε c 0 are the ultimate tensile strain and ultimate compressive strain at failure under tension and compression; ε t and ε c are the equivalent principal strains under tensile and compressive states.

Based on elastic damage theory, the effect of damage on the elastic modulus of coal can be expressed as: E=1−DE0(10) where, E 0 is the initial elastic modulus. When D = 1 , the elastic modulus becomes 0, which can severely impact computational convergence in numerical simulations. Therefore, a minimum residual value can be assigned to the elastic modulus in this formulation to avoid numerical instability.

Considering the influence of damage on permeability [30], Eq. (5) can be expressed as follows: κ=κ0ϕϕ03expλαD(11) where, λ α represents the damage coefficient.

3 Experimental Verification and Numerical Model Establishment

3.1 Experimental Procedure and Numerical Model Development

Under deep engineering conditions, triaxial compression tests and corresponding numerical simulations were conducted on gas-bearing coal. A comparative analysis of the failure process was performed to verify the ability of the constructed model to simulate the coupled seepage–stress–damage response of gas-bearing coal.

The coal blocks were machined into standard cylindrical coal samples with dimensions of 50 mm × 100 mm Φ × H , whose end flatness and perpendicularity complied with the ISRM recommended standards. After measuring the basic physical parameters of the coal samples (mass, P-wave velocity), specimens with similar parameters were selected for testing. The average mass of the prepared samples was 290.66 g, and the average P-wave velocity was 1.882 km/s.

The experiments were performed using an electro-hydraulic servo-controlled triaxial seepage testing system. The system consists of a rock mechanics loading system and a gas seepage system. The schematic diagram of the experimental setup is shown in Fig. 1.

• (1)The rock mechanics loading system has a maximum axial load of 2000 kN with a measurement accuracy of 1%, a maximum displacement of 100 mm, and a maximum confining pressure of 80 MPa. It offers multiple control modes (stress, strain, and displacement) and is equipped with a mechanical data acquisition system that collects axial and circumferential displacement and load data.
• (2)The gas seepage system is a customized gas seepage/pressure boosting system composed of a gas source, air compressor, control cabinet, flow meter, and booster pump. It features multiple special internal seals to prevent gas leakage effectively. The gas source pressure ranges from 0 to 5 MPa, and the working pressure ranges from 0 to 100 MPa. The system can operate with media such as carbon dioxide, methane, high-purity nitrogen, and air. It is equipped with a data acquisition system to monitor the inlet/outlet pressure and flow rate.

Figure 1: Experimental system diagram.

To simulate both shallow and deep mining environments, two confining pressure levels were set at 5 MPa and 25 MPa respectively (assuming an average overburden density of 2500 kg/m3), with a fixed gas pressure of 4 MPa. The coal sample assembly is illustrated in Fig. 2. The experimental procedure consisted of four steps:

• (1)Sample installation: The specimen was fixed between upper and lower seepage pressure heads and sealed with heat-shrink tubing. Measuring devices were installed on the specimen surface, which was then placed in the triaxial pressure chamber. Test signal lines and gas inlet/outlet pipelines were connected.
• (2)Confinement and gas pressure application: After installation, the pressure chamber was sealed. A 2 kN preload was applied, followed by oil filling into the chamber. The confining pressure was increased to the target value (5 MPa or 25 MPa) and maintained throughout the test. A constant gas pressure of 4 MPa was introduced through the seepage control system, with the inlet pressure maintained stable while monitoring the outlet gas flow rate.
• (3)Gas saturation and loading: After the outlet flow rate stabilized, the outlet was closed to allow 24 h for full gas adsorption (during which circumferential displacement changes became essentially constant, indicating adsorption completion). The loading program was then initiated under displacement control at 0.3 mm/min until specimen failure.
• (4)Post-test procedure: After specimen failure, the test was terminated by sequentially releasing gas pressure, confining pressure, and axial pressure. The specimen was then dismantled, data saved, and preparation made for the next test.

Figure 2: Coal sample installation diagram and geometric model for numerical simulation.

Numerical simulation of the triaxial compression process for gas-bearing coal was conducted using COMSOL Multiphysics software. The geometric model, boundary conditions, and loading process adopted in the simulation were all referenced from the experimental procedures. To characterize the heterogeneity of the coal matrix, the Weibull distribution (Eq. (12)) was employed to describe the mechanical parameters of the coal in its initial state. fu=mu0uu0m−1exp−uu0m(12) where, u represents the physics parameters, u 0 denotes the corresponding scale parameter, which reflects the average value of all discrete element parameters in the coal body; and m is the shape parameter, characterizing the heterogeneity of the coal mass.

The mechanical parameters of the coal and seepage characteristic parameters used in the numerical model are listed in Table 1. The physical properties of coal and gas were obtained from previous experimental studies [27,28]. The shape parameters were calibrated using the meso-element parameter calibration method reported by Zhou et al. [31], whereas the damage coefficient was adopted from Zhu’s research [30]. The numerical simulation was conducted in two steps: (1) the coal was completed saturated with gas and (2) the same loading rate as the experiment was used. The boundary conditions for the stress field were as follows: initial displacement was set to zero; constant confining pressure was applied to the left and right boundaries; a roller support condition was applied to the bottom boundary to limit horizontal displacement; and a prescribed axial displacement was applied to the top boundary to simulate the displacement-controlled loading process in the experiment. The boundary conditions for the seepage field were as follows: the initial pressure was set to atmospheric pressure; a specified gas pressure was applied to the bottom boundary throughout the simulation. Because the outlet was sealed and the sides were enclosed after gas flow began in the experiment, the top and side boundaries in the numerical model were established as zero-flux boundaries. Most current studies [32,33] rely on the secondary development of software, where damage is evaluated and damage variables are calculated in MATLAB. The degraded parameters are imported into COMSOL to calculate the damage. This method can only calculate the damage process within preset iterative steps. In contrast, in this study, the state variable function in COMSOL Multiphysics is used to treat the damage field as a state variable, enabling continuous tracking of its full-time evolution.

3.2 Comparative Analysis of Experimental and Numerical Results

The experimental results indicate that the peak load displacements of gas-bearing coal at confining pressures of 5 MPa and 25 MPa are 1.38 mm and 2.07 mm, respectively. These displacement values were set as the final loading magnitudes for the prescribed boundary displacement at the top boundary in the numerical simulation. The fragment size distribution and failure pattern characteristics of gas-bearing coal after loading at different confining pressures are shown in Fig. 3. After failure, the fragments were sieved and weighed, and the fractal dimension of the fragment size ( D f ) was calculated to characterize the distribution of failure fragments [34]. A higher D f of the fragment size indicates finer fragmentation and a greater degree of damage.

Figure 3: Comparison of damage evolution and failure pattern in triaxial compression of gas-bearing coal.

The numerical results indicate that damage in gas-bearing coal under triaxial loading evolves through the processes of crack initiation, propagation, coalescence, and localization. As the confining pressure increased to 25 MPa, these random damage points disappeared, corresponding to a decrease in the D f of the fragment size from 2.38 to 1.27 after failure, indicating a significant reduction in the degree of fragmentation. Under continued displacement loading, the damage zones progressively expanded and coalesced, resulting in the formation of macroscopically damaged areas, which aligned with crack propagation and penetration. At a confining pressure of 5 MPa, numerous randomly distributed damage points appear at the initial loading stage (displacement of 0.78 mm), indicating that local tensile–shear instability is activated at multiple locations. These dispersed damage nuclei promote crack branching and interaction, resulting in more tortuous crack propagation paths and ultimately finer fragmentation after failure. In contrast, at a confining pressure of 25 MPa, the randomly distributed damage points are markedly suppressed, and damage tends to localize along fewer dominant paths. As a result, a greater proportion of deformation energy is concentrated in the propagation and opening of the main cracks, resulting in more distinct macroscopic shear-dominated failure surfaces and a lower degree of fragmentation. Therefore, the macroscopic failure morphology observed in the experiments is an external manifestation of the spatial organization and localization degree of damage evolution in the numerical model.

Notably, the triaxial tests and the roadway-scale simulations correspond to different spatial scales and boundary conditions. In the present study, triaxial tests are used mainly to validate the ability of the model to reproduce the damage initiation, localization, and failure morphology evolution of gas-bearing coal under controlled isothermal conditions rather than to provide a strict one-to-one prediction of field excavation. Accordingly, the subsequent roadway-scale simulation results should be interpreted primarily as mechanism-oriented comparative analyses under simplified conditions. In actual engineering environments, temperature variation, structural discontinuities, excavation rate, support installation sequence, and active gas drainage may further affect the coupled response. These factors, in conjunction with field-based calibration, will be integrated into future studies.

The comparative analysis between the numerical simulations and the experimental results of triaxial compression of gas-bearing coal validates the ability of the model to simulate the characteristic damage evolution and failure morphology. This approach provides a basis for further investigation of the coupled seepage–stress–damage evolution process in deep coal mining.

4 Numerical Simulation of Deep Coal Mining Processes

4.1 Numerical Model Development

To investigate the controlling effect of the lateral pressure coefficient on the multifield chain evolution path at the engineering scale, a roadway excavation model was further established after the sample-scale response characteristics were examined. This section presents the model geometry, boundary conditions, excavation loading scheme, and parameter settings used for the comparative analysis. The main purpose of this model is to identify, under unified simplified conditions, how different lateral pressure coefficients affect the dominant direction and evolution characteristics of the stress, damage, seepage, and energy responses of the surrounding rock.

A schematic of the numerical model is shown in Fig. 4. The underground structure is simplified into overlying strata, a coal seam, and a floor stratum. The roadway is situated within the coal seam. For simplification, the analysis focuses on a cross-section of the roadway. The roadway is composed of a 2.5 m × 5 m rectangle and a semicircle with a radius of 2.5 m. The boundary conditions for the stress field are configured as follows: The left and bottom boundaries are roller supports. The right boundary is subjected to a lateral load, the magnitude of which is determined by ξ . The top boundary is under a vertical load that represents the overburden pressure. This study focuses on coal seam mining at kilometer depth, and the vertical load is calculated at average overburden density of 2500 kg/m 3 . The boundary conditions for the seepage field are defined as follows: All four lateral boundaries are zero-flux boundaries. The roadway perimeter is set to atmospheric pressure. The initial gas pressure within the coal seam is 4 MPa. The investigation targets the fluid‒solid damage coupling process. Therefore, during the simulation of coal seam extraction, no artificial pressure relief measures are applied to the coal seam gas, considering only its spontaneous seepage process. To simulate the roadway excavation process, during which the stress state evolves as the working face advances, the roadway excavation process is simulated by applying a linear load to the roadway cross-section. In the initial state, the roadway perimeter is under in situ stress. After excavation completion, the roadway surface becomes exposed with zero constraint.

Figure 4: Schematic diagram of numerical model.

To ensure model convergence, both mechanical constraints and atmospheric boundary conditions were applied as linear functions. The excavation process was configured for a duration of 30 min to monitor roadway stability during excavation and subsequent flow field evolution. ξ was assigned values of 0.5, 0.8, 1.0, 1.2, and 1.5. Considering the heterogeneity of the coal seam, its elastic modulus is assumed to exhibit a Weibull distribution (Eq. (12)), as shown in Fig. 5. The excavation process was numerically simulated using a fully coupled seepage‒stress‒damage model. All the parameters employed in the numerical model are listed in Table 1.

Figure 5: Distribution of elastic modulus in the coal seam around the roadway.

4.2 Stress Field Evolution during Mining

Roadway excavation disrupts the initial stress state, resulting in stress redistribution. Fig. 6 presents the evolution of the stress distribution during roadway excavation at different lateral pressure coefficients. The lateral pressure coefficient significantly affects the stress distribution patterns. At ξ = 0.5 and ξ = 1.5 , distinct stress concentration zones emerged after excavation was completed, whereas the other coefficients resulted in approximately annular high-stress zones around the roadway. At all the coefficients, the initial stress concentrations were concentrated at the junctions between the roof, floor and the sidewalls. At ξ = 0.5 , after the initial stress concentration, the high-stress zone progressively expanded horizontally along the sidewalls, and the maximum difference between the horizontal and vertical stresses was nearly one order of magnitude. In contrast, at ξ = 1.5 , the high-stress zone tended to develop vertically along the roof and floor. These findings indicate that pressure relief measures should be adapted to specific geological conditions during underground mining operations.

Figure 6: The stress distribution around the roadway under different lateral pressure coefficients (a) ξ = 0.5  (b) ξ = 0.8  (c) ξ = 1.0  (d) ξ = 1.2  (e) ξ = 1.5 .

Roadway deformation directly reflects the stress-induced deformation effects and collapse risk. Fig. 7 presents the deformation curves of the roof, floor, and the sidewalls during excavation at different lateral pressure coefficients (with the arc length as the unified coordinate). At ξ = 0.5 , both the roof, floor and the sidewalls exhibited maximum displacement in all the scenarios: the maximum roof deformation was 0.76 m at the center, the maximum floor deformation was 0.33 m at both ends, and the maximum sidewall deformation was 0.10 m in the middle–lower section. The deformation at different lateral pressure coefficients was generally similar, with the magnitude of deformation decreasing as the lateral pressure coefficient approached 1. Conversely, at ξ = 1 , significant deformation occurred mainly in the roof area, whereas the sidewall and floor deformations were relatively small, indicating a more uniform stress distribution.

Figure 7: Deformation curve of roadway subsidence under different lateral pressure coefficients.

4.3 Damage Field Evolution during Mining

The evolution of damage to coal seams induced by mining disturbance is an extremely critical parameter. The extent of damage not only affects the physical properties of the coal mass itself and results in the propagation of internal fracture networks but also significantly alters gas storage and migration characteristics. Fig. 8 shows the evolution of the damage distribution around the roadway during excavation at different lateral pressure coefficients.

At ξ = 0.5 , where the lateral pressure (12.5 MPa) is relatively close to the initial gas pressure (4 MPa), the coal is significantly affected by gas, resulting in substantial damage during the initial excavation stage. Analysis of the damage distribution at other coefficients reveals that, from the perspective of damage crack initiation, damage consistently initiates at the junctions between the roof and floor and the sidewalls.

At ξ = 0.8 , owing to the lower horizontal stress, damage initiates from the sidewalls and gradually propagates, with a specific degree of damage also present in the floor. At ξ = 1 , damage develops almost simultaneously in the sidewalls, roof, and floor, ultimately creating an annular damage zone around the roadway. At ξ > 1 , damage primarily initiates and propagates from the roof and floor, with only small damage zones emerging in the sidewalls.

With respect to damage characteristics, the patterns are similar across different lateral pressure coefficients. The damage distribution is relatively uniform around the roadway perimeter and exhibits a radial propagation pattern as it extends inward.

Figure 8: Distribution of damage area around roadway with different lateral pressure coefficients (a) ξ = 0.5  (b) ξ = 0.8  (c) ξ = 1.0  (d) ξ = 1.2  (e) ξ = 1.5 .

Global parameters of the model are calculated, the computation is characterized using the nodal average, with the calculation method as follows: Par=∑i=1nPari/n(13) where, P a r represents a category of global parameters; n is the node numbers.

The evolution of the total damage value in the computational model was analyzed, with a focus on the process before excavation was completed. As shown in Fig. 9, the ξ = 0.5 scenario was excluded from the analysis because of its substantially high initial damage value. The total damage value increases with increasing lateral pressure coefficient. The damage values at 30 min ( ξ = 0.8, 1.0, 1.2, and 1.5) are 0.023, 0.025, 0.029, and 0.036, respectively.

Figure 9: Damage evolution before the end of excavation under different lateral pressure coefficients.

Fig. 10 presents the evolution patterns of elastic energy around the roadway at different lateral pressure coefficients. The magnitude of elastic energy in the coal seam represents its energy storage capacity, with high-energy zones being more prone to outbursts. Therefore, we analyze both the distribution and total density of elastic energy in the coal seam. In terms of magnitude, the ξ = 0.5 and ξ = 1.5 scenarios exhibit the most extensive high elastic energy zones. When the lateral pressure coefficient deviates from 1, the stress concentration increases, increasing the energy accumulation capacity. In terms of spatial distribution, at low lateral pressure coefficients, high elastic energy zones are distributed predominantly in horizontal regions, whereas at high lateral pressure coefficients, they are concentrated around the roof and floor areas.

The total elastic energy density of the model was also calculated and analyzed. Fig. 11 presents the evolution curves of coal seam elastic energy density during excavation at different lateral pressure coefficients. Similar to the damage evolution process, the ξ = 0.5 scenario has the highest elastic energy density of 0.378 MJ / m 3 , which results from the competitive mechanism between gas pressure and low lateral pressure. As the ξ value increases from 0.8 to 1.5, the corresponding elastic energy density increases from 0.073 MJ / m 3 to 0.172 MJ / m 3 .

Figure 10: Evolution of elastic energy around roadway under different lateral pressure coefficients (a) ξ = 0.5  (b) ξ = 0.8  (c) ξ = 1.0  (d) ξ = 1.2  (e) ξ = 1.5 .

Figure 11: Elastic energy density evolution before the end of excavation under different lateral pressure coefficients.

4.4 Seepage Field Evolution during Mining

Gas pressure is a critical factor affecting coal and gas outbursts. Its distribution affects the stability of the surrounding rock and affects the selection of dynamic disaster prevention measures. Fig. 12 shows the gas pressure distribution around the roadway during coal seam extraction and after stabilization (12 h after excavation) at different lateral pressure coefficients. The reduction in gas pressure directly correlates to the volume of gas emission. In terms of gas pressure magnitude, at ξ = 0.5 , damage significantly improves coal seam permeability, resulting in rapid gas emission. Consequently, the gas pressure nearly equilibrates with atmospheric pressure within 12 h after excavation. Although the final gas pressure magnitudes are relatively consistent at other lateral pressure coefficients, the distribution of low-pressure zones differs significantly. With respect to the spatial distribution of gas pressure, at ξ = 0.8, 1.0, and 1.2, the gas pressure decreases more slowly, and its distribution pattern remains relatively consistent, gradually extending from the roadway periphery to the deeper coal seam during excavation. In contrast, at ξ = 1.5 , the low-pressure zone gradually expands vertically from the roof and floor into the coal seam, and only limited low-pressure areas develop along the sidewalls.

Figure 12: The gas pressure change diagram around the roadway after coal seam mining under different lateral pressure coefficients (a) ξ = 0.5  (b) ξ = 0.8  (c) ξ = 1.0  (d) ξ = 1.2  (e) ξ = 1.5 .

The average gas pressure in the coal seam after roadway excavation was calculated, with stabilization occurring 12 h after excavation. The pressure evolution is shown in Fig. 13. At ξ = 0.5 , the gas pressure decreases at a significantly rapid rate. During the stabilization phase, the pressure approaches zero, at which point the outburst risk decreases significantly, and the hazard mechanism shifts to rock burst. For lateral pressure coefficients of 0.8, 1.0, 1.2, and 1.5, the corresponding stabilized pressures are 3.57 MPa, 3.53 MPa, 3.49 MPa, and 3.39 MPa, respectively, and all the reductions remain within 20%.

Figure 13: The change of gas pressure after coal seam mining under different lateral pressure coefficients.

The distribution of coal seam permeability directly controls both the gas pressure gradient and flow pathways. High-permeability zones accelerate the reduction in the gas pressure gradient, whereas low-permeability areas can create high-pressure regions, increasing outburst risk. Fig. 14 presents the permeability evolution after coal seam extraction at different lateral pressure coefficients. Because permeability is affected by the extent of damage, a widespread high-permeability zone forms throughout the coal seam at ξ = 0.5 . With respect to the spatial evolution of permeability, at ξ = 0.8 and ξ = 1.0 , high-permeability zones initiate and expand from the vicinity of the sidewalls. After stabilization, these zones are concentrated near the sidewalls and floor. Conversely, at ξ = 1.2 and ξ = 1.5 , high-permeability zones initiate and propagate from the roof and floor areas. After stabilization, the high-permeability zones extend mainly to the roof and floor.

Figure 14: Permeability changes around the roadway under different lateral pressure coefficients (a) ξ = 0.5  (b) ξ = 0.8  (c) ξ = 1.0  (d) ξ = 1.2  (e) ξ = 1.5 .

5 Analysis of Multi-Field Coupled Chain Evolution Paths under Mining-Induced Disturbances

5.1 Multi-Field Chain Evolution Pathways

The coupled results indicate that the stress, damage, seepage, and energy fields during roadway excavation do not evolve through simple synchronous variation but instead exhibit a chain-like coupling relationship with asymmetric feedback characteristics. Stress redistribution induces this process because it determines the locations of stress concentration and the orientation of the principal stress difference, thereby controlling where damage initiates and is localized. Damage evolution is the key structural link through which the mechanical field is transmitted to the seepage field: as tensile or shear damage occurs, the local stiffness decreases, and the permeability increases, which promotes gas pressure redistribution and alters the effective stress state. The energy field reflects the variation in local energy storage conditions throughout this process. Notably, the elastic energy density considered in this study is used mainly as a comparative indicator of local energy accumulation capacity under different conditions rather than as a rigorous decomposition of fracture-related dissipation energy. Therefore, the spatial overlap among high-stress, high-damage, and high-energy zones suggests that energy accumulation and release participate in the coupled evolution process, whereas the exact partition between recoverable strain energy and damage dissipation remains underexplored. The chain evolution path revealed in this study can be summarized as follows: stress redistribution–damage initiation and localization–seepage pathway adjustment–gas pressure reconfiguration, accompanied by the synchronous evolution of local elastic energy density.

5.2 Control Effect of Lateral Pressure Coefficient on Multi-Field Evolution Paths

The lateral pressure coefficient, a key parameter that characterizes geological structures, significantly affects the multifield coupling process in deep coal mining. This coefficient regulates the dominant direction of multifield evolution by affecting the initial stress state and constraint conditions of the surrounding rock.

At low lateral pressure coefficients ξ < 0.8 , in this scenario, the difference in the horizontal load from the gas pressure (4 MPa) is relatively small, resulting in a weaker horizontal constraint on the surrounding rock. After excavation and unloading, the unloading effect in the horizontal direction becomes more distinct. The difference between the confining pressure and gas pressure results in a significantly steeper horizontal stress gradient. Consequently, damage initiation, crack propagation, the formation of high-permeability zones, and energy preferentially occur in the horizontal direction (sidewalls), as shown in Fig. 8a, Fig. 10a and Fig. 14a. The horizontally connected high-permeability channels provide efficient pathways for gas flow, enabling rapid gas pressure release (Fig. 12a). While this approach reduces the risk of outbursts, it simultaneously increases the horizontal deformation and instability of the surrounding rock.

At balanced lateral pressure coefficients 0.8 ≤ ξ ≤ 1.2 , with minimal structural heterogeneity and a relatively uniform in situ stress environment, the pressure differential between the confining pressure and gas pressure increases, whereas the constraining effect of the confining pressure is improved. Consequently, the development of damage and seepage fields exhibits notable symmetry, and the surrounding rock system maintains a state of relative equilibrium. The evolution of each physical field exhibits significant synergy, and the stress, damage, and seepage fields exhibit annular evolution patterns around the roadway, as shown in Fig. 8b–d, Fig. 10b–d and Fig. 14b–d, respectively.

At high lateral pressure coefficients ξ > 1.2 , the elevated horizontal load results in a relatively weaker vertical constraint. After excavation, the vertical unloading effect becomes dominant. The disturbance effect of gas pressure on the stress field is significantly decreased in a high in situ stress environment, thereby resulting in a vertical orientation of the maximum difference in principal stress. Consequently, both damage evolution and high-permeability zones develop primarily vertically along the roof and floor of the roadway (Fig. 8e and Fig. 14e).

While the vertically extended gas migration paths facilitate the fastest pressure drop rate compared with the other lateral pressure coefficients (except ξ < 0.8 ), the system tends to maintain a high-pressure state because of the strong geostress confinement. This approach promotes significant energy accumulation in the roof and floor areas (Fig. 10e), substantially increasing the potential risk of coal and gas outbursts.

The above results indicate that the lateral pressure coefficient not only changes the magnitude distribution of individual physical quantities but also changes the dominant direction and transmission efficiency of the coupled multifield process. The lateral pressure coefficient determines not only where damage and seepage preferentially develop but also how excavation-induced disturbance is transmitted among the stress, damage, seepage, and energy fields. This explains why different geological stress conditions may correspond to different dominant hazard patterns; therefore, differentiated prevention and control strategies are needed.

5.3 Differentiated Hazard Prevention and Control Strategies

The findings of this study clearly guide disaster prevention and control in deep coal mines. In practical engineering applications, the lateral pressure coefficient and gas pressure of a coal seam can be determined through field measurements. This approach enables the prediction of dominant disaster types and high-risk zones during mining operations, allowing for the development of corresponding prevention measures.

With respect to mining areas with a lateral pressure coefficient of ξ < 0.8 , prevention and control strategies should focus on stability issues such as significant deformation of the sidewalls. Support designs must increase the reinforcement strength of the sidewalls while leveraging the characteristics of rapid gas pressure release to achieve efficient gas extraction.

With respect to mining areas with a lateral pressure coefficient of 0.8 ≤ ξ ≤ 1.2 , where both the energy accumulation capacity and gas pressure maintain intermediate levels, prevention strategies should simultaneously address gas pressure relief and support reinforcement in both the sidewalls and roof/floor of the roadway. This can be achieved by implementing radially arranged extraction boreholes around the perimeter as well as uniformly strong bolting-mesh-cable support systems. These integrated measures ensure a uniform stress distribution in the surrounding rock while concurrently reducing the gas pressure around the roadway.

With respect to mining areas with a lateral pressure coefficient of ξ > 1.2 , particular attention must be given to the outburst risk in the roof and floor directions of the roadway. Targeted reinforcement measures must be implemented with priority given to releasing the accumulated energy in the coal‒rock mass. This approach involves implementing robust roof-floor reinforcement strategies vertically while utilizing long drilling boreholes for predrainage of gas or implementing fracturing-induced pressure relief in the sidewall direction, thereby effectively disrupting the high-energy accumulation and high gas pressure gradient environment that promotes disasters.

Support design for roadways in deep mines should consider the directional characteristics of in situ stresses by employing asymmetric support structures and differentiated pressure relief strategies, enabling prevention and control measures for site-specific conditions.

6 Conclusions

On the basis of a coupled seepage–stress–damage model for gas-bearing coal, this study investigated the evolution of the stress, damage, seepage, and energy fields during roadway excavation in deep coal seams at different lateral pressure coefficients. Under the isothermal assumption and simplified boundary conditions, the main conclusions are as follows:

• (1)During roadway excavation in deep coal seams, the stress, damage, seepage, and energy fields clearly exhibit a coupled evolution pattern. An excavation disturbance induces stress redistribution, and the resulting stress concentration further induces damage initiation and propagation. Damage evolution affects gas pressure distribution by altering permeability and seepage pathways, whereas gas pressure variation feeds back to the mechanical state of the coal through the effective stress mechanism. The spatial overlap among high-stress zones, high-damage zones, and high elastic-energy-density zones suggests that energy accumulation is critical in this coupled instability process.
• (2) By altering the initial stress distribution and unloading constraint conditions of the surrounding rock, the lateral pressure coefficient significantly affects the dominant direction and spatial pattern of the multifield evolution. At relatively low lateral pressure coefficients ξ<0.8, the damage zone, high-stress zone, and low-gas-pressure zone preferentially develop along the sidewalls. At near-balanced lateral pressure coefficients 0.8≤ξ≤1.2, these fields exhibit relatively coordinated annular expansion around the roadway. At relatively high lateral pressure coefficients ξ>1.2, the corresponding zones extend mainly to the roof and floor.
• (3) In terms of mechanism-oriented hazard prevention, roadway support and pressure-relief measures should be configured in a differentiated manner according to the lateral pressure coefficient. With respect to ξ<0.8, sidewall stability control should be emphasized. With respect to 0.8≤ξ≤1.2, balanced support and pressure-relief measures around the roadway are needed. At ξ>1.2, reinforcement and directional energy-release measures in the roof–floor direction deserve particular attention.
• (4) The present model is established under isothermal conditions and simplified boundary assumptions, and triaxial tests are mainly used to verify the characteristics of damage evolution and failure morphology. Therefore, the present results are more appropriately interpreted as a comparative mechanistic analysis rather than as a strict prediction for specific engineering conditions. Future work will further address thermodynamically consistent constitutive formulation, the calibration of coupled permeability parameters, sensitivity and uncertainty analysis, and field-scale validation.

The present findings provide a mechanism-oriented reference for understanding the coupled response of deep gas-bearing coal seams and for developing differentiated hazard-control strategies.