A numerical study of the mining‐induced energy redistribution in a coal seam adjacent to an extracted coal panel during longwall face mining: A case study

Mining of the coal seam adjacent to the extracted coal panel is widely practiced in China to improve coal resource recovery rates, but the energy change associated with longwall mining may present a new risk to ground stability and gateroad failure. To understand the energy change effectively, a meticulously validated numerical model, built using FLAC3D software, incorporating a double‐yield model for the gob materials and a strain‐softening model for the coal/rock masses, was developed to investigate the energy redistribution in coal seams and barrier pillars during the process of mining the coal seam adjacent to the extracted coal panel. The model results showed that the closer the region was to the LW 5301 gob, the higher the location and magnitude of the peak strain energy density. Consequently, the risk of rockburst in coal seams was greater than that in barrier pillars. Along the coal seam strike, with increasing advancing distance from the setup room, the magnitude of the peak strain energy density gradually increased to 1.88‐2.78 times the premining energy level. In addition, along the coal seam dip, the maximum distance from the peak strain energy density to the edge of the coal seam was approximately 16‐28 m, which is greater than that in the barrier pillar. The proposed numerical simulation procedure and calibrated method could be a viable alternative approach to evaluate longwall mining‐induced energy changes.


| INTRODUCTION
In recent years, an increasing number of deep coal mines in China have employed a single-entry system, with a small pillar (4-9 m wide) left between the extracted coal panel and the adjacent coal seam, to improve the coal resource recovery rates or eliminate the coal burst risk. [1][2][3] Some studies have shown that a longwall retreat adjacent to an extracted coal panel leads to stress concentrations and energy redistributions around the openings, including gateroads and pillars, which creates unfavorable conditions resulting in gateroad instability and rockburst development. [4][5][6][7] Rezaei et al 8 and Jiang et al 9 stated that the ground stability and failure mechanisms of gateroads vary depending on the mining-induced energy changes. Thus, it can be expected that understanding mining-induced energy redistribution around openings provides fundamental knowledge for stability analysis of gateroads and evaluating coal burst risk.
Numerous studies have been performed using different methods to estimate the energy evolution due to excavation. [10][11][12][13] Some researchers have analyzed the energy change with in situ investigations based on feedback analysis of acoustic emissions and microseismic and stresses. [14][15][16] Their results indicate that energy changes due to excavation are associated with in situ stress, geological and geotechnical conditions. These observations provide first-hand data for studying advanced case studies where the energy change is important to evaluate. Apart from fieldwork, many studies employ analytical methods for understanding the energy evolution due to excavation. As early as the 1960s, Cook 17 noted that for elastic mediums, the loss of potential energy was exactly twice the strain energy stored in stress concentrations around an excavation. Salamon 18 proposed a novel analytical procedure to quantify the relationship between the energy components (including strain energy, released energy, and stored energy). Dong et al 19 derived analytical expressions for the mining-induced energy redistribution around rectangular openings and then discussed the influence of the initial stress coefficient and excavation shape on the energy distribution. These studies have made certain meaningful achievements in quantifying energy changes. However, all of these studies are focused on the energy change around underground small-scale excavations. At present, many researchers have analyzed energy redistribution around underground mining excavations with numerical modeling methods. 10,20,21 Mitri 22 assessed the mining-induced energy storage rate with finite element software while assuming ideal elasticity of the surrounding rock. Some researchers combined numerical approaches with fieldwork to investigate the energy redistributions in circular openings. 23 For underground longwall face excavations, Wang et al 24 presented a numerical investigation to assess the energy redistribution in an island longwall panel. In his numerical model, very elastic soft material was used to simulate the gob material. Liu et al 25 discussed the accumulation and dissipation law of energy during the process of ascending mining using FLAC3D software with no consideration of the goaf model. Wang et al 26 also investigated the energy change in a coal seam using a numerical model. His research was limited to the longwall face where there was no goaf. These studies indicated that the energy change in a coal seam varies depending on the in situ stresses, overburden movements, mining methods, and conditions. Due to complex geological conditions, the mining-induced energy change in the longwall face is still unclear, especially for coal seams adjacent to extracted coal panels.
To solve such a problem, this paper presents a case study on the mining-induced energy redistribution in a coal seam adjacent to an extracted coal panel using a numerical method. In the numerical model, the double-yield model was used to simulate the gob material with consideration of goaf reconsolidation, and the results of the global model were examined by comparison with field measurements. Then, the distribution of strain energy density was analyzed in detail. Finally, the effect of the extracted coal panel on the extraction scale was discussed.

