# Understanding the Distribution Laws of In Situ Stress in Mines Based on Global Positioning System Data and Numerical Simulation Analysis

Authors: Wentao Liu; Zhongyou Zheng; Lei Zhu; Kai Xu; Chengyong Liu; Yuyi Wu; Shunzhen Hu
Journal: Geofluids
 (2023)

---

## Abstract

To ensure safe exploitation of the Muduchaideng Mine, we use a numerical simulation to invert the in situ stress field based on observation data from GPS stations, investigate differences in movement velocities at different locations, and analyze the lefthanded movement in the area around the wellfield. The relationship between the difference in motion velocity and the in situ stress field was explained through theoretical analysis. The law of in situ stress distribution was studied using data collected from the surrounding mines, and a regression equation as well as the relationship between the magnitudes of the three-way stresses was obtained. The inversion results were consistent with the measured data, indicating that inversion of the in situ stress field was accurate under the conditions found in the intact Ordos block. Additionally, the influence of tectonic stress on the key layer was larger than that on the bedrock of the coal seam roof, and the energy accumulated in the key layer was higher. The stress concentration coefficient of the key layer was similar to the variation law of the coal seam thickness. Therefore, the southwest area of the wellfield with a thicker coal seam is the key area for impact ground pressure prevention and control.


---

## Body

## 1. Introduction

The Muduchaideng Mine is located southeast of the Hujierte mining area in the central part of the Ordos block, which is surrounded by fault structures and fault basins. Presently, this area experiences frequent tectonic activity and dynamic tectonic stress. Many scholars have comprehensively studied the motion characteristics and in situ stress distribution law in this block [1, 2]. Li et al. analyzed the crust-mantle kinematic characteristics of the Ordos block and its surrounding areas using global positioning data (GPS), SKS shear-wave splitting results, and other geophysical data. They found that the Ordos block rotates counterclockwise relative to the Eurasian continent and that its internal structure is relatively intact; moreover, the western and eastern edges of the block show strong movement, while the southern and northern edges exhibit weak movement. Additionally, the boundary fault shows sinistral movement [3]. Li et al. analyzed the relationship between coal seams of different burial depths and in situ stresses and studied the distribution of the in situ stress field and its control effect on fractures of coal reservoirs in the Liulin area of the Ordos Basin. Furthermore, by combining the regional structural characteristics and the testing parameters of the CBM injection pressure drop, information on the evolution of the in situ stress field and its rule of change in the plane and vertical directions was obtained [4]. Liu et al. accurately predicted the spatial distribution of stress fields in heterogeneous continental reservoirs with a reservoir mechanics heterogeneity model using the Ordos Basin, China, as an example [5]. Similar studies have been conducted in other regions for the in situ stress field [6, 7]. Ganguli and Sen developed a comprehensive geomechanical model to evaluate the present-day stresses and pore pressures in a prolific onshore hydrocarbon-bearing field located in the southern Cambay Basin, western India [8]. Ma et al. proposed a hybrid neural network to predict horizontal in situ stresses based on logging data from X-1 and X-2 horizontal wells in the Sichuan Basin, southwestern China [9].Most of the mines in the vicinity of the Muduchaideng Mine have been newly constructed and just been opened. Various geological aspects remain to be improved, and studies on the in situ stress in the area have rarely been reported. Due to the complex arrangement of underground tunnels and significant stress concentration, implementation of in situ measurements of ground stress and their inversion on a large scale are difficult. Therefore, numerical simulation analysis based on more convenient motion measurement methods is more effective in studying the distribution pattern of in situ stress. This paper is organized as follows: Section2 describes the geological environment around the well field and the basic conditions of the Muduchaideng Mine. Section 3 focuses on GPS data analysis and introduces the underlying analysis conditions and data sources. Section 4 presents results of the inversion of the in situ stress field based on numerical simulation software. In Section 5, a theoretical analysis of the tectonic stress is presented. In Section 6, the measured data are applied to verify the simulation results, and Section 7 presents the conclusions.


## 2. Study Area

The general outline of the Ordos block is similar to a rectangle, with its tectonic form resembling a broad and gentle large oblique. The core of the oblique is to the west, and it is characterized by a broad and gentle eastern flank and a relatively steep western flank. The geological structures around the block are complex, with the Alashan, Qilain, and Qaidamu blocks to the west, the Yinchuan-Hetao rift system in the north, and the Shaanxi-Shanxi rift system and Datong volcanic group in the east (Figure1). However, the geological structures within the block are simple, without large faults and folds. The central and eastern areas consist of subhorizontal rock formations.Figure 1
Ordos block periphery tectonic distribution map.The Muduchaideng Mine is located in the central part of the Ordos block. It has a 6.0 Mt/a production capacity and includes 3-1 coal as the main product, which is buried at a depth of 617.7–869.6 m, with an average depth of 662.5 m. The depth increases from southeast to northwest, with the maximum depth observed in the recessed area of the northwest region. The thickness of the coal seam is 3.40–7.63 m, with an average value of 5.06 m. Its overall structure is monocline with a dip angle of 1°–6°, and it stretches out in a NNE direction. The coal seam roof is mainly composed of siltstone and fine-grained sandstone, and the lithology of the floor is mainly mudstone, sandy mudstone, and siltstone. The 3-1 coal is divided into two mining areas by the industrial site at the center of the wellfield: 301 mining area on the east and 302 mining area on the west. A main roadway in the east–west direction along the center of the wellfield divides the two mining areas into two parts each: a northern and a southern part. As the coal seam is buried deep, the underground stress concentration is relatively significant. Consequently, rock bursts have occurred several times in the adjacent mine, leading to roadway blasting and a collapse of the roof, which seriously affected mine production safety. There have been many cases of top and bottom slab settlement of more than 1 m, damage, deformation of the support structure, etc. Being a rock burst mine, rock burst prevention and control in the Muduchaideng Mine must therefore be planned and implemented.


## 3. GPS Data Analysis

