Overlying Strata Dynamic Movement Law and Prediction Method Caused by Longwall Coal-Mining: A Case Study

Abstract

The surface subsidence caused by underground mining is a spatiotemporal process. The impact of mining on surface structures (houses, highways, railways, dikes, etc.) and structures in rock strata (shafts, roadways, chambers, etc.) is a dynamic process. It is necessary to study the dynamic movement law of the surface and overlying strata in the mining process of the working face to predict the extent of the impact of mining on the aforementioned structures. It provides a reference for pre-reinforcement and post-mining treatment. This paper studies the variation of surface dynamic movement based on the survey line above the working face of Peigou Coal Mine. The numerical simulation model of the overlying strata dynamic movement is established to study the dynamic movement law of rock strata with different depths, and the fitting function of surface and overlying strata dynamic movement is determined. Finally, the subsidence velocity prediction function of the major section of the surface and overlying strata in the Peigou Coal Mine is established. The accuracy of this prediction function is demonstrated by contrasting the subsidence and subsidence velocity curves of the surface subsidence basin survey line with the numerical model. In this paper, a numerical simulation method for the dynamic movement of the surface and overlying strata and a function for predicting the subsidence velocity on the strike major section are established, which provides an important theoretical reference for the dynamic protection of the structures on the surface and in the overlying strata.

1. Introduction

The surface subsidence caused by underground mining is a spatiotemporal process. In the mining process, the impact of mining varies for various areas of the surface, and its influence is dominated by the relative position of the working face and the surface area. In production practice, it is nowhere near sufficient to solve the actual problems on site according to the final subsidence law [1,2]. It is often necessary to study the dynamic subsidence and master the variation of subsidence velocity of surface and overlying strata in order to judge the intensity and location of surface and overlying strata movement and deformation, as well as to protect and repair surface structures (houses, dams, roads, and railways, etc.) and structures in overlying strata (chambers, roadways, etc.) [3,4,5].

The dynamic movement of the surface and overlying strata caused by coal mining has been deeply studied in China and other countries. Considering the changes in stress-strain state of rock mass and coal mass around the goaf and the working face in the process of mining [6,7,8,9,10], there is a close relationship between surface movement and overlying strata movement. Shu [11,12] proposed a theoretical model of the relationship between sub-surface and surface subsidence movements, which can be used to predict sub-surface subsidence components at any point within the zone influenced by the extraction of a panel according to the subsidence components obtained from field measurements. Based on the subsidence data in the Southern Coalfield of the Sydney Coal Basin, the proposed prediction is examined. Amar [2] studied the behavior of dynamic active and residual subsidence for a few panels of Jharia coalfield, and found that the subsidence and slope were linearly related to time. Compressive and tensile strains showed a typical fluctuating characteristic behavior, and the rate of mining being a key and controlling parameter for the rate of subsidence and its inter-relationship was developed. The time effect is studied in the coalfield of Asturias, Northern Spain, with the aim of predicting subsidence phenomena and characterizing the trough in the different intermediate stages of the process of excavation and subsidence. In addition, the subsidence is predicted following the models of Knothe and Sroka–Schober and a new time function based on the normal distribution function [13].

According to the difference of surface subsidence velocity in the process of surface dynamic subsidence, Huang [14] divided the whole process of surface subsidence into three stages: Subsidence development, full subsidence, and subsidence attenuation. The development and variation laws of surface subsidence deformation in these three stages were respectively analyzed, and a new concept of “deformation velocity in surface subsidence” was also put forward. Chudek [15] analyzed the impact of overlying strata lithology on the time influence parameters of surface movement and determined the relationship between these parameters and mining depth according to a large amount of survey data. Knothe [16], a foreign scholar, introduced the time function into the prediction of dynamic surface subsidence. Numerous scholars have conducted extensive studies based on Knothe time function, and the findings have been able to accurately predict the entire process of surface subsidence. Scholars have worked to improve the Knothe time function in light of its drawbacks and defects, and they have provided the ideal distribution form of time function to make up for its shortcomings. Based on the static prediction model of probability integral, Zhang [17,18,19] discussed the prediction method of dynamic subsidence for the inclined main section using the “optimized segmented Knothe time function”. The outcome showed that regardless of which direction the mining begins from or whether the prediction direction is uphill or downhill, the results of dynamic subsidence and deformation of the inclined major section are always consistent with the rules of the theory. Guo [20] analyzed the dynamic surface movement characteristics and parameters affected by fully mechanized mining under thick hydraulic collapse, based on the data from the monitoring station. The results indicate that the surface movement is intensive, the surface subsidence velocity is rapid, and the surface deformation and damages are serious in this mining area with thick loess. Deng [21] used the relationship between the maximum subsidence velocity and the relative position of the working face to calculate the prediction formula of the subsidence velocity at any point and time during the mining process. For monitoring the movement of overlying strata, vertical boreholes were drilled from the ground into the rock stratum, Wang [22] adopted the distributed optical fiber sensor monitoring technology and multipoint borehole extensometer methods to the movement data of the strata. Guo [23] conducted field monitoring of overburden displacement, stress, and water pressure changes at the longwall panel in Anhui, China. A three-dimensional annular-shaped overlying zone along the perimeter of the longwall panel is identified for optimal methane drainage during mining. The above researches mainly studied the process and stage of surface subsidence induced by underground longwall coal mining, and obtained the prediction model of surface dynamic subsidence. However, for the protection of surface and underground structures, the distribution characteristics and prediction of subsidence velocity are crucial. In the existing researches, the impact of mining extraction area on surface dynamic movement characteristics is not considered, which makes it impossible to effectively predict the whole process of the working face mining. Moreover, due to the difficulty of overlying the strata field monitoring experiment, the research on overlying strata movement is mainly a static movement, and the relationship between the dynamic movement of the ground surface and the dynamic movement of the overlying strata is not established. As a result, there are few research results on the dynamic movement law and prediction theory. This paper first examines the relationship of movement and deformation between surface and overlying strata, and then analyzes the variation characteristics of subsidence value and subsidence velocity at the maximum subsidence point along the survey line with different goaf areas. The dynamic relationship between the surface dynamic movement parameters (maximum subsidence velocity and lag distance of maximum subsidence velocity) and the advancing distance of the working face is analyzed. The numerical simulation model of the dynamic movement of overlying strata is established, and the rock mechanics and dynamic excavation parameters are determined. When comparing the subsidence and subsidence velocity curves of the surface survey line and monitoring the line of the numerical model, it is considered that the parameters of the numerical model are accurate and can provide appropriate data for analyzing the movement of overlying strata. The variation functions and lag distance of the maximum subsidence velocity at different cover depths of rock strata during the advancement of the working face are determined based on numerical simulation results. The subsidence velocity prediction function on the major section of the surface and overlying strata subsidence basin is obtained based on the distribution function of surface subsidence velocity with supercritical extraction on the major section and the dynamic subsidence prediction method of the major strike section.