| Geology and geotechnical overview of the Xinhe mine
The site for the field test was located in the Xinhe coal mine, which is located in Shandong province, China. The Xinhe coal mine is well known for its large number of longwall panels adjacent to mined-out areas. There are three mineable coal seams in the Xinhe coal mine. The targeted longwall face 5302 (LW5302), as shown in Figure 1, was selected for this study. In LW 5302, the longwall top coal caving (LTCC) method was used to extract coal seam No. 3. The mining and caving heights of coal seam No. 3 were 3.9 and 6.2 m, respectively. In Figure 1, it can be seen that LW 5302 was 120 m wide by 754 m long. The northwest side was adjacent to the longwall face 5301 (LW 5301) gob. The longwall top coal caving method was also used in LW 5301 to extract coal seam No. 3. After LW 5301 was completely mined out, the extraction of coal seam No. 3 was started in LW5302. In the LTCC operation, LW 5302 and 5301 were advanced along the strike direction of the coal seam.
Coal seam No. 3 was mined approximately −980 m deep. It was approximately 10.2 m thick on average and dipped 2-6° with an average of 3°. Figure 2 presents the generalized stratigraphic column of the study site, which was constructed based on core logging. As shown in Figure 2, the immediate roof, approximately 29.8 m thick, was composed of sandstone and siltstone, fine sandstone, and mudstone. The main roof consisted of fine sandstone and siltstone with an average thickness of 41.8 m.

| Rock mass material properties
Laboratory experiments were conducted on rock/coal core samples collected from the floor, roof, and coal seam of LW 5302. These core samples were extracted according to the ISRM recommended standard. The laboratory tests were performed in a servo-controlled testing system developed by the State Key Laboratory of Mining Disaster Prevention and Control co-founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology. The properties of the coal and rock samples obtained from laboratory tests are listed in Table 1.
Due to the engineering discontinuities in the rock mass (such as joints, cracks, bedding planes, and different mineral compositions), the strength of the rock mass was much smaller than that of the intact rock. [27][28][29] Therefore, it was necessary to calibrate the rock mass parameters. In the past, many rock mass classification systems have been proposed and used to evaluate and classify rock mass, such as the Q system, RQD index, rock mass rating (RMR) method, geological strength index (GSI) method, and R mi system. 30,31 The GSI could be estimated from rock mass descriptions by field observations and provides a complete set of mechanical properties using the generalized Hoek-Brown failure criterion. 32 The generalized Hoek-Brown failure criterion is expressed as: where σ 1 and σ 3 are the major and minor effective principal stresses; σ c is the uniaxial compressive strength of the intact rock material; and m b , s, and a are material constants. When the GSI is known, m b , s and a are expressed as: where m i is the intact rock parameter and D is a disturbance factor that varies from 0 for undisturbed in situ rock masses to 1 for very disturbed rock masses. Based on the Hoek-Brown envelope generated by solving Equation (1), the equivalent Mohr-Coulomb strength parameters c and f could be evaluated. The rock mass deformation modulus E rm and strength σ cmass were evaluated from the following equations given by Hoek et al (2006) as follows: For convenience, the RocLab software program, based on the generalized Hoek-Brown failure criterion, was used to determine the rock mass parameters. In this study, m i and σ c were obtained from laboratory experiments. The GSI was obtained from field geological descriptions. Because significant floor heave occurred in the gateroad of LW 5302 during its service duration, the value of D was set to 0.5. Note that the failure envelope range option was selected as the "tunnel" with a depth of 980 m. The rock mass properties are listed in Table 2.

GENERATION
In the following section, a local 3D numerical model was developed to study the energy change acting on coal seam No. 3 using the finite-difference software FLAC3D.