Since the late Quaternary (before 100000–120000 a), when the active tectonic zone divided the Ordos block and formed a geological unit with a unified motion, the sliding velocity of the structural belt has maintained the same speed and mode of motion; that is, the fault slip rate remained unchanged for a long time [10, 11]. Therefore, the velocity boundary conditions of the block surrounded by active faults have remained constant since the late Quaternary, and the movements of all points inside the block have been relatively stable, ensuring the integrity of the internal block structure.The Ordos block is surrounded by fault zones and has been widely recognized as a simple integrated geological structure with weak geological activity [12, 13]. Thus, in general, GPS data can easily be used to represent the approximate mean motion velocity of the long-term observation point.Presently, the relative positioning accuracy of high-precision GPS is 10-8–10-9 [14]; that is, the positioning error for the baseline observation within a range of ~10 km of the Muduchaideng wellfield is ~10-2–10-3 mm, which can accurately describe the differences in the movement inside the wellfield. It provides an accurate and efficient method for surface surveys in conditions where the surface of the well field is deserted and the conditions are relatively vicious.Wang et al. collected 362 observations from GPS stations, which mainly depicted the velocity field of tectonic deformation within the Chinese continent [13]. Of these observations, we selected data from around the Muduchaideng wellfield and interpolated it using the Kriging method so that we could extract the velocity distribution pattern of the target area. Subsequently, the geographic coordinates were converted into plane rectangular coordinates according to the corresponding relations of flex points at the wellfield boundary, and the velocity distribution isograms along the eastern and northern directions within the wellfield range were constructed (Figures 2 and 3).Figure 2
Velocity distribution isogram along the eastern direction within the wellfield (mm/a).Figure 3
Velocity distribution isogram along the northern direction within the wellfield (mm/a).The eastward velocity isolines in the wellfield were parallel lines spaced equidistantly (Figure2). The moving velocity increased from the northeast to the southwest of the wellfield, and the velocity gradient direction was ~S15.61°W. Furthermore, the northward velocity curves in the wellfield (a negative value indicates southward velocity) were parallel and spaced equidistantly (Figure 3). The moving velocity increased from northwest to southeast of the wellfield, and the velocity gradient direction was ~S18.43°E–S35.00°E. These isoline curves gradually turned southward, and the movement characteristics of the interior layer within the wellfield were consistent with those previously observed in the Ordos block [15].


## 4. 3D Numerical Simulation

### 4.1. Judgment of the Key Stratum

Presently, the second working face of the 302 mining area is being mined. Borehole B22 is located approximately 280 m along the tendency to the east of the production working face and 1182 m along the strike distance from the open-off cut. The depth of the drill hole is 847.74 m, in which the burial depth of the 3-1 coal seam is 630.79 m. Borehole B22 was used to determine the key stratum according to the mechanical parameters of all strata in the nearby area, as it is closest to the mining position. The key layer can be determined in two steps: (i) Determine the location of the hard rock layer based on the deflection size of the rock deformation; (ii) calculate the breakage step of each hard rock layer and determine the sequence of rock layer breakage. The layer with the largest breakage step is the main key layer [16, 17]. The corresponding results are shown in Table 1.Table 1
Borehole B22 key layer discrimination results.
NumberStrataLithologyThickness (m)Distance to coal seam (m)Discriminatory results1QhWind-accumulated sand130.38500.412K1zSiltstone65.91434.5Key layer3Fine-grained sandstone29.75404.754Medium sandstone15.98388.775Fine-grained sandstone8.84379.936Silty mudstone14.05365.887Fine-grained sandstone3.3362.588Coarse-grained sandstone30.34332.24Key layer9Mudstone5.5326.7410Fine-grained sandstone16.31310.4311J2aCoarse-grained sandstone45.03265.412Silty mudstone40.05225.3513Medium sandstone7.2218.1514J2zSilty mudstone57.94160.2115Medium sandstone6.63153.5816Mudstone1.4152.1817Medium sandstone12.17140.0118Silty mudstone41.1298.8919J1−2y3aCoarse-grained sandstone61.5737.32Key layer20Fine-grained sandstone1.735.6221Mudstone0.435.2222Coal1.333.9223Fine-grained sandstone20.7613.16Key layer24Medium sandstone1.311.8625Fine-grained sandstone11.86026Coal5.327Silty mudstone8.88The results showed that the coal seam roof contained four key strata: fine-grained sandstone stratum (20.76 m thick), coarse-grained sandstone stratum (61.57 m thick), coarse-grained sandstone stratum (30.34 m thick), and siltstone stratum (65.91 m thick). These four strata were located 13.16 m, 37.32 m, 332.24 m, and 434.5 m away from the coal seam roof and were defined as low-location, middle-location, middle/high-location, and high-location key strata, respectively. According to the on-site data statistics, the middle-location key stratum occupied the highest ratio of total energy release of the roof, had the most serious impact on the rock burst of the working plane, and had the most stable occurrence. Therefore, in this study, the middle-location key stratum (hereinafter referred to as the key stratum) was considered for simulation and analysis.Based on the data of the 98 borehole columns in the wellfield, strata positions, and thicknesses of the overlying stratum of the 3-1 coal, the thickness of the topsoil stratum was determined to be 22.8–159.0 m, with an average of 103.5 m, while that of the key stratum was 11.6–85.3 m, with an average of 44.0 m. Most boreholes were interpreted as medium-grained sandstone and coarse-grained sandstone, while some were interpreted as fine-grained sandstone and siltstone. The spacing between the bottom of the key stratum and the roof of the 3-1 coal seam was 9.0–92.1 m, with an average of 43.1 m.


### 4.2. Model Development

To strengthen the influence of the key strata and weaken the influence of other strata, the rock strata in the borehole column were regrouped and sequentially divided from top to bottom into topsoil stratum, bedrock stratum 2, key stratum, bedrock stratum 1, coal seam, and bottom. After unifying the data density using software interpolation [18, 19], the interface of the rock strata was drawn, and subsequently, the FLAC3D software was imported to divide the cells. Consequently, 297,600 cells were divided. To reduce the influence of the boundary conditions on the simulation results, the model was expanded by 1 km around the wellfield boundary. The dimensions of the established model were 10000m×12400m (totallength×totalwidth) and 400–1299.18 m height.


### 4.3. Simulation Process and Result Analysis