In this paper, based on the surface survey data above the panel in Peigou Coal Mine, the variation characteristics of surface dynamic movement parameters, the maximum subsidence velocity (MSV), and the lag distance of maximum subsidence velocity (LDMSV) are obtained during different advancing distances of the working face. To explore the overlying strata movement law, a numerical simulation model is established, and variation functions of dynamic movement parameters for rock strata with different depths are obtained based on the results of the numerical simulation. The prediction model of subsidence velocity at any time and any position on the strike major section of the subsidence basin is established through the dynamic description of surface and overlying strata movement. The study’s findings not only provide guidance for surface and overlying strata structures protection and treatment technologies, but also further enrich the body of knowledge regarding the dynamic movement of the surface and overlying strata.

2. Relationship between Surface and Overlying Strata Movement and Deformation

In the past, the overlying strata movement mainly focused on the strata pressure around the working face, the fluid migration of the overlying strata aquifer, the mining-induced fracture evolution, and relieved methane delivery. These factors form an important basis for addressing a number of safety and environmental issues, such as ground pressure accidents, water and gas disasters, surface subsidence, etc. [24]. Studies on surface subsidence focus on environmental, facility, life, and ecological impacts. However, there is a link between overlying strata and surface subsidence, and mining excavation will inevitably cause the movement of the overlying strata. Considering the key layer as the main body of rock strata movement research, mechanical methods are used to explain changes in the structural form, stress field, and mining-induced fracture field of rock strata. Therefore, there is a relatively unified understanding and complete mechanical description of ground pressure, mining subsidence, and water and gas migration in mining-induced rock strata [25].

The relationship between surface subsidence and rock strata movement was studied by similarity simulation and numerical simulation tests. Adhikary [26] established a practical hydrogeological model based on the large coal seam sub-critical longwall panel in Australia, and stated that the rock collapse and fracture zones are trapezoidal large structures similar to the “inverted funnel” arch. Zuo [27] used the mechanical model to explain the complex mechanical relationship between strata movement and surface subsidence based on the hyperbola-like shape of the broken boundary of strata movement. The boundary of surface subsidence and strata movement was unified as a whole, and an overall strata movement model that can describe the movement of overlying strata to surface subsidence was established. The movement and deformation of overlying strata caused by coal mining are a mechanical phenomenon in rock mass. Due to the presence of rock mass and the complexity of mining boundary conditions, it is difficult to develop a complete mathematical model and its parameters in order to express the overlying strata damage and movement deformation caused by coal mining. According to the laws of strata subsidence and displacement at different buried levels [28,29], as shown in Figure 1, it is assumed that the subsidence curve of different strata considers the projection line of coal wall as the inflection point, and the relationship between overlying strata and surface movement and deformation lies in the following:

(1)

The formation of goaf space will cause the bending deformation of the lower rock strata. After the upper rock strata lose support, their weight and the load of the overlying strata will be transferred to the top of the coal pillar, forming an incremental stress (called abutment pressure) that is greater than the original stress and acting on the lower coal and rock mass to cause compression deformation. The subsidence value of the strata outside of the inflection point is caused by the compression deformation of overlying strata above the coal pillar, which accumulates layer by layer from the coal seam to the ground surface, and shows a gradual increasing trend as the distance between the strata and the coal seam increases.

(2)

Although the fractured rock above the goaf is recomputed by the upper pressure, the rock strata subsidence coefficient from the goaf to the surface gradually decreases due to the plastic expansion of the caved zone and fractured zone above the goaf and the elastic expansion effect of the continuous deformation zone [25]. Due to the weakening of the mining influence, the subsidence coefficient decreases more slowly over time.