| Numerical simulation model
In this study, it was decided to model only the selected section, which can be considered the most critical case regarding energy redistribution. Considering that the purpose of the model was to study the energy change in coal seam No. 3, a 3D numerical model composed of two longwall faces, LW 5301 and LW 5302, was generated, as shown in Figure 2. The size of the model was 800 m (width) × 455.2 m (width) × 170 m (height). Numerous trials were conducted with reduced computation time due to an efficient model with a suitable mesh/grid density or mesh size around LW 5302. The overburden was assumed to be composed of two parts, a rock part 150 m in height and an ice part that was 980 m thick. The rock part of the overburden was established in the model, while the ice part was replaced by a constant vertical stress of 23.25 MPa at the top boundary of the model. No displacement in any side boundary was allowed, and the bottom boundary was fixed. The in situ stresses in this study were based on past experiences with the simulation of energy change during LW5301 mining. The in situ stresses were applied in the form of initial stress, and the at-rest pressure coefficients along the x and y direction were set as 1.15 and 0.95, respectively.

| Simulation plans
In this simulation study, the numerical model was solved using the following steps: (a) apply the in situ stress state and reach the initial equilibrium, (b) excavate the gateroad of LW 5301 along the strike and apply the primary support, (c) retreat LW 5301 along the strike until the excavation of LW 5301 is complete, and (d) retreat LW 5302 along the strike for each cycle once the analysis reaches equilibrium. The directions of the gateroad development and retreat mining are shown in Figure 3.

| The constitutive model and parameters for rock materials
In a numerical model, the input parameters, constitutive model, and yield criterion decisioning are essential. In previous studies, three groups of constitutive models (elastic, pseudoelastic, and plastic) and different yield criteria (eg, Mohr-Coulomb criterion, Hoek-Brown Criterion) were used to model coal and rock materials. 33,34 Based on past experiences in simulating coal and rocks, in this study, the strain-softening and Mohr-Coulomb models were adopted to simulate the mechanical behavior of the coal and roof/ floor rock layers. The rock mass mechanical properties for   the numerical simulation are listed in Table 2. Most of these properties were evaluated from intact rock properties, and the calculation method is presented in Section 2.2. Note that in the strain-softening model embedded in FLAC3D, the cohesion, friction, and dilation were defined as piecewise-linear functions of a softening parameter measuring the plastic shear strain. A plastic shear strain threshold of rock materials, which is required for the reduction of cohesion and friction from the peak to residual values, ranges from 0.01% 2 to 0.1%. 35 In general, there exist a few guidelines for the determination of the plastic shear strain threshold. In the current study, the plastic shear strain threshold was assumed to be 0.1%. Because it is generally difficult to estimate the residual cohesion strength, the peak cohesion value of the coal and roof/floor rock layers was scaled to the residual strength using a reduction factor of 0.1 based on past experiences in rock simulation. This factor had been commonly applied to most coal measure lithologies. In the postfailure phase, the friction angle and tensile strength remained unchanged based on a similar assumption in several previous studies. Table 2 presents the strain-softening parameters for the modeling.

| Double-yield model for the gob material
In longwall mining, the roof strata of the panel cave and fall into the gob. The caved-in materials are compacted and consolidated afterward, thereby altering the abutment load in the surrounding rock because a portion of the vertical load will be supported by the consolidated materials. 2 Therefore, the selection of a realistic gob-compaction model is an essential part of modeling longwall mining.