Tectonic stress is formed gradually during long-term crustal movement. When the tectonic belt at the block boundary changes, the internal stress distribution is disturbed and reaches a new equilibrium.To avoid plastic damage and stress concentration caused by cell failure, the model was developed based on an elastic constitutive model with initial equilibrium under the influence of self-weight stress. Since the wellfield is in the middle of the Ordos block and the surrounding strata are intact structures, the boundary conditions have less influence on the calculation results. Therefore, the boundary conditions were as follows: the model edges and bottom were fixed in the normal direction, and the upper surface was a free surface. Its mechanical parameters are shown in Table2.Table 2
Rock mechanics parameters.
NumberLithologyElastic modulus of rock, GPaPoisson’s ratio of rock1Baseboard9.080.32Coal2.680.263Bedrock30.120.34Key layer39.520.35Loose bed0.10.2


#### 4.3.1. Initial Equilibrium of the Model

After the initial equilibrium, the vertical stress was the maximum principal stress in the model, and the intensity of the maximum principal stress in the key stratum was 13.4–23.9 MPa (Figure4). The distribution law indicated that the stress increased gradually from the southeast to northwest of the wellfield, with the highest stress near the recessed area of the northwestern part and low stresses in the southern and eastern parts. This was consistent with the change law of the burial depth of the coal seam.Figure 4
Maximum principal stress distribution of key strata after the initial equilibrium (Pa).


#### 4.3.2. GPS Data Import

After the model reached equilibrium, the analyzed GPS data were discretized and projected onto the model. Subsequently, the velocity vector was loaded into each cell by the built-in language programming application in the numerical simulation software, and the simulated velocity amount and direction in the study area were verified, which confirmed that they were highly consistent with the measured data. Additionally, in line with the crustal movement law, the boundary constraints around the model were removed, and the boundary conditions were changed to the measured velocity data. After importing the data, the east and north velocities were calculated as shown in Figures5 and 6.Figure 5
East velocity isogram in the model after applying tectonic stress (m/a).Figure 6
North velocity isogram in the model after applying tectonic stress (m/a).


#### 4.3.3. Tectonic Stress Loading

As the block boundary conditions in different geological periods are unclear, the following approximations were made for block movement in the simulation:(i)
The late Quaternary (~100000–120000 a, considering its average of 110000 a) was considered the loading time of tectonic stress(ii)
During the initial stage of tectonic stress loading, the tectonic stress accumulated in the previous period was fully released(iii)
The calculation time step corresponded to the movement period of the blockAfter considering these approximations, the movement velocity of cells was reset before each operation to ensure uniform movement velocity and boundary conditions inside the block. The distributions of the maximum principal stress and movement velocity of the key stratum after the simulation are shown in Figure7. The distribution of the maximum principal stress direction of the key stratum is shown in Figure 8.Figure 7
Isogram of maximum principal stress and moving velocity of the key stratum after loading tectonic stress (Pa, m/a).Figure 8
Isogram of the southeast angle of the maximum principal stress in the key stratum after loading tectonic stress (°).The maximum principal stress was 21.3–26.4 MPa after applying tectonic stress. The maximum principal stresses in the northwestern, recessed, and southwestern parts of the wellfield were comparatively higher than those near the eastern boundary. Compared to the stress distribution after the initial equilibrium, the area under low maximum principal stress decreased after the tectonic stress was applied, while that under moderate principal stress increased. In general, these changes were more and less prominent in the western and eastern areas, respectively, whereas no significant changes were observed in the central areas. The maximum principal stress direction is S60.07°E-S64.48°E, with a slightly larger offset angle in the southwestern part of the wellfield and a gradually decreasing angle toward the northeastern part.The strata within the wellfield range moved from northwest to southeast in a direction of S72.16°E. The velocity range was 9.40–9.57 mm/a, decreasing from the southwest to northeast of the wellfield, thus suggesting that the wellfield gradually rotated counterclockwise.


#### 4.3.4. Calculation of the Stress Concentration Coefficient

To assess the impact risk, Dou et al. proposed the use of a relative stress concentration coefficient composition that considers the ratio of the maximum principal and gravitational stresses caused by a certain factor as the relative stress concentration coefficients since they affect rock burst in a certain target area [20]. The coefficient was used to evaluate the influences of the factor on rock burst.The maximum principal stress at the center point of each cell in the key stratum of the model was exported, and the quotient of the maximum principal stress at the initial equilibrium was calculated to obtain the relative stress concentration coefficient (hereinafter referred to as the stress concentration coefficient). Its distribution isoline is shown in Figure9.Figure 9
Distribution layout of the stress concentration coefficient of the key stratum.The distribution of the stress concentration coefficient differed greatly from that of the maximum principal stress in the key stratum but was similar to the variation law of coal seam thickness to some extent (Figure9). In most areas of the region, the stress concentration coefficient gradually increased from the northeast to the southwest of the wellfield. Generally, the mining height is limited to a certain range as the thickness of the coal seam changes; thus, the greater the thickness of the coal seam is, the higher the mining height. Accordingly, the coal seam in the southwest of the wellfield was found to be a key area for the prevention and control of rock bursts due to its large thickness, high stress concentration, large postmining disturbance range, and high energy release.To compare the differences in tectonic stress in different regions and strata, four measuring points were established in the model: No. 1 was located in the central part of the northwestern recessed area, Nos. 2 and 4 were located in the same monitoring position (different strata) in the center of the wellfield, and No. 3 was located near the eastern boundary of the wellfield. The measuring points for Nos. 1–3 were located in the key stratum, while the measuring point for No. 4 was located in the bedrock of the coal seam roof. The variations in the maximum principal and intermediate principal stresses in the simulation process are shown in Figure10.Figure 10
Stress variation curves of the measuring points.As seen in Figure10, after application of the tectonic stress, the intermediate principal stress, which was oriented in the same direction as the resultant force, increased rapidly. Consequently, it became the maximum principal stress by replacing the vertical stress, while the vertical stress remained unchanged. This indicated that the tectonic stress did not influence the vertical stress but significantly influenced the horizontal stress. Moreover, after the tectonic stress was applied, the increase rate of the horizontal stress in the key stratum was significantly higher than that in the bedrock of the coal seam roof, indicating that the tectonic stress greatly affected the key stratum and that higher energy is stored in the key stratum.Furthermore, the recessed area in the northwestern part of the wellfield exhibited the largest maximum principal stress but smallest stress concentration coefficient. This was primarily due to extensive variation in vertical and horizontal stresses near the recessed area. As a result, the horizontal stress required a longer time to reach and surpass the vertical stress and become the maximum principal stress. Finally, at the end of the simulation, the difference between the two stresses at the recessed area was significantly lower than at other positions, thereby resulting in the smallest stress concentration coefficient.