(3)

Considering the cantilever effect of the strata at the side of the coal wall, the inflection point of the subsidence curve (the subsidence value is 0.5 W₀) is projected onto the calculation boundary on the coal seam. The distance between the actual mining boundary and the calculation boundary becomes the offset of inflection point. The surface and the strata outside of the inflection point produce tensile deformation in the horizontal direction, and the position of the maximum displacement is located at the inflection point. On the other hand, compressive deformation of the surface and upper strata occurs in the horizontal direction within the inflection points on both sides. However, for the lower strata, they are compressed both horizontally and vertically outside of the inflection point. Within the inflection point, they are subject to tensile deformation in the horizontal direction due to the accumulative effect of broken rock mass. From the perspective of rock strata mining influence range, a subsidence basin that is significantly larger than the goaf area appears on the surface with the expansion of the damage range of the rock strata.

In summary, due to the complex mechanical behavior of the mining-induced rock mass, rock strata movement is closely related to the surface movement. The surface subsidence basin is the external manifestation of the rock strata movement, and the rock strata movement is the fundamental reason for the surface subsidence. Due to the complexity of the occurrence of the rock mass itself and the mechanical boundary condition, it is impossible to fix many in-situ monitoring points for the movement of boreholes in the overlying strata as the surface survey stations [30,31]. Therefore, based on the measured surface data, this paper establishes a numerical simulation model of overlying strata dynamic movement and studies its law. Combined with the aforementioned information, the dynamic movement prediction model of the surface and overlying strata is established.

3. Surface Dynamic Movement Characteristics

3.1. Geological Setting and Survey Lines Plan

The Peigou Coal Mine is located in Zhengzhou mining area of the Henan Province, China (Figure 2a,b). The panel 31,071 has a width of 130 m, a strike length of 1100 m (Figure 2c), an average mining height of 7.5 m, an average incline angle of 15°, and an average burial depth of 300 m. The mining method is a fully mechanized top-coal caving. The roofs of the mining coal seam are mainly fine sandstone and sandy mudstone. The sandstone fissure aquifer of lower Shanxi formation is the upper main aquifer in the local coal mine, and this aquifer is widely distributed in the entire coal mine with an average thickness of 37.13 m. The distance between the bottom interface of this aquifer and the coal seam is approximately 6 m. The aquifer consists of sandstones with medium and coarse grain sizes. Moreover, fractures have been found in certain local regions of the aquifer. In general, the aquifer in Shanxi formation exhibits poor permeability and a weak water yield property. However, there is a reservoir called Modongwang Up Neath (panel 31071), where the risk of reservoir water inrush is high. In order to accurately evaluate the impact of mining damage on Modongwang reservoir and its dam, a surface movement survey station is established on the surface above the panel. Before mining, movement survey lines in the strike and the inclined directions are arranged on the surface above the panel (Figure 2c).

3.2. Subsidence and Subsidence Velocity Curves of Survey Station with Maximum Subsidence

Figure 3 shows the variation curves of the subsidence velocity and subsidence value at the survey station QZ in the surface strike survey line. Due to survey error, surface subsidence also occurs when the working face is distant. Therefore, this paper mainly considers the increase in subsidence velocity as the basis for judging the beginning of the mining influence. When the distance between the working face and the survey station is about 54 m, this station begins to subside. It can be determined that the working face’s leading distance is 54 m, and the advance effect angle is equal to 79.27°. When the working face passes over this station, the subsidence velocity increases gradually, and then at a distance of 76 m, the subsidence velocity reaches its maximum value of 38.36 mm/d. The surface position subsides most violently at this point. Then, the subsidence velocity gradually decreases, and the subsidence curve of the maximum subsidence position essentially flattens out. When the working face passes over this station by roughly 296 m, the subsidence velocity of this station is 1.67 mm/d, the active phase of the station subsidence ends, and the recession phase starts. Then, the subsidence velocity of this station slows down, and its subsidence value becomes stable at 2140 mm.

3.3. Variation of Surface Dynamic Movement Parameters in the Working Face Mining Process

• The MSV of surface subsidence basin

To express the influence induced by underground mining to the surface with different mining area sizes, subsidence velocity is considered as an index to measure the severity of mining influence. The relationship between the maximum subsidence velocity on the major section of the surface subsidence basin and the advancing position of the working face is depicted in Figure 4.

Based on the Knothe time function of the surface dynamic subsidence, the MSV(V_max) function formula is obtained by curve fitting:

$$V_{\max }=50.40(1-e^{-0.0028d})$$

(1)

In Equation (1), d is the mining distance of the working face, m.

It can be seen from Figure 4 that with the advancement in the working face, the goaf area increases, and the influence of underground mining on the surface intensifies. As a result, the MSV on the major section of the surface trend gradually increases. When the mining distance of the working face reaches 663.8 m, the V_max (MSV) of surface survey station S9 increases from 0 to 45.43 mm/d. However, when the mining distance exceeds 400 m (i.e., 1.4 times the mining depth), it reaches supercritical mining in the strike direction, and the increase amplitude of the MSV station gradually decreases, reaching a stable value of 50.40 mm/d.

2.

Variation of LDMSV