| Approach for gob modeling based on the double-yield model
Currently, there are two approaches to numerically simulating the gob material: elastic and plastic. 1,2 In the elastic approach, it is assumed that the gob material behaves as either a linear or nonlinear elastic material. For example, Jiang et al 2 and Cheng et al 4 used very elastic soft material (the assigned deformation modulus and Poisson's ratio were constants) to simulate gob material. Because of their simplicity, these models ignored the strain-hardening behavior of the compacted gob material. Saeedi et al 36 developed a numerical gob model based on the Salamon model (nonlinear elastic model) to capture the gob compaction and response. The formulation of this proposed method should be implemented by means of FISh (a FLAC3D programming language) in every calculation step. In the plastic approach, the double-yield (DY) model is commonly used to simulate the strain-stiffening behavior of materials, in which the support capacity increases as its volume is gradually compressed. Yavuz 37 commented that the DY model, which is available in FLAC3D, is capable of simulating the irreversible compaction of gob material in addition to shear yielding. Therefore, it is reasonable to employ the DY model to simulate the compaction of the gob material.
According to the available literature, the input parameters required for the DY model can be divided into two groups: the cap pressure and material properties. The cap pressure could be calculated using Salamon's model, which is expressed by the following equation: where σ v is the vertical stress applied to the gob material; ε v is the current vertical strain under the applied vertical stress; E 0 is the initial tangent modulus, E 0 = 10.37 1.042 c b 7.7 , ε m is the maximum volumetric strain of the gob material, m = (b − 1) ∕b; σ c is the compressive strength of the rock pieces; and b is the bulking factor of the caved rock, which is merely dependent on the mining height and caved zone height; b can be determined as follows: where h c is the mining height and h cr is the height of the caved zone, which is approximately 2-8 times the mining height. For LW 5302, the total mining height was 8.67 m, while the height of the caved zone was assumed to be 25.3 m. According to Equation (8), the bulking factor, maximum volumetric strain, and initial tangent modulus were 1.34, 0.25 and 18.34 GPa, respectively. Hence, the calculated cap pressure parameters for the DY model are listed in Table 3.

| Verification of the calibrated doubleyield model for gob modeling
Investigations have shown that it is difficult to estimate material properties, that is, the shear modulus, friction angle, and cohesion, for gob materials because these properties could not easily be obtained from field surveys or laboratory experiments. Therefore, a trial-and-error method was used to determine the unknown parameters through numerical model calibration against Salamon's model. For this purpose, a single-element model with dimensions of 1 m × 1 m × 1 m was developed to validate the material properties for the DY model. Both sides of the model were confined, and a constant velocity was applied on its top surface to model the vertical loading. Based on the DY model, a calibration procedure in line with that suggested by Yavuz 37 was followed to match the material properties needed as input parameters. Figure 4 shows the numerical and analytical stress-strain curves of the loading tests of gob samples. The final input parameters for the DY model are presented in Table 4. Using these parameters, it was possible to reasonably reproduce the strain-stiffening behavior of gob materials.
To evaluate the reliability of the calibrated input parameters for gob modeling, the vertical stress in the caved zone was calculated and plotted in Figure 5. The results show that the vertical stress at the gob was approximately zero and increased to 90.25% of the virgin stress approximately 300 m behind the face line. Past studies have reported that the vertical stress in the gob returns to its virgin stress level at a distance of approximately 0.2-0.3 times the overburden depth. The model results of gob modeling are in agreement with previous studies. The FLAC3D model also shows that the stress in the gob center was approximately 19.1 MPa, while the gob side showed a stress of 14.8 MPa after the longwall face was advanced. This difference was probably partly due to the side coal pillar acting as an additional support and altering the load on the gob side.

| Mesh dependency
As discussed in Section 3.4, the coal seam and roof/floor rock layers were modeled as a strain-softening material. Trueman and Mawdesley 38 suggested that numerical methods that use a strain-softening approach are not robust because the postpeak response is highly sensitive to the mesh size. Therefore, the simulation results might be sensitive to the mesh density of the numerical model, especially for the coal seam involved in the  To account for the mesh dependency, five models that only differed in the mesh size/mesh density used to the model the coal seam involved in LW 5302 mining as presented in Section 3.1 were also investigated. Figure 6 demonstrates the simulated principal stress monitored in the zone at the same depth as the coal seam at the last mining step with different mesh sizes/densities. The data in Figure 6 show that the gateroad deformation of the numerical model had a clear relationship with the mesh density. With a mesh density lower than 1 168 830, the simulated gateroad deformation increased quickly with mesh density. As the mesh density increased to 1 407 600, the gateroad deformation approached a constant level, and a further increase in the mesh density did not significantly improve the gateroad deformation. Because there was no marked difference between model 5 and model 4, the latter of which had a smaller element size, model 4 was considered to be the optimum mesh density for the model configuration to produce reasonably accurate simulation results.

| Validation with field monitoring
The coal/rock strength, gob model, and mesh dependency were calibrated, and they were then combined into the global model. To verify the global model, the simulated deformation of the gateroad was compared with field measurements.  Figure 7 shows the comparison of the LW 5302 tailgate deformation for the numerical results and field measurements at different stages of LW 5302 mining. Both the simulated results and field measurements increased as LW 5302 approached the monitoring location. The numerical simulation predicted final roof-to-floor and rib-to-rib convergences of 656 and 896 mm, respectively, while the measured roof-tofloor and rib-to-rib convergences were 710 and 924 mm, respectively. The numerical deformation agreed well with the field measurements. The maximum difference in the deformation between the model results and field measurements was 54 mm for roof-to-floor convergences and 28 mm for rib-to-rib convergences. Those differences might be due to an inherent defect of the continuum numerical model, by which it is difficult to capture fracture generation, propagation, and coalescence; these discontinuous fractures lead to lower deformation modules, which result in severe tailgate deformation in the field. However, for a large part, the numerical results shown in this work reproduced the "real" behaviors of the rock mass.

| MODEL RESULTS AND ANALYSIS
The strain energy parameters were defined in FLAC3D to capture the mining-induced energy distribution in the coal seam. The elastic strain energy is given by Miao et al 39 as follows: where W e is the elastic strain energy of the surrounding rock. In the present study, Equation (9) was implemented by means of FISH. In this manner, the mining-induced energy distribution can be illustrated.   Figure 8 shows that before coal excavation, the strain energy was uniformly distributed in the coal seam. Duo to the coal extraction, the strain energy density was redistributed and changed in a complex manner in front of the face. The strain energy was relatively low near the coal face and rapidly increased to a peak value in front of the coal face. The strain energy then decreased to the original energy level approximately 200 m in front of the coalface. The model results indicate that the energy distribution caused by the longwall mining was generally similar to the front abutment pressure distribution. The similarity was because the stress distribution in the coal seam had to be readjusted for a new state of equilibrium, which could drive the energy redistribution. Figure 9 shows the principal stress tensors in the roof above the mined-out area simulated in the numerical model. As seen in Figure 9, after the longwall was sufficiently excavated, a stress arch could be clearly identified, which transferred the overburden load to the unmined area, forming an overstressed zone. The most significant influence of this high stress was that it would cause the coal seam to absorb strain energy, leading to the formation of an energy concentration area.
The magnitude of the peak strain energy density (PSED) is presented in Figure 10 each time the coalface was advanced 12 m. Figure 10 shows that the PSED was approximately 1.18 and 2.19 times the premining energy level. Simultaneously, the location of the PSED reached 16-28 m ahead of the working face. As presented above, the computed energy distribution may be similar to the stress distribution. In past field measurements, the peak front abutment pressure was 2.0-2.5 times the initial field stress. In the current model, the ratio of the peak strain energy density to the premining energy was 2.19 times, which is in agreement with previous studies. As shown in Figure 10, as the face advanced, the PSED increased in magnitude while its location moved farther away from the face line. At the initial coalface advance, the magnitude of the PSED increased gradually. After approximately 204 m of advance, the magnitude increased to 176.7 kJ/m 3 and remained stable with a slight change. To display the dependence on the excavation scale, an approximate relation between the PSED and the advancing distance is expressed as: where E p1 is the PSED in the unmined coal seam adjacent to the opening, and A 1 , t 1 and y 0 are fitting constants. Studies show that the more the strain energy accumulated, the more severe the rockburst potential. Hence, the rockburst potential was more severe at the later stage than at the beginning, and more attention should be paid to the rockburst potential at the stable mining stage. Rockburst prevention techniques, such as destress blasting, water infusion, and directional hydraulic fracturing, might be carried out in advance to decrease the strain energy stored within the coal.

| Energy distribution at the barrier pillar along the strike
The barrier pillars protect the border area from disturbance by longwall mining. The energy distribution in the barrier pillar is associated with the mining stage. The strain energy distribution in the barrier pillar at different mining stages is displayed in Figure 11.
It is observed from these figures that the extraction of coal resulted in energy concentration in the barrier pillar. The strain energy in the barrier pillar adjacent to the 5302 headgate was relieved, which is described in blue in Figure 11. The area of the energy concentration increased with the face advance, while the area of energy relief remained almost unchanged. The changes were due to the transfer of overburden stress. The PSED along the strike reached 10-14 m ahead of the face line and was approximately 0.98-1.17 times the mining height. Clearly, the energy concentration zones in the barrier pillar were closer to the working face than those in the mined seam ahead of the face line. To analyze the energy evolution in detail, the PSED with face advance is plotted in Figure 12. As shown in Figure 12, after approximately 204 m of full extraction mining, the PSED in the barrier pillar increased to 1.88 times the initial energy level and then stabilized as the face advanced. A similar relation between the PSED and the advancing distance can also be expressed as: where E p2 is the PSED in the unmined coal seam adjacent to the opening and A 1 , t 1 and y 0 are fitting constants. As shown in Figure 13, the predicted strain energy density along the dip varied significantly in the coal seam ahead of the working face. A similar horse saddle-shaped horizontal energy distribution was predicted across the panel width, where the curve exhibited a "double peak" shape. The strain energy density in the panel's central region was approximately 1.14-1.38 times the original strain energy level. However, the peak energy region near the 5302 tailgate corresponded to approximately 1.60-2.78 times the original strain energy, while near the 5302 headgate, the corresponding value was approximately 1.21-2.12 times, which was lower compared to the 5302 tailgate. In past field observations, it has also been found that high-stress-induced guttering failure and rib spalling have always been common in the 5302 headgate. There may be complex reasons for this phenomenon, but the main reason is that the effect of side abutment stress due to the mined-out LW 5301 causes the stress to concentrate rapidly near the 5302  Closer to the location of the "double peak," when the section was not less than 125 m ahead of the working face, the location of the peak energy region moved further away from the gateroad rib with the face advancing and then stabilized. After 204 m of extraction, the peak energy region near the 5302 tailgate was located at a depth of approximately 26.0 m from the 5302 tailgate rib. However, near the 5302 headgate, the corresponding value was 9.16 m, which was also lower than that from the 5302 tailgate. Based on the model results, it can be deduced that in LW5302 of the Xinhe coal mine, the destress blasting and hydraulic fracturing that are used to relieve the stress or energy concentration should be applied at different designed depths. Figure 14 also illustrates the mining-induced energy distribution at the unmined coal seam along the dip. Figure 14 shows that the distribution of the strain energy density in the barrier pillar has a "single peak" shape. The strain energy density was approximately zero at the 5302 headgate rib. It increased rapidly and then decreased exponentially to its original value at a distance of approximately 0.22 times the overburden depth inward away from the ribs of the 5302 headgate. Further investigation into the magnitude and location of the PSED at the barrier pillar along the dip revealed that the corresponding value increased when the face advanced in the same section. For a section 25 m ahead of the working face, the PSED gradually increased and then stabilized at approximately 1.97 times the original strain energy level, occurring at approximately 8.4 m inward away from the ribs of the 5302 headgate, while for a section 175 m in front of the working face, the corresponding values were 1.09 times and 6.85 m, respectively. This indicated that the closer to the working face, the greater the rockburst risk. Figure 16 illustrates the magnitude and location of the PSED at the barrier pillar and coal seam in front of the working face.

| Energy distribution at the barrier pillar along the dip
As shown in Figures 13 and 14, the magnitude and location of the PSED at the barrier pillar along the dip were always lower than those at the coal seam in front of the working face. For example, the ratio of the PSED in mined seams near the 5302 headgate to that in barrier pillars was approximately 1.04-1.25. One reason is that they are caused by the front abutment pressure. The coal seam near the 5302 headgate is closer to the LW 5302 gob than the barrier pillar. The front abutment pressure will cause more stress concentration in the mined seam than in the barrier pillar. As a result, the energy concentration increased in the coal seam near the 5302 tailgate.

| Mining-induced energy distribution in the plane of the coal seam
Based on the FLAC3D modeling results, a conceptual model is constructed to demonstrate the energy distribution in the No. 3 coal seam in front of the working face (see Figure 15). The general results can be summarized as follows: 1. Both the coal seam and barrier pillar ahead of the working face can be divided into three regions along the strike. In front of the mining face along the strike, the strain energy density was approximately zero near the coal face and increased quickly with increasing distance from the coal face. At some distance ahead of the working face, the PSED in the coal seam along the strike was located at the boundary between region I and region II, while the PSED in the barrier pillar along the strike was located at the boundary between region IV and region V. Regions II and IV are the de-energized zones in the coal seam and barrier pillar. Because shearer cutting was implemented behind this region, the coal mass failed when the energy was released during mining, and its stored energy was lower than the original value. In regions I and V, because the stress exceeded the original stress, the coal mass absorbed more strain energy. Therefore, regions II and IV were the plastic-softening-strengthened zones. In regions III and IV, when the distance from the working face increased, the strain energy density decreased gradually and finally approached its original value U y , as shown in Figure 15. 2. Considering that the stress state in the coal seam ahead of the working face was seriously affected by the side abutment pressure of LW 5301, four cross-sections (cross-sections a-a, b-b, c-c, and d-d) in Figure 15 show the position of the PSED along the strike in the coal seam and barrier pillar. L 1 , L 2 , L 3 , and L 4 are the corresponding distances between the working face and the PSED. The model results showed that L 4 ≈ L 3 < L 2 < L 1 . This indicates that the closer the region was to the LW 5301 gob, the closer the PSED to the working face. With the face advancing, L 1 varied from 16 to 28 m; L 3 varied from 12 to 20 m; L 4 ≈ L 3 varied from 10 to 14 m. 3. The four blue dotted lines in Figure 15 show the position of the PSED along the dip in the coal seam and barrier pillar.  4 . In addition, as for the effect of the excavation scale, the growth of the PSED in coal seams is likely to be more severe than that in barrier pillars.

| The effect of the extracted coal panel on the energy distribution
The closer the region was to the LW 5301 gob, the larger the magnitude of the strain energy density. The energy distribution in the coal seam and barrier pillar was affected by the interaction between the front abutment pressure (formed by LW 5302 extraction) and side abutment pressure (formed by LW 5301 extraction). This interaction can be interpreted using the bearing beam analogue, which has been generally accepted as a simplified tool for the analysis of overburden load transfer during underground space development. [40][41][42] Gao et al 43 stated that the strata located in the lower portion of the continuous deformation zone can be considered a bearing beam. This bearing beam complements the "stress arch" developed in the roof above the unmined coal, forming an entire compressive stress arch around the goaf.
Diederichs and Kaiser 44 also suggested that this compressive arch transfers the weight of overlying strata to the abutments. As a result, abutment pressures are formed around the edges of the gob. Figure 16 shows the principal stress tensors in the laminated roofs above the mined-out area simulated in the numerical model after sufficient extraction of LW 5301. It can be seen that a compressive arch, which causes the side abutment pressure to act along the side of LW 5301, can be clearly identified in the cross-section along the dip. As demonstrated by previous studies, the magnitude of the side abutment pressure changes in the influenced zone of the side abutment. Assuming that extraction of LW 5301 had little influence on the movement of overlying strata, a continuous deformation zone above LW 5302 gob was also formed in all cross-sections along the direction of LW 5302 mining, and a compressive stress arch could be formed. Its skewback lay in the coal seam (see Figure 16) and formed the front abutment pressure ahead of the LW 5302. Ahead of the face, a "stress-strengthening" effect occurred in the coal seam where the skewbacks of the side and front abutment pressure intersected and were superimposed on each other. Due to the nonuniform side abutment pressure along the dip, this "stress-strengthening" effect seemed to decrease. As a result, the farther outward away from the LW 5301 goaf, the smaller the stress. Considering a unit zone within the infinite in situ rock mass, any stress increase within the zone must be accompanied by energy stored within the zone. Therefore, the closer the region was to the LW 5301 gob, the higher the strain energy density.
It would be easy to anticipate that the nonuniform distribution of energy in the coal seam adjacent to the extracted coal panel generated by the overburden movement in the whole goaf area is difficult to change due to geotechnique limitations, once certain geological and mining conditions are determined. As a principal part of the energy storage unit, the energy storage capacity of the coal seam can be reduced by destroying its structural integrity through certain special geotechniques. When the coal seam was pre-fractured ahead of the coal face along the interedge of the gateroad, the width of the elastic core could be reduced, and the whole coal seam even yielded or failed, causing the overburden load to be transferred into neighboring mine structures.

| The effect of the extraction scale on the energy distribution
As discussed in Sections 4.1 and 4.2, it is easy to note that with increasing extraction scale, the PSED increased until the advancing distance reached the sum of the widths of LW 5301 and 5302. It is well known that the mining-induced stress or energy distribution depends on the strata strength, stiffness, overburden depth, retreating rate, and extraction scale (advancing distance, mining height, and longwall width). 3,8,45 Once certain geological and mining parameters of the coal face are known, including the overburden depth, retreating rate, mining height, and longwall width, it is easy to understand the relationship between the Gob Headgate Tailgate L 4

F I G U R E 1 6
The principal stress tensors in the laminated roofs above the mined-out area mining-induced stress or energy distribution and the advancing distance. Based on previous field investigations and numerical modeling, Xie et al 46 stated that the voussoir beam formed by the main roof lies in the fractured zone, forming a "fracture arch" with the fracture positions of each critical stratum as the boundary. In Chinese underground coal mining practices, 46 the height of the "fracture arch" increases gradually at the beginning of the mining process and reaches a maximum at approximately half of the longwall face width when the longwall face advances the width of the longwall face. As described in Section 2.1, LW 5302 employs the single-entry system, with a small pillar left between the adjacent LW 5301 gob and LW 5302. After LW 5302 had advanced from the setup room, the failure of the yield pillar occurred in the gob, and the gob of LW 5301 and LW 5302 had amalgamated into a vast gob, with a width of 190 m (the sum of the two longwall face widths). Based on the interpretation above, once the LW 5302 had advanced approximately 190 m, the height of the "fracture arch" could reach its maximum, and the height of the "stress arch" was also unchanged. As a result, the "overburden load transfer" effect seemed to stabilize. Thus, it could be deduced from the model results that the PSED gradually increased and then stabilized as the face advanced by approximately 204 m.

| CONCLUSIONS
The aim of this case study conducted at the Xinhe coal mine, Shandong province, China, was to investigate the mininginduced energy redistribution based on numerical modeling. The main conclusions presented in this paper can be summarized as follows: 1. Numerical modeling using FLAC3D software was used to evaluate the mining-induced energy change level. To improve the appropriate representation of the gob-compaction process when modeling longwall mining, a double-yield model was used to model the gob material response, and its input parameters were meticulously calibrated based on a back-analysis procedure. Furthermore, this paper also provided a basic set of calibration procedures (including mesh dependency and validation of the global model) that can be used to evaluate the energy, stress and deformation associated with longwall mining. 2. The model results predict that in front of the longwall face, there is a range (approximately 200 m) that is seriously affected by longwall mining. In this range, the strain energy is relatively low at the coalface and rapidly increases to a peak value. In the coal seam, the peak strain energy density was approximately 16-28 m ahead of the coalface, while in the barrier pillar, the peak strain energy density was located approximately 10-14 m in front of the coalface. The closer the region was to the LW 5301 gob, the higher the location and magnitude of the peak strain energy density. As a result, the risk of rockburst in coal seams is greater than that in barrier pillars. In addition, with increasing advancing distance, the magnitude of the peak strain energy density gradually increased and then stabilized where the coalface had advanced approximately 200 m. This result indicated that the rockburst potential was more severe at the later stage than at the beginning of longwall face mining. 3. The form of the energy distribution along the dip is almost the same when the distance in front of the coal face is different. The difference is the magnitude of the strain energy density. In the coal seam, the energy distribution has a similar horse-saddle shape, where the curve exhibits a "double peak" shape. In the barrier pillar, the energy distribution has a "single peak" shape. The peak strain energy density in the coal seam was approximately 16-28 m inward away from the ribs of the 5302 tailgate, and in the barrier pillar, the peak strain energy density occurred approximately 8.4 m inward away from the ribs of the 5302 headgate. As a result, the destress blasting and hydraulic fracturing that are used to relieve the stress or energy concentration should be applied at different designed depths.