## 4.1. Judgment of the Key Stratum

Presently, the second working face of the 302 mining area is being mined. Borehole B22 is located approximately 280 m along the tendency to the east of the production working face and 1182 m along the strike distance from the open-off cut. The depth of the drill hole is 847.74 m, in which the burial depth of the 3-1 coal seam is 630.79 m. Borehole B22 was used to determine the key stratum according to the mechanical parameters of all strata in the nearby area, as it is closest to the mining position. The key layer can be determined in two steps: (i) Determine the location of the hard rock layer based on the deflection size of the rock deformation; (ii) calculate the breakage step of each hard rock layer and determine the sequence of rock layer breakage. The layer with the largest breakage step is the main key layer [16, 17]. The corresponding results are shown in Table 1.Table 1
Borehole B22 key layer discrimination results.
NumberStrataLithologyThickness (m)Distance to coal seam (m)Discriminatory results1QhWind-accumulated sand130.38500.412K1zSiltstone65.91434.5Key layer3Fine-grained sandstone29.75404.754Medium sandstone15.98388.775Fine-grained sandstone8.84379.936Silty mudstone14.05365.887Fine-grained sandstone3.3362.588Coarse-grained sandstone30.34332.24Key layer9Mudstone5.5326.7410Fine-grained sandstone16.31310.4311J2aCoarse-grained sandstone45.03265.412Silty mudstone40.05225.3513Medium sandstone7.2218.1514J2zSilty mudstone57.94160.2115Medium sandstone6.63153.5816Mudstone1.4152.1817Medium sandstone12.17140.0118Silty mudstone41.1298.8919J1−2y3aCoarse-grained sandstone61.5737.32Key layer20Fine-grained sandstone1.735.6221Mudstone0.435.2222Coal1.333.9223Fine-grained sandstone20.7613.16Key layer24Medium sandstone1.311.8625Fine-grained sandstone11.86026Coal5.327Silty mudstone8.88The results showed that the coal seam roof contained four key strata: fine-grained sandstone stratum (20.76 m thick), coarse-grained sandstone stratum (61.57 m thick), coarse-grained sandstone stratum (30.34 m thick), and siltstone stratum (65.91 m thick). These four strata were located 13.16 m, 37.32 m, 332.24 m, and 434.5 m away from the coal seam roof and were defined as low-location, middle-location, middle/high-location, and high-location key strata, respectively. According to the on-site data statistics, the middle-location key stratum occupied the highest ratio of total energy release of the roof, had the most serious impact on the rock burst of the working plane, and had the most stable occurrence. Therefore, in this study, the middle-location key stratum (hereinafter referred to as the key stratum) was considered for simulation and analysis.Based on the data of the 98 borehole columns in the wellfield, strata positions, and thicknesses of the overlying stratum of the 3-1 coal, the thickness of the topsoil stratum was determined to be 22.8–159.0 m, with an average of 103.5 m, while that of the key stratum was 11.6–85.3 m, with an average of 44.0 m. Most boreholes were interpreted as medium-grained sandstone and coarse-grained sandstone, while some were interpreted as fine-grained sandstone and siltstone. The spacing between the bottom of the key stratum and the roof of the 3-1 coal seam was 9.0–92.1 m, with an average of 43.1 m.


## 4.2. Model Development

To strengthen the influence of the key strata and weaken the influence of other strata, the rock strata in the borehole column were regrouped and sequentially divided from top to bottom into topsoil stratum, bedrock stratum 2, key stratum, bedrock stratum 1, coal seam, and bottom. After unifying the data density using software interpolation [18, 19], the interface of the rock strata was drawn, and subsequently, the FLAC3D software was imported to divide the cells. Consequently, 297,600 cells were divided. To reduce the influence of the boundary conditions on the simulation results, the model was expanded by 1 km around the wellfield boundary. The dimensions of the established model were 10000m×12400m (totallength×totalwidth) and 400–1299.18 m height.


## 4.3. Simulation Process and Result Analysis

Tectonic stress is formed gradually during long-term crustal movement. When the tectonic belt at the block boundary changes, the internal stress distribution is disturbed and reaches a new equilibrium.To avoid plastic damage and stress concentration caused by cell failure, the model was developed based on an elastic constitutive model with initial equilibrium under the influence of self-weight stress. Since the wellfield is in the middle of the Ordos block and the surrounding strata are intact structures, the boundary conditions have less influence on the calculation results. Therefore, the boundary conditions were as follows: the model edges and bottom were fixed in the normal direction, and the upper surface was a free surface. Its mechanical parameters are shown in Table2.Table 2
Rock mechanics parameters.
NumberLithologyElastic modulus of rock, GPaPoisson’s ratio of rock1Baseboard9.080.32Coal2.680.263Bedrock30.120.34Key layer39.520.35Loose bed0.10.2


### 4.3.1. Initial Equilibrium of the Model

After the initial equilibrium, the vertical stress was the maximum principal stress in the model, and the intensity of the maximum principal stress in the key stratum was 13.4–23.9 MPa (Figure4). The distribution law indicated that the stress increased gradually from the southeast to northwest of the wellfield, with the highest stress near the recessed area of the northwestern part and low stresses in the southern and eastern parts. This was consistent with the change law of the burial depth of the coal seam.Figure 4
Maximum principal stress distribution of key strata after the initial equilibrium (Pa).


### 4.3.2. GPS Data Import

After the model reached equilibrium, the analyzed GPS data were discretized and projected onto the model. Subsequently, the velocity vector was loaded into each cell by the built-in language programming application in the numerical simulation software, and the simulated velocity amount and direction in the study area were verified, which confirmed that they were highly consistent with the measured data. Additionally, in line with the crustal movement law, the boundary constraints around the model were removed, and the boundary conditions were changed to the measured velocity data. After importing the data, the east and north velocities were calculated as shown in Figures5 and 6.Figure 5
East velocity isogram in the model after applying tectonic stress (m/a).Figure 6
North velocity isogram in the model after applying tectonic stress (m/a).