On the surface subsidence velocity curve, the position of the MSV always lags behind the working face by a certain distance. This phenomenon is called the lag phenomenon of MSV. Knowing the LDMSV of the surface is useful in identifying the area where the surface moves violently during the working face mining process and the time when the maximum subsidence velocity occurs, which is significant as a guide for the protection of surface buildings. It is well known that when the working face starts mining from the open-off cut, each station on the surface undergoes a process from subcritical mining to supercritical mining, and LDMSV will be a dynamic process. Therefore, studying the variation of LDMSV can dynamically determine the area with the most violent surface movement in the mining process.

To describe the variation of LDMSV on the strike major section of the surface subsidence basin during the mining process, the nonlinear functional relationship between the LDMSV(L) and the mining distance is obtained by curve fitting:

$$L=95.60(1-e^{-0.11d})$$

(2)

Figure 5 illustrates the relationship between the LDMSV and the mining distance as an exponential function. The LDMSV significantly increases with the mining distance before the latter is advanced to 400 m. When the working face reaches supercritical mining, the LDMSV gradually flattens, increases to a certain extent, and then essentially does not increase, maintaining stability at 95.50 m; namely, the lag angle of MSV is 71.47°. It shows that in the mining process in this coal mine, the proportion of mudstone and sandstone in the overlying strata is large, and the overlying strata structure is weak. As the goaf area increases, the strata movement is quickly transmitted to the surface, resulting in a large lag angle for the MSV on the surface.

4. Numerical Simulation of Dynamic Movement of Overlying Strata

4.1. Establishment of Rock Strata Movement Numerical Model

The strike section of the panel is selected as the calculation object in the numerical model. The dimensions of the model are 350 × 1000 m, the mining distance is 440 m, and the coal pillars of the numerical model are 275 m and 285 m, which can ensure that the left and right boundaries of the model are outside of the range of the influence angle of strata movement and avoid the influence of boundary conditions on the movement caused by the excavation. As for the boundary conditions of the model, the left and right sides of the model have roller constraints, the bottom has full constraints, and the top has free boundary conditions. The input parameters of the numerical model rock mass are repeatedly checked by adjusting the geological strength index (GSI) [32,33,34], in order that the numerical model can truly reflect the strata movement. These parameters are obtained by the inversion method. The physical and mechanical parameters of coal and rock mass and the parameters related to Hoek-Brown strength criterion are shown in

Table 1. Mechanical property parameter of rock and coal in numerical model.

Position	lithology	Density/(kg·m−3)	Elastic Modulus/GPa	Cohesion/MPa	Internal Friction Angle/°	Uniaxial Compressive Strength/MPa	Tensile Strength/MPa	GSI	mi
Roof	Mudstone	2300	2.97	1.13	21.56	35	0.07	38	5
	Sandy Mudstone	2300	3.98	1.78	23.49	50	0.09	40	6
	Sandstone Interbedding	2400	4.09	1.63	25.26	42	0.07	42	7
	Medium Grained Sandstone	2550	6.91	4.14	30.53	85	0.11	45	12
	Fine Grained sandstone	2550	7.50	4.87	29.04	106	0.16	45	10
Coal Seam	Coal	1400	1.58	0.88	24.60	25	0.01	30	10
Floor	Mudstone	2300	2.97	1.13	21.56	35	0.07	38	5
	Medium Grained Sandstone	2300	6.91	4.14	30.53	85	0.11	45	12
	Limestone	2500	5.71	3.82	27.00	92	0.12	41	9

.

4.2. Verification of the Numerical Model of Strata Movement

The rationality of the numerical model directly determines the accuracy of the numerical simulation analysis. In this paper, the accuracy of the numerical model is assessed through the comparative analysis of the survey data of the surface movement survey line, the theoretical calculation value of the abutment pressure, and the theoretical calculation value of the stress-recovery distance in the goaf with the numerical simulation.

(1) The uniaxial compressive strength of intact rock sample σci is 25 MPa, the GSI is 30, and the coal mechanical parameters are obtained by Willson’s theoretical formula [35].

As shown in Figure 6, the peak abutment pressure of the coal seam is 17.50 MPa, the stress concentration factor K is 2.49, the distance between the peak abutment pressure and the coal wall is 16.21 m, and the influence distance before the peak abutment pressure is 65.86 m. According to the numerical simulation, the abutment pressure peak value is located 15 m near the front of the coal wall, and its peak value is 18.52 MPa. The errors of peak stress value and peak position are 5.51% and 8.06%, respectively. Moreover, according to the fact that the ground pressure of the headentry and tailentry of the panel is apparent in the area of 14~17 m, the correctness of the numerical model for the calculation value of abutment pressure is verified.

(2) According to the numerical model stress monitoring in the goaf, it is found that the stress in the goaf at 101.5 m behind the coal wall reaches the original rock stress, and the corresponding stress-recovery distance in the goaf is 101.5 m (i.e., 0.6 times the mining depth), conforming to the laws of domestic and foreign experience, which is the mining depth of 0.3~0.4 times the stress-recovery distance in the goaf.

(3) The parameters of the surface probability integral method are obtained based on research on the strike survey line over the panel: The subsidence coefficient q = 0.80, the tangent of influence angle tan β = 1.88, and the offset of inflection point is 0.15 H, where H is the average depth of coal seam.

The subsidence data of the surface are extracted with the mining distance of 440 m, and the numerical model subsidence curve is obtained and used for comparison with the data calculated by the probability integral method. The subsidence curves of the measured data and numerical model are shown in Figure 7. It can be surveyed that the subsidence coefficient q is 0.80 with supercritical mining. The maximum subsidence obtained by the numerical model is 5.78 m, i.e., the subsidence coefficient is 0.77, and thus the error rate is 3.66%. It can be considered that the measured results and the numerical model results are consistent in terms of accuracy and distribution law. This shows that the numerical model can accurately inverse the rock strata movement and can be used as an effective means to study the law of overlying strata and surface movement [36].

4.3. Determination of Dynamic Excavation Parameters for the Numerical Model

In the discrete element numerical simulation software, the calculation of the model is mainly controlled by two different methods [37]. The first method is by setting the number of the software’s operation cycle steps, and stopping the operation before moving on to the next stage. Moreover, the operation will stop or move on to the next stage if the maximum unbalance force ratio (the ratio of the maximum unbalanced force to the initial unbalanced force) is less than the default value even if the number of calculation steps has not reached the set cycle steps. The second method is by setting the maximum unbalanced force ratio (MUFR). When the maximum unbalanced force ratio calculated by the software is less than the aforementioned set value, the operation will stop or move on to the next stage. The calculation principle is shown in Figure 8.

For the dynamic representation of coal seam excavation, a certain excavation step is generally set, and a certain cycle step is given, or the value of MUFR is changed. When the calculation is completed, the model enters the next excavation operation. However, in the continuous excavation, the deformation and failure of the overlying strata are a constantly changing process. It is necessary to set appropriate excavation steps and cycle steps (or MUFR) to accurately simulate an actual mining operation. Therefore, it is necessary to check the parameters used in the numerical simulation through relevant data measured on site, in order that the numerical simulation can be used as a reliable means to study the movement of the surface and overlying strata.

(1)

Determination of excavation step

Kbj-60 Ⅲ rock pressure observation system is installed in Peigou Coal Mine to accurately grasp the ground stress and other data of roof strata. During the mining, the monitoring system continuously logs the stress of the hydraulic support for the underground working face in the support extension. The extension data are gathered by specialized personnel and uploaded to a specialized computer software. Based on the statistical data of the mine pressure monitoring system during the mining, it is known that the average periodic stress step distance of the panel is 18.5 m, and the numerical simulation excavation step distance is determined to be 20 m.

(2)

Setting of MUFR

According to the schematic diagram of the dynamic calculation scheme in Figure 9, an initial value is first set, and the numerical simulation calculation is compared with the measured surface subsidence curve. Then, the maximum unbalanced force ratio is continuously adjusted in order that the numerical simulation results are consistent with the field-measured results. It is found that when the excavation step is 20 m, the maximum unbalanced force ratio is 4.8 e⁻⁵. After the numerical simulation of dynamic excavation, the dynamic subsidence data of the surface during the mining process are obtained, and the subsidence curve of survey station Z8 and the corresponding numerical model positions are drawn, as shown in Figure 10. The subsidence curves of the survey stations on the surface and the corresponding points in the numerical model are very similar. The relationship curve between the survey station subsidence and mining distance presents an asymmetric “s” shape with a small front end and a large rear end. The two curves can show three stages of dynamic development of the surface.

Figure 11 illustrates a comparison of measured data and numerical simulation results of subsidence velocity. It can be noticed that when two adjacent survey stations with large measurement error are excluded with a mining distance of about 400 m, the curve of subsidence velocity obtained by measurement is similar to the one obtained from the numerical model in curve shape, the MSV on the curve are 124.32 mm/d and 132.78 mm/d, and the relative error rate is 6.4%. This demonstrates that the numerical model can accurately calculate the surface dynamic movement and can therefore be used to study the dynamic movement of the surface and overlying strata.

4.4. Dynamic Movement Characteristics of Overlying Strata

(1)

Variation of overlying strata MSV

Considering the surface subsidence curve of the numerical model as an example, the MSV of the surface is fitted according to the analysis method of the surface dynamic movement, as shown in Figure 4. The functional relationship between the surface MSV and the mining distance is obtained:

$$V_{\max }=137.96(1-e^{-0.0083d})$$

(3)

The subsidence data of rock stratus are then selected to study the dynamic movement, which are 250 m, 210 m, 170 m, 130 m, 90 m, 50 m, and 10 m away from the coal seam. Due to its discreteness, the MSV of rock strata under different mining distances is fitted and defined by referral to the form of the relationship function between the MSV of the surface and the mining distance (Equation (3)):

$$V_{\max }^i=V^{\prime }(1-e^{-bd})$$

(4)

In Equation (4), Vⁱ_max is the MSV of the strata with the i depth, mm/d; V ’is the stable value of the MSV, mm/d; and b is the growth coefficient.

By substituting the MSV of different depths of strata in the mining process into Equation (4), the parameter of the dynamic change function of the MSV of strata is obtained as shown in

Table 2. The MSV parameter of the overlying strata.

Distance to Coal Seam (m)	300	250	210	170	130	90	50	30
V′ (mm/d)	137.96	139.04	141.48	162.45	165.5	188.41	231.3	261.58
b	0.0083	0.0085	0.0093	0.0083	0.0093	0.0095	0.0127	0.01