### 4.3.3. Tectonic Stress Loading

As the block boundary conditions in different geological periods are unclear, the following approximations were made for block movement in the simulation:(i)
The late Quaternary (~100000–120000 a, considering its average of 110000 a) was considered the loading time of tectonic stress(ii)
During the initial stage of tectonic stress loading, the tectonic stress accumulated in the previous period was fully released(iii)
The calculation time step corresponded to the movement period of the blockAfter considering these approximations, the movement velocity of cells was reset before each operation to ensure uniform movement velocity and boundary conditions inside the block. The distributions of the maximum principal stress and movement velocity of the key stratum after the simulation are shown in Figure7. The distribution of the maximum principal stress direction of the key stratum is shown in Figure 8.Figure 7
Isogram of maximum principal stress and moving velocity of the key stratum after loading tectonic stress (Pa, m/a).Figure 8
Isogram of the southeast angle of the maximum principal stress in the key stratum after loading tectonic stress (°).The maximum principal stress was 21.3–26.4 MPa after applying tectonic stress. The maximum principal stresses in the northwestern, recessed, and southwestern parts of the wellfield were comparatively higher than those near the eastern boundary. Compared to the stress distribution after the initial equilibrium, the area under low maximum principal stress decreased after the tectonic stress was applied, while that under moderate principal stress increased. In general, these changes were more and less prominent in the western and eastern areas, respectively, whereas no significant changes were observed in the central areas. The maximum principal stress direction is S60.07°E-S64.48°E, with a slightly larger offset angle in the southwestern part of the wellfield and a gradually decreasing angle toward the northeastern part.The strata within the wellfield range moved from northwest to southeast in a direction of S72.16°E. The velocity range was 9.40–9.57 mm/a, decreasing from the southwest to northeast of the wellfield, thus suggesting that the wellfield gradually rotated counterclockwise.


### 4.3.4. Calculation of the Stress Concentration Coefficient

To assess the impact risk, Dou et al. proposed the use of a relative stress concentration coefficient composition that considers the ratio of the maximum principal and gravitational stresses caused by a certain factor as the relative stress concentration coefficients since they affect rock burst in a certain target area [20]. The coefficient was used to evaluate the influences of the factor on rock burst.The maximum principal stress at the center point of each cell in the key stratum of the model was exported, and the quotient of the maximum principal stress at the initial equilibrium was calculated to obtain the relative stress concentration coefficient (hereinafter referred to as the stress concentration coefficient). Its distribution isoline is shown in Figure9.Figure 9
Distribution layout of the stress concentration coefficient of the key stratum.The distribution of the stress concentration coefficient differed greatly from that of the maximum principal stress in the key stratum but was similar to the variation law of coal seam thickness to some extent (Figure9). In most areas of the region, the stress concentration coefficient gradually increased from the northeast to the southwest of the wellfield. Generally, the mining height is limited to a certain range as the thickness of the coal seam changes; thus, the greater the thickness of the coal seam is, the higher the mining height. Accordingly, the coal seam in the southwest of the wellfield was found to be a key area for the prevention and control of rock bursts due to its large thickness, high stress concentration, large postmining disturbance range, and high energy release.To compare the differences in tectonic stress in different regions and strata, four measuring points were established in the model: No. 1 was located in the central part of the northwestern recessed area, Nos. 2 and 4 were located in the same monitoring position (different strata) in the center of the wellfield, and No. 3 was located near the eastern boundary of the wellfield. The measuring points for Nos. 1–3 were located in the key stratum, while the measuring point for No. 4 was located in the bedrock of the coal seam roof. The variations in the maximum principal and intermediate principal stresses in the simulation process are shown in Figure10.Figure 10
Stress variation curves of the measuring points.As seen in Figure10, after application of the tectonic stress, the intermediate principal stress, which was oriented in the same direction as the resultant force, increased rapidly. Consequently, it became the maximum principal stress by replacing the vertical stress, while the vertical stress remained unchanged. This indicated that the tectonic stress did not influence the vertical stress but significantly influenced the horizontal stress. Moreover, after the tectonic stress was applied, the increase rate of the horizontal stress in the key stratum was significantly higher than that in the bedrock of the coal seam roof, indicating that the tectonic stress greatly affected the key stratum and that higher energy is stored in the key stratum.Furthermore, the recessed area in the northwestern part of the wellfield exhibited the largest maximum principal stress but smallest stress concentration coefficient. This was primarily due to extensive variation in vertical and horizontal stresses near the recessed area. As a result, the horizontal stress required a longer time to reach and surpass the vertical stress and become the maximum principal stress. Finally, at the end of the simulation, the difference between the two stresses at the recessed area was significantly lower than at other positions, thereby resulting in the smallest stress concentration coefficient.


## 4.3.1. Initial Equilibrium of the Model

After the initial equilibrium, the vertical stress was the maximum principal stress in the model, and the intensity of the maximum principal stress in the key stratum was 13.4–23.9 MPa (Figure4). The distribution law indicated that the stress increased gradually from the southeast to northwest of the wellfield, with the highest stress near the recessed area of the northwestern part and low stresses in the southern and eastern parts. This was consistent with the change law of the burial depth of the coal seam.Figure 4
Maximum principal stress distribution of key strata after the initial equilibrium (Pa).


## 4.3.2. GPS Data Import

After the model reached equilibrium, the analyzed GPS data were discretized and projected onto the model. Subsequently, the velocity vector was loaded into each cell by the built-in language programming application in the numerical simulation software, and the simulated velocity amount and direction in the study area were verified, which confirmed that they were highly consistent with the measured data. Additionally, in line with the crustal movement law, the boundary constraints around the model were removed, and the boundary conditions were changed to the measured velocity data. After importing the data, the east and north velocities were calculated as shown in Figures5 and 6.Figure 5
East velocity isogram in the model after applying tectonic stress (m/a).Figure 6
North velocity isogram in the model after applying tectonic stress (m/a).


## 4.3.3. Tectonic Stress Loading

As the block boundary conditions in different geological periods are unclear, the following approximations were made for block movement in the simulation:(i)
The late Quaternary (~100000–120000 a, considering its average of 110000 a) was considered the loading time of tectonic stress(ii)
During the initial stage of tectonic stress loading, the tectonic stress accumulated in the previous period was fully released(iii)
The calculation time step corresponded to the movement period of the blockAfter considering these approximations, the movement velocity of cells was reset before each operation to ensure uniform movement velocity and boundary conditions inside the block. The distributions of the maximum principal stress and movement velocity of the key stratum after the simulation are shown in Figure7. The distribution of the maximum principal stress direction of the key stratum is shown in Figure 8.Figure 7
Isogram of maximum principal stress and moving velocity of the key stratum after loading tectonic stress (Pa, m/a).Figure 8
Isogram of the southeast angle of the maximum principal stress in the key stratum after loading tectonic stress (°).The maximum principal stress was 21.3–26.4 MPa after applying tectonic stress. The maximum principal stresses in the northwestern, recessed, and southwestern parts of the wellfield were comparatively higher than those near the eastern boundary. Compared to the stress distribution after the initial equilibrium, the area under low maximum principal stress decreased after the tectonic stress was applied, while that under moderate principal stress increased. In general, these changes were more and less prominent in the western and eastern areas, respectively, whereas no significant changes were observed in the central areas. The maximum principal stress direction is S60.07°E-S64.48°E, with a slightly larger offset angle in the southwestern part of the wellfield and a gradually decreasing angle toward the northeastern part.The strata within the wellfield range moved from northwest to southeast in a direction of S72.16°E. The velocity range was 9.40–9.57 mm/a, decreasing from the southwest to northeast of the wellfield, thus suggesting that the wellfield gradually rotated counterclockwise.


## 4.3.4. Calculation of the Stress Concentration Coefficient

To assess the impact risk, Dou et al. proposed the use of a relative stress concentration coefficient composition that considers the ratio of the maximum principal and gravitational stresses caused by a certain factor as the relative stress concentration coefficients since they affect rock burst in a certain target area [20]. The coefficient was used to evaluate the influences of the factor on rock burst.The maximum principal stress at the center point of each cell in the key stratum of the model was exported, and the quotient of the maximum principal stress at the initial equilibrium was calculated to obtain the relative stress concentration coefficient (hereinafter referred to as the stress concentration coefficient). Its distribution isoline is shown in Figure9.Figure 9
Distribution layout of the stress concentration coefficient of the key stratum.The distribution of the stress concentration coefficient differed greatly from that of the maximum principal stress in the key stratum but was similar to the variation law of coal seam thickness to some extent (Figure9). In most areas of the region, the stress concentration coefficient gradually increased from the northeast to the southwest of the wellfield. Generally, the mining height is limited to a certain range as the thickness of the coal seam changes; thus, the greater the thickness of the coal seam is, the higher the mining height. Accordingly, the coal seam in the southwest of the wellfield was found to be a key area for the prevention and control of rock bursts due to its large thickness, high stress concentration, large postmining disturbance range, and high energy release.To compare the differences in tectonic stress in different regions and strata, four measuring points were established in the model: No. 1 was located in the central part of the northwestern recessed area, Nos. 2 and 4 were located in the same monitoring position (different strata) in the center of the wellfield, and No. 3 was located near the eastern boundary of the wellfield. The measuring points for Nos. 1–3 were located in the key stratum, while the measuring point for No. 4 was located in the bedrock of the coal seam roof. The variations in the maximum principal and intermediate principal stresses in the simulation process are shown in Figure10.Figure 10
Stress variation curves of the measuring points.As seen in Figure10, after application of the tectonic stress, the intermediate principal stress, which was oriented in the same direction as the resultant force, increased rapidly. Consequently, it became the maximum principal stress by replacing the vertical stress, while the vertical stress remained unchanged. This indicated that the tectonic stress did not influence the vertical stress but significantly influenced the horizontal stress. Moreover, after the tectonic stress was applied, the increase rate of the horizontal stress in the key stratum was significantly higher than that in the bedrock of the coal seam roof, indicating that the tectonic stress greatly affected the key stratum and that higher energy is stored in the key stratum.Furthermore, the recessed area in the northwestern part of the wellfield exhibited the largest maximum principal stress but smallest stress concentration coefficient. This was primarily due to extensive variation in vertical and horizontal stresses near the recessed area. As a result, the horizontal stress required a longer time to reach and surpass the vertical stress and become the maximum principal stress. Finally, at the end of the simulation, the difference between the two stresses at the recessed area was significantly lower than at other positions, thereby resulting in the smallest stress concentration coefficient.


## 5. Theoretical Analysis of Tectonic Stress Distribution

### 5.1. Determining Tectonic Stress Expression

Suppose that the velocity at any Point A in the wellfield isvA (mm/a) and the velocity at Point K, which is at a distance L from Point A, is vK (mm/a); then, the strain at Point K (εK) in unit time is given as:
(1)εK=limK⟶AvK−vAL.Subsequently, the strain at Point K is given as:(2)σK=E×εK,where σk is the strain of Point K and E is the elastic modulus of the stratum.The relative positions of Points K and A are given as:(3)vK−vA=∂vE∂xi+∂vN∂xjdx+∂vE∂yi+∂vN∂yjdy,where vE is the velocity function of the eastward movement and vN is the velocity function of the northward movement.By inserting Equations (1) and (3) into Equation (2), we obtain:
(4)σK=Elimdx⟶0dy⟶0∂vE∂xi+∂vN∂xjdx+∂vE/∂yi+∂vN/∂yjdyd2x+d2y=E∂vE∂xcosα+∂vE∂ycosβi+∂vN∂xcosα+∂vN∂ycosβj,where α is the acute angle formed by the eastward velocity gradient and x-axis and β is the acute angle formed by the northward velocity gradient and y-axis.


### 5.2. Determining the Parameters Related to Tectonic Stress

#### 5.2.1. Eastward Velocity Analysis

The isolines of the eastward velocity were fitted as shown in Figure11. The eastward velocity distribution function was as follows:
(5)vE=−3.692×10−6x−1.324×10−5y+68.33.Figure 11
Fitting of the eastward velocity isoline (mm/a).Then(6)∂vE∂x=−3.692×10−6,(7)∂vE∂y=−1.324×10−5.