.

The MSV of different depths of strata gradually increases with the mining distance. When the working face advances to a certain distance, the MSV of strata tends to have a fixed value. It can be seen from Table 2 that the MSV of the rock strata gradually increases with the distance, which is between the rock strata and the coal seam. The MSV of the surface increases from 137.96 to 261.58 mm/d when the distance from the coal seam is 30 m. To express the relationship between the MSV after the rock strata is stabilized and the distance between the rock strata and the coal seam, the MSV of the rock strata in Table 2 is nonlinearly fitted. In addition, the fitting formula (Equation (5)) is obtained, as well as the variation curve of the MSV at various distances between the rock stratus and the coal seam (Figure 5):

$${V^{\prime }}_{\max }(z)=V_0{(\frac{H_0-\mathrm{z}}{H_0})}^{-0.30}$$

(5)

In Equation (5), z is the depth of the rock strata, m; H₀ is the depth of the coal seam, m; and V₀ is the MSV of the surface, mm/d, i.e., 132.75 mm/d.

It can be seen from Table 2 that the growth parameter B in the fitting parameter of the MSV of the rock strata increases with the decrease in the distance from the coal seam, indicating that when the MSV of the rock strata close to the coal seam is stable, the mining distance is small. According to the curve fitting of the above data, the relationship between the growth parameter b and the depth of rock strata z is obtained:

$$b(z)=0.008{(\frac{H_0-\mathrm{z}}{H_0})}^{-0.12}$$

(6)

Substituting Equations (5) and (6) into Equation (4), the relationship between the MSV of rock strata and the mining distance is obtained:

$$V_{}^i=V_0{(\frac{H_0-z}{H_0})}^{-0.30}[1-e^{-0.008{(\frac{H_0-z}{H_0})}^{-0.12}d}]$$

(7)

(2)

Variation of overlying strata LDMSV

Due to the different distances between the rock strata and the coal seam, the mining-induced impact and the dynamic movement of the strata are also different. This paper selects the numerical calculation data with strata depths of 250 m, 210 m, 170 m, 130 m, 90 m, 50 m, and 10 m in order to study the variation of the LDMSV of rock strata in the mining process.

Similar to the form of Equation (4), L’ is defined as the LDMSV of rock strata, and the parameter c represents the growth trend of the LDMSV with the mining distance, as shown in

Table 3. The LDMSV parameter of the overlying strata.

Distance to Coal Seam (m)	300	250	210	170	130	90	50	30
L′ (mm/d)	79.51	76.53	59.92	55.13	55.13	49.89	45.68	40.36
c	0.014	0.015	0.014	0.013	0.013	0.014	0.011	0.013

.

When the distance between the rock strata and the coal seam gradually decreases, L’ gradually decreases; namely, from 79.51 to 40.36 m when the distance between the rock strata and the coal seam is 30 m. However, L’ of the rock strata has little change with the growth parameter c, which is between 0.011 and 0.015, and its average value can be considered as 0.013. Similarly, the stable value L’_max of the lag distance from the MSV of the rock strata can be determined as shown in Equation (8). Therefore, the calculation formula of L’ with different mining depths can be obtained (Equation (9)):

$${L^{\prime }}_{\max }(z)=41.84\frac{H_0-\mathrm{z}}{H_0}+36.32$$

(8)

$$L^i=(41.84\frac{H_0-\mathrm{z}}{H_0}+36.32)(1-e^{-0.013d})$$

(9)

5. Calculation Method of Subsidence Velocity

5.1. Calculation Method of Subsidence Velocity in Strike Major Section of Surface Subsidence Basin

According to the variation of the subsidence velocity of the survey stations, the study in [13] shows that the distribution curve of the surface subsidence velocity is similar to the “quadratic curve” distribution. If the projection position on the surface at a certain time is considered as the coordinate origin, the direction of the working face on the strike is the positive direction of the x axis, and the subsidence velocity of the station is the y axis. Therefore, the relative relationship between the subsidence velocity curve of the ground surface and the working face position can be expressed. The subsidence velocity formula at one position on the strike major section is:

$$V(x)=\frac{V_{\max }}{1+{(\frac{x+L}{a})}^2}$$

(10)

In Equation (10), x is the distance to the working face, a is the morphological parameter, which indicates the steepness of the curve.

For the solution of the shape parameter a, this paper assumes that the a calculation point on the surface starts when it is significantly distant from the working face, and its subsidence tends to have a stable value when the working face pushes over significantly away from the point. The following results are obtained:

$$\Delta W={\displaystyle {\int }_{+\infty }^{-\infty }\frac{V_{\max }}{1+{(\frac{x+L}{a})}^2}}d(-\frac{x}{c})$$

(11)

Then:

$$a=\frac{W_{\max }}{V_{\max }}\cdot \frac{c}{\pi }$$

(12)

In Equation (12), W_max is the maximum surface subsidence and can be obtained by $W_{\max }=qm\sqrt[3]{n_1{n}_2}\cos \alpha$; m is the thinness of coal seam, n₁ and n₂ are respectively equal to the ratio of mining size and coal seam burial depth on the strike and incline direction, and c is the average mining velocity of the working face, m/d. Moreover, V_max can be calculated according to Equation (1).