#### 5.2.2. Northward Velocity Analysis

The isolines of the northward velocity were fitted as shown in Figure12. The northward velocity distribution function was:
(8)vN=−1.717×10−10x2+1.202×10−4x−2.298×10−5y+2.317×10−11xy+31.22.Figure 12
Fitting of the northward velocity isoline (mm/a).Then(9)∂vN∂x=−3.434×10−10x+2.317×10−11y+1.202×10−4,(10)∂vN∂y=2.317×10−11x−2.298×10−5.


### 5.3. Expression of the Tectonic Stress

By inserting Equations (6)–(10) into Equation (4), we obtained:
(11)σK=E∙mi+k1x+k2y+nj,where m=−3.692×10−6cosα−1.324×10−5cosβ, k1=−3.434×10−10cosα+2.317×10−11cosβ, k2=2.317×10−11cosα, and n=1.202×10−4cosα−2.298×10−5cosβ.Asα=74.39° and β=18.43°−35.00°, k1<0 and k2>0. m and n are constants; thus, σk decreases with an increase in x and increases with an increase in y. That is, the stress in the northwestern part of the wellfield was the largest, while that in the southeastern part was the smallest. Moreover, the stress changed rapidly along the x-axis but gradually along the y-axis.The elastic modulus of a rock mass is greatly affected by factors such as joints, fractures, and burial depths [21, 22]. Therefore, as it can be influenced by the structure, rock strata dip angle, or local heterogeneity, the horizontal stress in different areas is not uniformly distributed and differs regionally as the elastic modulus changes.


## 5.1. Determining Tectonic Stress Expression

Suppose that the velocity at any Point A in the wellfield isvA (mm/a) and the velocity at Point K, which is at a distance L from Point A, is vK (mm/a); then, the strain at Point K (εK) in unit time is given as:
(1)εK=limK⟶AvK−vAL.Subsequently, the strain at Point K is given as:(2)σK=E×εK,where σk is the strain of Point K and E is the elastic modulus of the stratum.The relative positions of Points K and A are given as:(3)vK−vA=∂vE∂xi+∂vN∂xjdx+∂vE∂yi+∂vN∂yjdy,where vE is the velocity function of the eastward movement and vN is the velocity function of the northward movement.By inserting Equations (1) and (3) into Equation (2), we obtain:
(4)σK=Elimdx⟶0dy⟶0∂vE∂xi+∂vN∂xjdx+∂vE/∂yi+∂vN/∂yjdyd2x+d2y=E∂vE∂xcosα+∂vE∂ycosβi+∂vN∂xcosα+∂vN∂ycosβj,where α is the acute angle formed by the eastward velocity gradient and x-axis and β is the acute angle formed by the northward velocity gradient and y-axis.


## 5.2. Determining the Parameters Related to Tectonic Stress

### 5.2.1. Eastward Velocity Analysis

The isolines of the eastward velocity were fitted as shown in Figure11. The eastward velocity distribution function was as follows:
(5)vE=−3.692×10−6x−1.324×10−5y+68.33.Figure 11
Fitting of the eastward velocity isoline (mm/a).Then(6)∂vE∂x=−3.692×10−6,(7)∂vE∂y=−1.324×10−5.


### 5.2.2. Northward Velocity Analysis

The isolines of the northward velocity were fitted as shown in Figure12. The northward velocity distribution function was:
(8)vN=−1.717×10−10x2+1.202×10−4x−2.298×10−5y+2.317×10−11xy+31.22.Figure 12
Fitting of the northward velocity isoline (mm/a).Then(9)∂vN∂x=−3.434×10−10x+2.317×10−11y+1.202×10−4,(10)∂vN∂y=2.317×10−11x−2.298×10−5.


## 5.2.1. Eastward Velocity Analysis

The isolines of the eastward velocity were fitted as shown in Figure11. The eastward velocity distribution function was as follows:
(5)vE=−3.692×10−6x−1.324×10−5y+68.33.Figure 11
Fitting of the eastward velocity isoline (mm/a).Then(6)∂vE∂x=−3.692×10−6,(7)∂vE∂y=−1.324×10−5.


## 5.2.2. Northward Velocity Analysis

The isolines of the northward velocity were fitted as shown in Figure12. The northward velocity distribution function was:
(8)vN=−1.717×10−10x2+1.202×10−4x−2.298×10−5y+2.317×10−11xy+31.22.Figure 12
Fitting of the northward velocity isoline (mm/a).Then(9)∂vN∂x=−3.434×10−10x+2.317×10−11y+1.202×10−4,(10)∂vN∂y=2.317×10−11x−2.298×10−5.


## 5.3. Expression of the Tectonic Stress

By inserting Equations (6)–(10) into Equation (4), we obtained:
(11)σK=E∙mi+k1x+k2y+nj,where m=−3.692×10−6cosα−1.324×10−5cosβ, k1=−3.434×10−10cosα+2.317×10−11cosβ, k2=2.317×10−11cosα, and n=1.202×10−4cosα−2.298×10−5cosβ.Asα=74.39° and β=18.43°−35.00°, k1<0 and k2>0. m and n are constants; thus, σk decreases with an increase in x and increases with an increase in y. That is, the stress in the northwestern part of the wellfield was the largest, while that in the southeastern part was the smallest. Moreover, the stress changed rapidly along the x-axis but gradually along the y-axis.The elastic modulus of a rock mass is greatly affected by factors such as joints, fractures, and burial depths [21, 22]. Therefore, as it can be influenced by the structure, rock strata dip angle, or local heterogeneity, the horizontal stress in different areas is not uniformly distributed and differs regionally as the elastic modulus changes.


## 6. Analysis of the Measured Data

### 6.1. Analysis of the In Situ Stress Data of the Surrounding Wellfield