From Equation (12), it can be seen that when the mining velocity of the working face is constant, the value of the shape parameter a changes with the change in the maximum subsidence and maximum subsidence velocity values. Moreover, it can be judged that before reaching the supercritical mining, the shape parameter a and parameter L of the subsidence velocity curve will change with the increase in the goaf area, and the change law can be obtained from the relationship between the MSV, LDMSV, and mining distance. After a mining distance of 400 m, the maximum subsidence W_max no longer increases due to supercritical mining in the strike direction of the working face, while the MSV of Equation (1) increases, but the increasing range gradually decreases. Therefore, it can be considered that the parameter a of the subsidence velocity curve increases slightly after supercritical mining, and the shape of subsidence velocity curve becomes steeper.

In summary, the variation of the subsidence velocity curve with the mining process can be described as follows: With the increase in the goaf area, the MSV gradually increases, the subsidence velocity curve shape gradually becomes steep, and the LDMSV gradually increases. When the working face reaches supercritical mining; namely, the working face mining distance exceeds 400 m, the subsidence velocity curve of the strike major section maintains a certain lag distance with the working face in a fixed form and moves forward with the mining. Figure 12 shows the relative position of the subsidence and subsidence velocity curves of the strike major section of the surface, when the working face mining distance is 200m, 400m, 600m, and 800 m, respectively.

The mining and geology parameters of the panel of Peigou Coal Mine are brought into Equations (1), (2), and (12), and the formula for calculating the subsidence velocity at any position on the strike major section at any time of the working face is obtained by bringing in Equation (10):

$$V(x,d)=\frac{50.40(1-e^{-0.0028d})}{1+(\frac{x+95.60(1-e^{-0.11d})}{qm\sqrt[3]{n_1{n}_2}\cos \alpha \cdot c}\cdot 50.40(1-e^{-0.0028d})\cdot \pi {)}^2}$$

(13)

When the working face mining distance is equal to 360 m, the subsidence velocity calculation points on the strike major section of subsidence basin at different mining distances are obtained according to Equation (13) and compared with the measured values, as shown in Figure 13.

It can be seen from Figure 13 that when the mining distance of the working face is 360 m, the deviation between the measured and calculated values of the subsidence velocity at each point is small, and the average error is 1.57 mm/d, indicating that the predicted results of dynamic subsidence velocity can meet engineering needs.

5.2. Prediction Model of Overlying Strata Dynamic Movement

For the correction of the LDMSV of rock strata, a correction factor a_L can be added before Equation (9), and its value is equal to Ls/L´(z)|z = H₀. Therefore, the calculation formula of LDMSV of overlying strata based on the measured LDMSV of the surface is obtained:

$$L^i=a_L(41.84\frac{H_0-z}{H_0}+36.32)(1-e^{-0.013d})$$

(14)

Assuming that the subsidence velocity distribution curve on the major section of the surface and overlying strata is similar to the form of a quadratic curve during the working face mining, referring to the subsidence velocity function of the strike major section of the surface (Equation (10)), the subsidence velocity formula of the overlying strata is:

$$V^i(x)=\frac{V^i{}_{\max }}{1+{(\frac{x+L^i}{a^i})}^2}$$

(15)

According to Equation (11), the parameter a of the subsidence velocity curve will change with the increase in mining distance. When the working face moves from subcritical mining to supercritical mining, the shape of the subsidence velocity curve becomes steeper. The shape parameter aⁱ of the distribution function of rock strata subsidence velocity is mainly calculated according to the maximum subsidence value Wⁱ_MAX, the Vⁱ_max (MSV), and the working face mining velocity c of the rock strike under this level. The maximum subsidence velocity Vⁱ_max and the rock strata can be referred to Equation (5), and the maximum subsidence Wⁱ_MAX of the rock strata can be determined according to the mining distance and the distance between the rock strata and the coal seam; namely:

$$W_{\max }^i=M{q}_i\sqrt[\beta ]{n_1^i{n}_2^i}\cos \alpha$$

(16)

In Equation (16), for rock strata i, the mining coefficients nⁱ₁ and nⁱ₂ of the working face are equal to ηD/(H₀-z), where η is a coefficient generally considered as 0.7~0.9, D is the penal length of incline or strike (m), and z is the depth of rock strata i (m).

Therefore, Equations (7), (14), and (16) are combined with Equation (12) and brought into Equation (15) to obtain the calculation formula of subsidence velocity of the surface and overlying strata on the strike major section at different mining distances:

$$V(x,d,z)=\frac{V_{\mathrm{S}}{(\frac{H_0-z}{H_0})}^{-0.30}[1-e^{-0.008{(\frac{H_0-z}{H_0})}^{-0.12}d}]}{1+{\{\frac{x+a_L(41.84\frac{H_0-z}{H_0}+36.32)(1-e^{-0.013d})}{M{q}_i\sqrt[3]{n_1^i{n}_2^i}\cos \alpha \cdot c}\cdot V_{\mathrm{S}}{(\frac{H_0-z}{H_0})}^{-0.30}[1-e^{-0.008{(\frac{H_0-z}{H_0})}^{-0.12}d}]\cdot \pi \}}^2}$$

(17)