Both vertical and horizontal stresses are significantly linearly related to depth [23]. In situ stress data of four mines around the wellfield were collected, and a linear regression analysis was conducted to investigate the relationship between 3D stress and depth. The results are shown in Figure 13.Figure 13
Linear changes in the principal stresses with burial depth along the three directions.Figure13 shows a good fitting effect of vertical stress and burial depth, indicating a high correlation between vertical stress and burial depth. However, the fitting effect of horizontal stress and burial depth was poor, indicating that the horizontal stress was affected not only by the burial depth but also by other factors such as geology, which can explain the high dispersion seen in the data.The regression equations were as follows:(12)σH=0.0249H+3.4206σV=0.0244Hσh=0.007H+5.771.Based on the data point collection, the maximum and minimum principal stresses were both horizontal stresses (σH and σh, respectively), while the intermediate principal stress (σV) was the vertical stress. The relationship between the three stresses was as follows: σH>σV>σh. With increasing burial depth, the maximum principal stress increased the most, followed by the vertical stress and then the minimum principal stress.Based on the regression analysis results of the in situ stress data obtained from the “China Coal Mine Stress Environment Database,” the changes in the maximum and minimum principal and vertical stresses with burial depth were as follows:(13)σH=0.0215H+3.267σV=0.0245Hσh=0.0113H+1.954.Compared with the database statistical results, the regression gradient of vertical stress acquired through the measured data was more accurate, and the regression gradient of the maximum principal stress was slightly larger than the statistical results. The general trend of the regression line was consistent, and the relationship between the stresses of the three directions was the same as that observed for the measured data, that is,σH>σV>σh.


### 6.2. In Situ Stress Test of the Muduchaideng Mine

An in situ stress measuring station was developed while constructing the Muduchaideng Mine. The hole in the measuring point roof was 3.2 m from the roadway floor, and the depth of the top hole was ~25 m. The rock layer (11.9–13.0 m), which is the hydraulic fracturing section, was relatively complete and composed of medium-grained sandstone. The measuring point was at a burial depth of 641 m and an elevation of 649.7 m.The in situ stress test showed that the maximum horizontal principal stress was 13.92 MPa, and the maximum vertical stress was 15.72 MPa. The maximum horizontal principal stress direction of the measuring point was S57.3°E. Using the measured data of the surrounding mines, a multiple regression analysis was conducted to correct the in situ stress test results of the mine [24]. After correction, the maximum principal stress was 20.10 MPa, and the vertical stress was 15.64 MPa.In the cell corresponding to the in situ stress test in the numerical simulation model (measurement point No. 4 in Figure9), the simulation results of the principal stress indicated that the maximum principal stress was 20.12 MPa, the vertical stress was 15.55 MPa, and the direction was S62.7°E. The magnitude is similar to the postcorrection results of the measured data, and the directional deviation rate is less than 10%.


## 6.1. Analysis of the In Situ Stress Data of the Surrounding Wellfield

Both vertical and horizontal stresses are significantly linearly related to depth [23]. In situ stress data of four mines around the wellfield were collected, and a linear regression analysis was conducted to investigate the relationship between 3D stress and depth. The results are shown in Figure 13.Figure 13
Linear changes in the principal stresses with burial depth along the three directions.Figure13 shows a good fitting effect of vertical stress and burial depth, indicating a high correlation between vertical stress and burial depth. However, the fitting effect of horizontal stress and burial depth was poor, indicating that the horizontal stress was affected not only by the burial depth but also by other factors such as geology, which can explain the high dispersion seen in the data.The regression equations were as follows:(12)σH=0.0249H+3.4206σV=0.0244Hσh=0.007H+5.771.Based on the data point collection, the maximum and minimum principal stresses were both horizontal stresses (σH and σh, respectively), while the intermediate principal stress (σV) was the vertical stress. The relationship between the three stresses was as follows: σH>σV>σh. With increasing burial depth, the maximum principal stress increased the most, followed by the vertical stress and then the minimum principal stress.Based on the regression analysis results of the in situ stress data obtained from the “China Coal Mine Stress Environment Database,” the changes in the maximum and minimum principal and vertical stresses with burial depth were as follows:(13)σH=0.0215H+3.267σV=0.0245Hσh=0.0113H+1.954.Compared with the database statistical results, the regression gradient of vertical stress acquired through the measured data was more accurate, and the regression gradient of the maximum principal stress was slightly larger than the statistical results. The general trend of the regression line was consistent, and the relationship between the stresses of the three directions was the same as that observed for the measured data, that is,σH>σV>σh.


## 6.2. In Situ Stress Test of the Muduchaideng Mine

An in situ stress measuring station was developed while constructing the Muduchaideng Mine. The hole in the measuring point roof was 3.2 m from the roadway floor, and the depth of the top hole was ~25 m. The rock layer (11.9–13.0 m), which is the hydraulic fracturing section, was relatively complete and composed of medium-grained sandstone. The measuring point was at a burial depth of 641 m and an elevation of 649.7 m.The in situ stress test showed that the maximum horizontal principal stress was 13.92 MPa, and the maximum vertical stress was 15.72 MPa. The maximum horizontal principal stress direction of the measuring point was S57.3°E. Using the measured data of the surrounding mines, a multiple regression analysis was conducted to correct the in situ stress test results of the mine [24]. After correction, the maximum principal stress was 20.10 MPa, and the vertical stress was 15.64 MPa.In the cell corresponding to the in situ stress test in the numerical simulation model (measurement point No. 4 in Figure9), the simulation results of the principal stress indicated that the maximum principal stress was 20.12 MPa, the vertical stress was 15.55 MPa, and the direction was S62.7°E. The magnitude is similar to the postcorrection results of the measured data, and the directional deviation rate is less than 10%.


## 7. Conclusion

To improve the prevention and control of rock burst, this study examined the in situ stress in the Muduchaideng Mine by numerical simulation using GPS data and investigated the distribution pattern of stress concentration coefficients in key layers. The conclusions are as follows:(1)
Inversion of in situ stress using numerical simulation based on the difference in motion velocity at different locations is feasible under conditions similar to those found in the intact Ordos block(2)
The difference in the movement velocity of the block is the direct cause of the counterclockwise rotation of the Ordos block(3)
Tectonic stress has a large influence on the key layer, and the energy accumulated within the key layer is high; once the key layer is broken, it is easy to trigger impact ground pressure(4)
The stress concentration coefficient in the key layer has some similarity with the change law of the coal seam thickness. The southwestern part of the wellfield is affected by the thickness of the coal seam, which is the key area of impact pressure prevention and control. As such, we should increase our efforts to prevent the negative effects of impact pressure on mine safety and production



---

*Source: 8530609-2023-02-03.xml*