According to the above analysis of the MSV and the LDMSV on the strike major section of the surface in Peigou Coal Mine, it is known that in Equation (17), V₀ is equal to 50.40 mm/d, and the correction coefficient a_L of the LDMSV is equal to 1.22. Bringing the above parameters into Equation (17) results in the calculation formula of the subsidence velocity of the surface and overlying strata. When the mining velocity is 2.3 m/d and the mining time is selected as 20d, 40 d, 60 d, 100 d, and 140 d, the subsidence velocity distribution curves of the rock strata with the depths of 100 and 200 m, respectively, are shown in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13.

It can be seen from Figure 14 that the subsidence velocity of rock strata shows the following laws: First, the peak value on the subsidence velocity curve gradually increases with the increase in mining distance, and when the working face mining reaches supercritical mining, the increase amplitude of the MSV gradually decreases, finally stabilizing at a certain value. In addition, the shape of the subsidence velocity curve of the rock strata becomes steeper with the increase in mining distance, and finally tends to have a stable shape. Finally, the LDMSV of the subsidence velocity curve increases gradually with the increase in mining distance. In general, the peak value of the rock subsidence velocity curve gradually increases, and the curve shape gradually becomes steep during the mining. When the mining distance reaches a certain value, the peak value and shape of the rock strata subsidence velocity curve tend to be stable and maintain a certain lag distance from the working face, moving forward with the mining of the working face. By comparing Figure 14a,b, it can be found that with the same mining distance, the mining-induced influence of the strata close to the coal seam is greater, resulting in a larger MSV of the rock strata, a smaller LDMSV, and a steeper shape of the subsidence velocity curve.

The dynamic impact of mining on structures (shafts, roadways, and other structures) in overburden can be predicted and evaluated according to the variation rule of overlying strata dynamic subsidence and movement. During the mining process, the position and time of the most severe rock strata movement can be determined to reinforce, prevent in advance, and ensure the safe operation of underground mining.

6. Conclusions and Discussion

In this paper, the dynamic movement law and prediction method of the surface and overburden are studied with the stratum structure in Peigou Coal Mine. Since the dynamic movement in the overlying strata is not easy to be observed, this paper innovatively put forward a numerical simulation method to study the dynamic movement, which is used to establish the relationship between surface and overlying strata dynamic movement parameters. Therefore, the dynamic movement law and prediction method of overlying strata and surface are obtained. The research results can be used to predict the dynamic subsidence impact on the structures on the surface and in the overlying strata during the mining process, in order to obtain the mining time or position of working face, which induces the strongest impact on the these structures, and provides a reference for the protection of the structures. The key conclusions from this study are as follows:

(1)

In order to study the movement of overlying strata, a numerical simulation model was constructed, and the accuracy of the numerical model was demonstrated by verifying the surface subsidence data, the law of ground pressure, and the recovery distance of the goaf. Moreover, the dynamic excavation parameters of the numerical simulation were determined, and a simulation method of dynamic movement of overlying strata was provided.

(2)

By analyzing the simulation results of the numerical model, the dynamic movement characteristics of the rock mass with different depths during the mining of the working face were obtained, including the variation laws of the MSV and LDMSV. The dynamic movement parameter functions of rock strata with different depths were determined by curve fitting.

(3)

According to the variation law of dynamic movement parameters of surface and rock mass during mining, the prediction formula of subsidence velocity at different positions and mining distances on the strike major section was established. Considering the mining distance of the working face as 360 m and comparing the predicted and measured values of the surface subsidence velocity, the deviation between the measured and calculated values of the subsidence velocity at each point on the surface strike section was small, and the average error was 1.57 mm/d, indicating that the predicted results of dynamic movement can meet the needs of the project.

(4)

At the same depth, the subsidence velocity of rock strata is distributed in the form of a quadratic curve. During the mining process of the working face, the MSV of rock strata subsidence velocity curve gradually increases, and the shape gradually becomes steep. When supercritical mining is reached, the MSV and shape of the rock strata subsidence velocity curve tend to be stable and maintain a certain lag distance from the working face, moving forward with the mining of the working face. Under the condition of the same mining distance of the working face, the MSV of the subsidence velocity curve of the rock strata close to the coal seam is larger, the LDMSV is smaller, and the shape of the subsidence curve is relatively steeper. The position and time of the most severe rock strata movement can be determined to reinforce and prevent in advance and ensure the safe operation of underground mining.

Due to the complexity and diversity of stratum structure, the strata with different structures show different movement laws induced by underground mining, but entire or partial overburden strata movement after mining is controlled by strata structures of “key stratum”. Due to the different locations of the key strata in the overlying strata, the propagation law of the movement and deformation in the rock stratum is different. Another important factor influencing the dynamic movement of overlying strata is the bulking factor of the roof strata in the goaf [38,39,40,41,42,43]. The stratum immediately above the seam bends downward and later breaks into different size blocks depending on their lithological nature. This phenomenon continues upward until the lowest uncaved stratum receives the support of the caved and bulked rock pile. During this process, the strata are broken into blocks by vertical and/or sub-vertical fractures and horizontal cracks due to the vertical bed separation. Vertical and/or sub-vertical fractures andhorizontal crackswill be re-compacted by the stress induced by the weight of overlying strata, and the surface and rock stratum dynamic movement is closely related to the re-compaction of goaf. The above two aspects are the focus of further research.