*Article* **Prediction Method for Surface Subsidence of Coal Seam Mining in Loess Donga Based on the Probability Integration Model**

**Bingchao Zhao 1,2, Yaxin Guo 1,2,\*, Xuwei Mao 1,2, Di Zhai 1,2, Defu Zhu 2,3, Yuming Huo 4, Zedong Sun 3,5 and Jingbin Wang 1,2**


**Abstract:** The accurate prediction of surface subsidence is a significant foundation for the damage assessment of coal seam mining and ecological environment reclamation in loess donga. However, conventional models are very problematic, and the reliability of prediction is usually low. Therefore, we propose a method for predicting surface subsidence of coal seam mining in loess donga that is based on the probability integration model, combined with the movement principle of rock and soil layers in the respective study area, and considering the influence of slope stability and additional mining slip on mining subsidence. The feasibility of our new method was verified by a case study in the N1114 working face of the Ningtiaota coal mine (China) that is situated in an area with abundant loess dongas. The results show that slope slippage is the source of error in the prediction of subsidence in loess donga. The prediction idea of "dividing the surface of loess donga into horizontal strata area and slope sub-area, and predicting the subsidence value of the two areas, respectively" is put forward. A method for predicting the subsidence value of two regions is established. First, based on the theory of probability integral and rock formation movement, the probability integral parameters of the horizontal stratum area are determined, and the subsidence basins in the area are superimposed and calculated. Secondly, according to the slope stability and slip principle, the additional displacement of subsidence in the slope area with mining instability coefficient *G*cs > 0.87 is calculated. Finally, combined with the subsidence prediction results of the strata area and the slope sub-area, and the position of the slope, the accurate prediction of the surface subsidence in loess donga is realized. Our results show that the agreement between the curves predicted from our calculations and from the measured data are between 88.7–97.8%. The calculated error of the additional displacement of slope mining slip is between 1.0–9.8%. The excellent correlation between the modelled and measured data documents that our method provides, demonstrated a new efficient and valuable tool for the precise prediction of damages induced by mining of underground coal seams in loess donga.

**Keywords:** surface subsidence; probability integration; loess donga; superimposed calculation; additional displacement of slope mining slip

#### **1. Introduction**

The Yushenfu coalfield [1] is the largest coalfield in China and one of the seven largest coalfields in the world. Loess dongas are widely distributed in this area. Major coal mining

**Citation:** Zhao, B.; Guo, Y.; Mao, X.; Zhai, D.; Zhu, D.; Huo, Y.; Sun, Z.; Wang, J. Prediction Method for Surface Subsidence of Coal Seam Mining in Loess Donga Based on the Probability Integration Model. *Energies* **2022**, *15*, 2282. https:// doi.org/10.3390/en15062282

Academic Editors: Nikolaos Koukouzas and Adam Smoli ´nski

Received: 15 February 2022 Accepted: 19 March 2022 Published: 21 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

activity in this area has caused an increasing amount of subsidence areas. The subsidence has generated major damages to the cultivated land, water bodies, and buildings to varying degrees within the area affected by the mining, and also caused serious environmental problems [2], such as land desertification and soil erosion. An accurate prediction of surface subsidence is a prerequisite to reduce or even avoid such environmental problems. At present, among many surface subsidence prediction models, the probability integration model [3] is one of the most widely used methods for surface subsidence prediction in mining areas.

Although the model is sufficient several limitations exist: on the one hand, the parameter selection is the key to the probability integration method. However, the specific geological structure of loess donga may lead to large errors [4] in the predicted parameters, influencing the effective constraints of predicted parameters. On the other hand, the stability of the slope and the additional displacement [5] will deviate from the predicted results. We propose a prediction method of surface coal seam mining subsidence in loess donga by considering the essential causes of surface subsidence, i.e., the movement of underground rock layers and the stability of the slope [6], as important factors affecting the prediction accuracy, and combining the probability integral method model, the movement principles of rock and soil layers, and the influence of slope stability on mining subsidence.

The probability integration model [7] is a mining subsidence prediction method based on the random medium theory, with the predicted parameters as the key. Litwiniszyn [8] proposed the model in the 1950s and it has been further developed into a mature application model since then. The effective predicted parameters of the horizontal surface can be obtained in various ways. Therefore, the accuracy and reliability of the predicted results of the probability integral model are assured. However, concerning the prediction of the surface subsidence in donga areas, the model is affected by the probability integral parameters and the slope sliding effect [9], which leads to a decline of the application effect. Several studies have been conducted to solve this problem.

Field measurements [10] are a common method to obtain detailed and reliable predicted surface parameters. However, an accurate measurement is very time-consuming and requires a long period (at least 1 or 2 years), which squanders a lot of manpower and material resources. Moreover, the obtained parameters are only applicable to working faces under similar geological conditions, and the scope of application is limited [11]. Similarity simulation [12] is an effective method to acquire predicted parameters. A formation model, resembling the actual project, is constructed in the laboratory, according to the similarity principle, and the predicted parameters are inferred by monitoring the changes of the model. This method has major advantages, such as intuition, simplicity, and short experimental period. However, the complexity of mining geological conditions, material strength, and human factors have a significant influence on the experimental results [13]. Currently, predicted parameters can be acquired at low cost through numerical simulation [14,15] due to the rapid progress of the computer technology. However, the results are more random because of the influence of the simulation unit and value parameters of rock formation [16]. In addition, based on a large amount of measured field data, several further methods have been proposed to determine predicted parameters, including neural network method [17], support vector machine [18], and genetic algorithm [19]. These methods provide additional ideas for determining predicted parameters, but none of them consider the essential causes of mining subsidence, i.e., the movement of the subterranean strata, and the influence of the strata distribution on the predicted parameters in the donga areas. The stratum control theory [20] predicts that the surface subsidence will change periodically with the periodic breaking of the main key stratum in the overlying rocks. As the complex stratum distribution conditions, such as burial depth and soil layer thickness on the surface, vary in the donga area, a significant change of the parameters is expected [21]. Thus, studying the influence of the rock formation movement and formation distribution on the predicted parameters is of major significance for determining reasonable probability integral parameters in loess donga.

Mining operations below loess donga areas can easily induce geological disasters, such as loess slope slippage and collapse, which not only aggravates the destruction of the ecological environment, but also increases the difficulty of mining subsidence prediction and environmental disaster assessment. According to the topographic features in the special areas, Guo et al. [22] subdivided mining subsidence disasters into three types: collapse disasters, slump disasters, and landslide disasters. Stead, D et al. [23] used the theory of fracture mechanics and plastic mechanics to study the mechanical mechanism of landslides caused by mining. Wang et al. [24] applied numerical simulations to study the stability evaluation model of the mining slope body and to refine the influence principles of the slope stability. Luo et al. [25] established a mathematical model of slope stability under the influence of longwall mining subsidence and constrained a value of 1.5 as the critical value of slope stability. Instead of studying the influence of the slope stability in the loess donga on the prediction of mining subsidence, the previous studies mainly focused on the stability of the mining slope. However, the additional mining slip caused by slope instability in the loess donga is the main source for uncertainties in the subsidence prediction.

Considering findings from previous studies, we propose a method for predicting surface subsidence in loess donga based on a probability integration model, which resolves the deficiencies of the conventional probability integration model in predicting subsidence. The method combines the probability integral model with the influence of rock movement and stratum distribution on estimated parameters and considers the slope stability and the influence of the slip additional displacement on the subsidence results. The plausibility of the method is evaluated by a field test.

#### **2. Prediction Method for Surface Subsidence of Coal Seam Mining in Loess Donga**

Various subsidence prediction models have been proposed to predict surface subsidence, including the probability integration model [8], the Weibull distribution model [21], and the influence function model [7]. With the rapid development of computer technology, numerical simulation technology [14–16] became progressively used in mining subsidence prediction. The probability integration model is the numerical simulation with the longest application time and the widest application range. Moreover, the probability integration method is numerical simulation of the subsidence prediction models with the longest application time and the broadest scope of reasonable application. Unfortunately, it is not suitable for the application in loess donga and unsatisfactory results are expected due to the influence of strata distribution and topography. Considering these problems, we therefore improved the application of the traditional probability integral model in the donga area. Our study will expand the application range of the probability integral model and also has important significance for solving the limitations of the model.

#### *2.1. The Probability Integration Model*

The basic principle of the probability integration model is to superimpose the subsidence probability of countless mining units (*ds*) to form the surface subsidence curve *w*(*x*) and the horizontal movement curve *u*(*x*). The integration model of the coal seam unit mining is shown in Figure 1 [21].

Once the coal seam inclination attains full mining, the values of the subsidence and of the horizontal movement of any point on the main section in the strike direction can be calculated by Equation (1) [7].

$$w(\mathbf{x}) = \frac{m\eta}{2} \left[ \text{erf}\left(\frac{\sqrt{\pi}x}{r\_0}\right) - \text{erf}\left(\frac{\sqrt{\pi}x - l + 2d}{r\_0}\right) \right]$$

$$u(\mathbf{x}) = bm\eta \begin{bmatrix} -\pi \frac{r\_0^2}{r\_0^2} & -\pi \frac{\left(x - l + 2d\right)^2}{r\_0^2} \\ \mathbf{e} & \mathbf{0} \end{bmatrix} \tag{1}$$

*m* is the mining height, *η* the subsidence coefficient, erf the error integral function, *r*<sup>0</sup> the mining influence radius, *r*<sup>0</sup> = *h*/tan*β*, *h* the buried depth of the coal seam, *β* the comprehensive moving angle, *l* the mining size, *d* the deviation of inflection point, and *b* is the horizontal movement coefficient.

Figure 1 and Equation (1) show that the subsidence coefficient, the horizontal movement coefficient, the comprehensive movement angle, and the deviation of inflection point are the four predicted parameters of the probability integration model. According to the calculation principle of Figure 1, the one-dimensional mining parameters are subjected to two-dimensional normalization processing. In order to simplify the calculation, based on the principle of two-dimensional lattice calculation, two-dimensional lattice calculation software of probability integration has been independently developed [26]. By inputting the predicted parameters and mining parameters into the calculation software, the predicted results of the subsidence basins of the traditional horizontal surface can be obtained. However, in loess donga, influenced by stratum structural factors, including burial depth, soil layer thickness, slope angle, and slope height, the predicted parameters of different areas of the surface are not only different, but also difficult to obtain. Furthermore, the predicted results often have large deviations, as they are affected by slope slippage.

#### *2.2. Prediction Process of Surface Subsidence in Loess Donga*

Considering the characteristics of the probability integration model and the limitations of prediction process of surface subsidence in loess donga, we propose a prediction method of surface mining subsidence based on the probability integration model. Figure 2 illustrates the application of the method. At first, the loess donga is subdivided into several horizontal stratigraphic regions with the same characteristics, in light of the stratigraphic distribution, according to the stratigraphic histogram and the surface contour map. Among the horizontal stratigraphic area, slope sub-areas were separated according to the features of the slope (angle and height). The parameters of the stratigraphic structure determine the values of the predicted parameters. The movement of the subterranean strata is the critical cause for the movement of the surface. According to the stratum distribution and the movement principle of rock-soil layers, we calculated the probability integration parameters of different horizontal stratum regions. Based on the superposition calculation principle of the probability integration model, we subsequently calculated the subsidence basins formed by different horizontal strata regions. Finally, the subsidence basin formed by the superposition of the horizontal area is corrected according to additional displacement of slope sub-region mining slip.

**Figure 2.** LDS application process.

#### **3. Determination of Probability Integration Parameters in LDS**

Various strata distribution and rock-soil movement principles lead to variations in the estimated subsidence parameters in differing regions. The subsidence coefficient is the key parameter of the probability integration model. The comprehensive movement angles, deviation of inflection points, and the horizontal movement coefficient control the influence range of the sinking basin and the horizontal movement value of the basin. The determination of the predicted parameters of the different regions is the main challenge of the coal seam mining in loess donga (LDS) prediction method. Previous studies have shown that the distribution and movement principles of rock-soil layers are closely related to the probability integration parameters [6]. Therefore, we studied the influences of the distribution and movement principles of the rock-soil layers on the probability integration parameters.

#### *3.1. Calculation of Subsidence Coefficient*

The key stratum theory in rock formation control implies that the movement of the underground rock formation is the critical cause of mining subsidence, and the mechanical state of the main key formation in the formation directly determines the migration principles of the overlying rock [20]. The main key stratum is located in the caving zone, the fracture zone or the bending subsidence zone of the mining overburden. The main stratum can be determined from the formation borehole histogram and the key stratum identification method [27]. The height of the caving zone (*h*c) and the height of the fracture zone (*h*f) of the mining overburden can be determined by various methods, including on-site measurement, empirical formula, engineering analogy, and theoretical calculations [28]. Based on the rock movement theory, we therefore propose a method for calculating the subsidence coefficient, according to the position of the main key stratum in the mining overburden.

Continuous deformation characteristics of the surface movement are the premise of predicting subsidence basins. If the height of the caving zone *h*<sup>c</sup> is larger than the height of the main key stratum Σ*h* (*h*<sup>c</sup> ≥ Σ*h*), the main key stratum is located in the caving zone. Figure 3 shows the discontinuous deformation characteristics of the surface when the main key layer is located in the caving zone. Under the influence of underground mining, the main key is completely broken, resulting in the loss of the support of the main key layer to the overlying rock and soil layer. Upward cracks formed by the rupture of rock and soil strata develop through the main key layers to the surface. Multiple vertical cracks cut the continuously deformed surface, showing the feature of stepped subsidence.

**Figure 3.** The main key stratum is located in the caving zone.

Therefore, we focus on the calculation method of the surface subsidence coefficient of continuous deformation for the condition that the main key layer is located in the curved subsidence zone or the fracture zone. Figure 4 illustrates the overlying strata movement principles of the continuous land surface deformation. In either of those cases, the main key stratum did not break, or was broken to form a stable masonry beam structure, which continued to undertake the supporting role of the overlying rock and soil strata. The upward cracks formed by the rupture of the rock and soil strata did not develop on the surface, and the surface curved and subsided to form a continuous subsidence basin. If the fracture zone height is smaller than *h* (*h*<sup>f</sup> ≤ Σ*h*), the main key stratum is in the bending subsidence zone (Figure 4a). In this case, the development height of the overburden mining fracture zone is below the main key stratum, not affecting the soil layer, and the overlying soil layer has no effect on the crack development. Based on the data of field measurements, Geweltzmann established Equation (2), which describes the link between the height of the fracture zone and the subsidence coefficient for mining a single horizontal coal seam by fallen method [29].

$$h\_{\rm f}^2 = \frac{7.25 m \eta}{\left(\cot \beta\_{\rm r} + \cot \psi\right)^2 K\_{\rm r}} \tag{2}$$

*ψ* is the full mining angle, which is related to the full mining degree (usually 55–59◦) and *K*<sup>r</sup> is the limit curvature of the top rock formation. The value of *K*r is related to the ratio *A* of the clay rock (mudstone and argillaceous sandstone) in the rock formation [30]. According to the ratio of the thickness of the bedrock constituted by the two rock types, *K*r is calculated by Equation (3).

$$K\_{\rm r} = 0.002 + 0.04A \tag{3}$$

If *h*<sup>c</sup> < Σ*h*< *h*<sup>f</sup> ≤ *h*, the main key stratum is located in the fracture zone. Figure 4b shows that the fracture zone is well developed through the bedrock and continues upward. The development height of the fracture zone is influenced by the thickness and properties of the soil layer. Since the clay layer thickness is not considered in Equation (2), the ultimate curvature *K*<sup>0</sup> (*K*<sup>0</sup> = *K*<sup>r</sup> + Δ*K*r) of the clay layer is introduced into Equation (4) to replace *K*<sup>r</sup> in Equation (2), and the soil layer movement angle is used to calculate the height of crack development [31].

$$h\_{\rm f}^2 = \frac{7.25 m \eta}{\left(\cot \beta\_{\rm s} + \cot \psi\right)^2 \left(K\_{\rm r} + \Delta K\_{\rm r}\right)}\tag{4}$$

In Equation (4), Δ*K*<sup>r</sup> is the ultimate curvature increment of the soil layer, Δ*K*<sup>r</sup> = 0.1 *h*s.

**Figure 4.** Migration principles of overlying rock under continuous deformation: (**a**) the main key stratum is located in the curved subsidence zone; and (**b**) the main key stratum is located in the bending zone.

The calculation method of the subsidence coefficient, derived from Equation (2) and Equation (4), is shown in Equation (5) in Table 1 for the scenario when the main key stratum is located at different positions.

**Table 1.** The subsidence coefficient of the main key stratum at different positions.


#### *3.2. Comprehensive Movement Angle, Horizontal Movement Coefficient, and Deviation of Inflection Point*

The size of the comprehensive movement angle is mainly related to the thickness of the rock (soil) layer and the rock movement parameters (Figure 1). The comprehensive moving angle is calculated from Equation (6):

$$\beta = \arctan\left(\frac{h\tan\beta\_{\rm s}\tan\beta\_{\rm r}}{h\_{\rm s}\tan\beta\_{\rm r} + h\_{\rm r}\tan\beta\_{\rm s}}\right) \tag{6}$$

with *h*<sup>s</sup> as the thickness of the soil layer, *h*<sup>r</sup> as the thickness of the rock layer, *β*<sup>s</sup> as the displacement angle of the soil layer, and *β*<sup>r</sup> as the displacement angle of the rock layer. *β*<sup>s</sup> and *β*<sup>r</sup> can be selected according to Tables 2 and 3 [21].


**Table 2.** Displacement angle (*β*s) of loose layer.

**Table 3.** Displacement angle (*β*r) of strata.


For the horizontal movement, coefficient *b* values of 0.2–0.3 are usually applied for China's coal mining fields. By relating the ratio of the thickness of the clay layer to the mining depth, *b* is calculated from Equation (7) [21].

$$b = 0.3 - 0.1\frac{h\_{\rm s}}{h} \tag{7}$$

The deviation of inflection point d is related to *f*, *l*, and *h*, *d* of the near-horizontal ore seam is calculated by Equation (8) in Table 4 [21].

**Table 4.** Displacement angle (*β*s) of loose layer.


Combined with the regional strata occurrence characteristics and mining rock movement principles (development height of caving zone and fracture zone) in loess donga, probabilistic integration parameters can be calculated from Equations (5)–(8) for each horizontal stratigraphic region. Subsequently, input of the constrained block parameters and predicted parameters of the subdivided horizontal stratigraphic areas into a selfdevelopment probability integration calculation software [26] enables the recognition of subsidence basins in the horizontal stratigraphic area after calculation.

#### **4. Slope Stability and Slip Principles in Loess Donga**

The subsidence basin calculated according to the superposition of the regional probability integration parameters of different horizontal strata only considers the influence of the horizontal strata in dongas, but excludes the effect of the mining-slip of the slope body. The main characteristics of the Yushenfu mining field (China): a plenitude of nearly horizontal coal seams, thick loose layers, high mining height, thin rock layers, and widespread distribution on the land surface of loess dongas. In the loess donga, the height of the slope body is usually 10–20 m, and the angle 10–30◦ [1]. Under specific circumstances, coal seam mining can easily lead to mining slip and instability of slopes, resulting in the different principles between mining subsidence and horizontal surface [6]. We studied the stability of the slope body and the additional displacement of mining slip to decipher the principles of the influence of the slope body on the subsidence in the loess gully area.

#### *4.1. Analysis of Slope Stability in Loess Donga*

Currently, the limit equilibrium method, the basic method for stability analysis of soil slope, is realized by calculating the safety factor of the slope to analyze the stability of the slope. Based on the Swedish strip method of Fellenius, we extended a series of other methods, including the Bishop method, the Janbu method, the Morgenster-Price method, and the Spencer method for the calculations [22,23,32].

We combined the analysis of the stability of the slope in the loess donga with the characteristics of the strip method. The slope is regarded as a uniform soil slope and its shear strength follows the Mohr-Coulomb failure criterion. The influence of the force between the soil strips on the stability of the soil slope is excluded. To simplify the calculation, we assume that the slip line of the soil slope is an arc. Figure 5 shows the mechanical model of the loess donga slope.

**Figure 5.** Mechanical model of slope body in loess dongas.

Figure 5 shows that the sliding slope oAB is approximately divided into parallelogram abcd vertical soil strips with number *n*. The safety factor of the slope in loess donga is set as *K*S, which is the ratio of the sliding moment *M*<sup>S</sup> to the anti-sliding moment *M*T. The sliding moment is caused by the shear component induced by the gravity of the soil strip. The anti-slip strength of the soil strip generates anti-slip moments. According to the Mohr-Coulomb strength theory, *K*<sup>S</sup> can be calculated from Equation (9) [32].

$$K\_{\rm S} = \sum\_{i=1}^{n} \frac{M\_{\rm S\bar{i}}}{M\_{\rm Ti}} = \sum\_{i=1}^{n} \frac{r \left(\gamma\_{\rm si} \mathbf{x}\_{i} h\_{\rm Pi} \tan \phi\_{\rm i\bar{i}} \cos a\_{\rm i} + c\_{\rm i\bar{i}} l\_{\rm i} \right)}{r \gamma\_{\rm si} \mathbf{x}\_{i} h\_{\rm pi} \sin a\_{\rm i}} = \sum\_{i=1}^{n} \frac{\gamma\_{\rm si} \mathbf{x}\_{i} h\_{\rm i} \tan \phi\_{\rm i\bar{i}} \cos a\_{\rm i} + c\_{\rm i\bar{i}} l\_{\rm i}}{\gamma\_{\rm si} \mathbf{x}\_{i} h\_{\rm pi} \sin a\_{\rm i}} \tag{9}$$

with *n* as the number of soil strips, r as the radius of the slip line, *γ*si as the bulk density of the soil strip, *x*<sup>i</sup> as the width of the soil strip, *h*pi as the height of the soil strip, *α*<sup>i</sup> as the angle between the direction of the gravity of the soil strip and the normal force *F*Ni, *φ*si as the internal friction angle of the soil strip, *c*si as the cohesive force of the soil strip, and *l*<sup>i</sup> as the arc length of the soil strip slip line. Since the soil strip element is *x*<sup>i</sup> ≈ *l*<sup>i</sup> and the soil slope is a homogeneous soil mass, the soil strip parameters *γ*si, *h*pi, *α*i, *φ*si, and *c*si can be converted into the bulk density *γ*<sup>s</sup> of the loess slope, the height of the loess slope *h*p, the angle of the loess slope *δ*p, the internal friction angle *φ*s, and soil cohesion *c*s. Therefore, Equation (9) is further simplified, and Equation (10) is obtained:

$$\mathcal{K}\_{\rm S} = \frac{2\left(\gamma\_{\rm s} h\_{\rm P} \cos^{2} \delta\_{\rm P} \tan \phi\_{\rm s} + c\_{\rm s}\right)}{\gamma\_{\rm s} h\_{\rm P} \sin 2\delta\_{\rm P}} \tag{10}$$

Equation (10) documents for a height of the slope, a bulk density of the slope or an angle of the slope of 0, the denominator becomes 0 and the equation meaningless. To analyze the relationship between slope stability and influencing factors in the loess donga, the slope instability coefficient *G*s was constrained in the experiment [6]. An increase in *G*s reduces the stability of the slope, and the risk of slippage caused by the instability of the slope increases. *G*s is calculated from Equation (11).

$$G\_{\rm S} = \frac{1}{K\_{\rm S}} = \frac{\gamma\_s h\_{\rm P} \sin 2\delta\_{\rm P}}{2\left(\gamma\_s h\_{\rm P} \cos^2 \delta\_{\rm P} \tan \phi\_{\rm P} + c\_{\rm P}\right)}\tag{11}$$

Equation (11) shows that the slope height, slope angle, soil bulk density, soil cohesion and soil internal friction angle are the main influencing factors of slope stability. To analyze the extent of the influence of every factor on the slope stability, the five factors are assumed to be independent of each other. Previous studies have shown that for *h*<sup>p</sup> > 30 m or *δ*<sup>p</sup> > 50◦, the slope will show collapse-type failure [32]. Therefore, we studied a range of the slope height from 0–30 m and a range of slope angle from 0–50◦. The test results show that the bulk density of the loess layer ranges from 16,300–18,600 N/m3, the cohesive force from 38,000–101,000 Pa, and the internal friction angle from 27.9–33.8◦ [29]. We studied the relationship between those factors and the stability of the slope by controlling the variable range of a single factor.

Figure 6a–c shows fitting degrees of 0.9860 for the *G*-*h*<sup>p</sup> power index function curve, of 0.9988 for the *G*-*δ*<sup>p</sup> proportional curve, and of 0.9992 for the *G*-*γ*<sup>s</sup> linear curve. *G* is positively correlated with *h*p, *δ*p, and *γ*s. With the increase in slope height, the slope angle, and the bulk density, the increasing speed of the slope sliding component force exceeds that of the anti-sliding component force, and the degree of the slope instability increases continuously. Figure 6d,e shows a fitting degree of 0.9998 for the *G*-*c*<sup>s</sup> linear curve, of 0.9994 for the *G*-*φ*<sup>s</sup> linear, and negative correlation of *G* with *c*<sup>s</sup> and *φ*s. With the increase in soil cohesion and internal friction angle, the shear strength of the soil slope, the anti-sliding component, and the stability of the slope increase. The results of grey correlation degree and orthogonal test analysis show that the sensitivity of the parameters affecting the slope stability in the loess donga is ranked from large to small, and the order is *δ*<sup>p</sup> > *h*<sup>p</sup> > *γ*<sup>s</sup> > *φ*<sup>s</sup> > *c*<sup>s</sup> [1]. The main reason is that the slope angle and slope height in this area have a large variation range, while the fluctuation range of the soil physical and mechanical parameters is small, which has little influence on the slope stability [22].

In conclusion, the slope instability coefficient in loess donga is positively correlated with slope height, slope angle, and soil bulk density, but is negatively correlated with soil cohesion and internal friction angle. According to Equation (11), the extent of instability of the slope body in donga areas can be determined prior to the impact of the mining activity.

**Figure 6.** Relationship between different influencing factors and slope instability coefficient: (**a**) Slope height (*h*p), (**b**) Slope angle (*δ*p), (**c**) Soil bulk density (*γ*s), (**d**) Soil cohesion (*c*s), and (**e**) Internal friction angle of soil (*φ*s).

#### *4.2. Additional Displacement of Slope Mining Slip*

The reliable understanding of striking differences in the movement principles of the surface and the horizontal surface is critical for mining in donga areas. The main movement and deformation characteristics of surface mining subsidence are subsidence, horizontal movement, and the influence range of the mining activity. The model of the movement and deformation of the inclined surface in the horizontal mining seam is shown in Figure 7 [6].

**Figure 7.** Deformation model of inclined surface in horizontal seam mining.

Figure 7 indicates an asymmetry of the inclined surface subsidence curve *w* (*x*, *z*) and the horizontal movement curve *u* (*x*, *z*) for a horizontal mine with the size *l*. The surface movement deformation curve is related to the stability of the slope, but also to the mining direction of the coal seam [21]. Therefore, the research progress should be combined with the influence of the stability of the mining slope in loess donga on the surface movement and deformation. Assuming that the slope slip is *R*(*x*), the additional subsidence slip Δ*w*(*x*) and the additional horizontal displacement slip Δ*u*(*x*) caused by the slope slip can be expressed in Equation (12). <sup>Δ</sup>*w*(*x*) <sup>=</sup> *<sup>R</sup>*(*x*) sin *<sup>δ</sup>*<sup>p</sup>

$$\begin{aligned} \Delta w(\mathbf{x}) &= \mathcal{R}(\mathbf{x}) \sin \delta\_{\mathbf{P}} \\ \Delta u(\mathbf{x}) &= \mathcal{R}(\mathbf{x}) \cos \delta\_{\mathbf{P}} \end{aligned} \tag{12}$$

The slippage of the slope is related to the mining instability degree of the slope. If the slope is unstable, it tends to slide. The higher the degree of instability, the higher the affinity for sliding. The slope can be subdivided into a first-class slope, a second-class slope, and a third-class slope according to the hazard degree of the slope body after sliding damage. The normative standards for the safety factor of slope stability are summarized in Table 5 [32].

**Table 5.** Slope stability safety factor.


The loess donga is a highly fragile ecological environment, which is attributed to the first-class slope, and mining in this area is regarded as an earthquake condition. According to Table 5, the safety factor *K*<sup>S</sup> in this area is 1.15, and *G*<sup>S</sup> = 0.87 as calculated by Equation (11).

The results of the numerical simulation and the field test show that, apart from the stability of the slope itself, mining has a larger impact on the stability of the slope. Compared with mining along the slope, the slope body is more prone to instability if mining occurs on the reverse slope. Taking *K*<sup>C</sup> as the influence coefficient of the slope mining direction on the slope instability coefficient for mining along the slope *K*<sup>C</sup> is 1.5, but for mining against the slope *K*<sup>C</sup> becomes 2.0 [24]. The mining instability coefficient *G*CS of the slope is expressed by Equation (13).

$$\mathbf{G\_{CS}} = \mathbf{K\_C} \mathbf{G\_S} \tag{13}$$

To sum up, a value of 0.87 is constrained as the critical coefficient that defines the slip instability in loess donga. At *G*CS ≤ 0.87 the slope is stable, and no slip occurs. At 0.87 < *G*CS, the slope is slippery and unstable, and the slip displacement *R*(*x*) is calculated from Equation (14).

$$R(\mathbf{x}) = G\_{\rm CS} \left[ w(\mathbf{x}\_0) \sin \delta\_{\rm P} + u(\mathbf{x}\_0) \cos \delta\_{\rm P} \right] \tag{14}$$

with *w*(*x*0) and *u*(*x*0) as the initial subsidence and initial horizontal movement, respectively. From Equations (12) and (14), Equation (15) for the additional displacement of slope slip is derived.

$$\begin{aligned} w(\mathbf{x}) &= w(\mathbf{x}\_0) + \mathbf{G}\_{\text{CS}} \left[ w(\mathbf{x}\_0) \sin \delta\_\mathbb{P} + \mu(\mathbf{x}\_0) \cos \delta\_\mathbb{P} \right] \sin \delta\_\mathbb{P} \\\ u(\mathbf{x}) &= u(\mathbf{x}\_0) + \mathbf{G}\_{\text{CS}} \left[ w(\mathbf{x}\_0) \sin \delta\_\mathbb{P} + \mu(\mathbf{x}\_0) \cos \delta\_\mathbb{P} \right] \cos \delta\_\mathbb{P} \end{aligned} \tag{15}$$

To summarize, the mining instability coefficient of the slope is calculated according to the slope parameters and the mining direction of the coal seam in loess donga. For conditions with 0.87 < *G*CS the slope is slippery and unstable, and additional mining slip will be generated as the slope subsides. First, the predicted basin of horizontal formation subsidence is obtained according to the superposition calculation of the software, and the initial subsidence curve and initial horizontal movement curve of the surface are extracted from it. Secondly, combined with the initial subsidence value and horizontal movement value of the slope area that will generate slip, as well as *δ*<sup>p</sup> and *G*CS, the additional displacement of slope mining slip is calculated by Equation (14). Finally, the movement and deformation curve of the initial horizontal stratum area is corrected by the additional displacement mining slip, and the prediction of the surface subsidence in the loess donga is finally realized.

#### **5. Case Study**

We selected the N1114 working face as a case study to verify the predicted effect of LDS on surface mining subsidence in loess donga. The methods of LDS and the field test were used to predict the subsidence basin.

#### *5.1. Engineering Case*

The Ningtiaota Coal Mine is a typical mine in the loess donga of the Yushenfu mining field in China. The N1114 working face mines 1–2 near-horizontal coal seams with an average thickness of 1.85 m. The working face has a width of 245 m and a recoverable length of 1922 m. The working face is limited to north and south by the unmined N1116 working face and the N1112 working face, respectively, and to the west by the 135 m stopline coal pillar and to the east by the 60 m mine field boundary coal pillar. An overview of the strata occurrence of N1114 working face is summarized in Table 6, that is constrained by combining drilling data with physical and mechanical parameters [33]. The key layers of the overlying rock are determined from the key layer theory and the data in Table 6. A 13.70 m thick horizon of fine sandstone is the main key layer and an 8 m horizon of fine sandstone is the sub-key layer. A 16.10 m layer of fine sandstone is a thick hard rock stratum, which is broken along with the breaking of the main key layer.

To fulfill the complete mining demand, the actual measurement and LDS subsidence comparison area is the section of N1114 working face from west to east, 300 m distant from the mining stop line. Full extraction has been achieved prior to the mining of the target area. Therefore, we added an additional 150 m computing area for the LDS to simulate actual mining on site. According to the contour distribution of the N1114 working face in Figure 8a, the trend direction with obvious slope distribution in the study area is selected for research. The profile A–A is constructed along the direction of the center of the working face, as shown in Figure 8b.


**Table 6.** Stratigraphy of the N1114 working face with physical and mechanical parameters.

**Figure 8.** Overview of the N1114 working face: (**a**) work surface layout; and (**b**) A–A stratigraphic profile.

#### *5.2. LDS Prediction*

We modelled the surface of loess donga according to the LDS application process that is subdivided into following steps:

Step 1: Figure 8b shows that the coal seam mining area and overlying strata are subdivided into four horizontal areas (A, B, C, and D), according to the distribution characteristics of geotechnical layers. According to the slope characteristics of dongas, the slope is further divided into six slope sub-areas: A1, B1, B2, C1, C2, and D1.

Step 2: The probability integration parameters of the A, B, C, and D horizontal formation regions are calculated, respectively.

Field measurement shows that *h*<sup>c</sup> is 7.5 m, *h*<sup>f</sup> is 24.5 m, and 24.5 m = *h*<sup>f</sup> < Δ*h* = 25.5 m [33]. The main key layer is in the curved subsidence zone and the surface forms a continuously deformed subsidence basin. From Figure 8 and the data of Table 6 it follows that the ratio of clay rock to the rock layer is 0.236 and the average firmness coefficient of the rock layer is 3.56. The mining height m of the working face is 1.85 m, the width is 245 m, and the calculated length is 450 m. The calculated length of the working face is larger than 1.4 *h* and the full extraction angle is taken is 55◦. The results of the probability integral parameters of the individual horizontal formation areas, calculated from Equations (5)–(8), are summarized in Table 7.

**Table 7.** Calculated probability integration parameters of the individual areas.


Step 3: Based on the probability integration model, we insert the parameters of the different regions and the coordinates of the inflection point of the block in the self-developed two-dimensional lattice probability integration calculation software [26] and calculated the value of the subsidence basin according to the mining order. The initial subsidence isoline cloud and horizontal movement cloud chart for the N1114 working face is shown in Figure 9.

Step 4: Subsequently, we calculated the slope instability coefficients of the six slope sub-regions (A1, B1, B2, C1, C2, and D1). According to Table 6 and Figure 8, the soil cohesion is 59 kPa, the friction angle in the soil is 28.2◦, and the bulk density of the soil is 18,600 N/m3. The results of the calculations are shown in Table 8.

**Table 8.** Slope instability coefficient of slope sub-regions.


**Figure 9.** Probability integration calculation result: (**a**) subsidence isoline cloud; and (**b**) horizontal movement isoline cloud.

The *G*CS-B1 relationship of the B1 slope is 0.87 < *G*CS-B1 = 0.97 (Table 8) indicating that the B1 slope will generate instability and slip during the coal seam mining process. Therefore, we constructed a cross section along the strike direction in the center of the working face (Figure 9) and the surface curve of dongas is obtained in combination with Figure 9 Subsequently, we refined the sinking and horizontal movement curves according to the calculated data in Table 8, using Equation (15). Figure 10 shows the predicted surface movement curve in the strike direction of the N1114 working face.

Figure 10 shows that the trend of the surface subsidence and horizontal movement curve in loess donga essentially conforms to the principles of surface mining movement and deformation. The maximum subsidence value reaches 1096.23 mm, close to the middle of the goaf. The minimum value of the horizontal movement of the surface is in the center of the goaf. The surface subsidence values on both sides of the goaf are 505.86 mm and 566.81 mm, respectively, hence almost half the maximum subsidence. The maximum values of the horizontal movement of the ground surface are 270.25 mm and 274.09 mm, which are located above both sides of the goaf. Affected by loess dongas, the surface subsidence value in the center of the goaf ranges between 1039.62 mm to 1096.23 mm. The horizontal movement in the center of the goaf fluctuates around 0, and the floating range is 14.45 mm to 20.51 mm. The subsidence and horizontal movement correction curves of the dongas area show that the slope at B1 induces different degrees of additional slip. The additional displacement of subsidence and horizontal movement slip generated by the B1 slope are 16.97 mm to 33.88 mm and 0.64 mm to 5.97 mm, respectively.

**Figure 10.** Surface movement prediction curve in dongas: (**a**) sinking curve; and (**b**) horizontal displacement curve.

#### *5.3. The Field Test*

A high-precision total station was used in the target research area for the field test to evaluate the application effect of LDS and refine the method. At first, two measurement control points were installed, and subsequently, the survey line A was arranged along the direction of the research area. From the west side, it was arranged 100 m away from the mining stop line of N1114 working face, including the 300 m research area. A total of 41 measuring points were arranged in the area. The distance between the measuring points was 10 m and the monitoring distance was 400 m. Since the mining speed of the working face was 20 m/d [33], until the mining of the working face was completed, the monitoring was performed five times once a week. After the mining stop of the working face, according to the change rate of the subsidence value, the monitoring time was 1 to 3 months and the monitoring time was five times until the surface movement stopped. The overall monitoring time was 1a, with a total of 10 monitoring times [21]. Figure 11 summarizes the measured curve of line A and the predicted curve of LDS in the target area.

**Figure 11.** Surface movement deformation curve of mining under loess donga: (**a**) sinking curve, (**b**) horizontal movement curve, and (**c**) legend.

Figure 11a shows that the interval between the 8th and 10th monitoring is more than 6 months. The fluctuation range of the subsidence value between the measuring points is 0.96–28.56 mm. The surface mobile basin is considered to be stable if the subsidence increment is less than 30 mm for six consecutive months [21]. The 10 monitoring results basically conform to the principles of surface movement and deformation, indicating that the monitoring is reasonable and effective. Comparing the measurement point data A11–A41 (Figure 11a) shows that the theoretically predicted subsidence value in the study area deviates from the actual subsidence value by 0.74–77.07 mm, with an average of 7.11 mm. The relative deviation is 0.1–12.9%, with an average of 1.6%. The deviation between the theoretical predicted horizontal movement value and the actual horizontal movement value is 0.14–17.07 mm (Figure 11b), with an average of 2.96 mm. The relative deviation ranged from 0.8% to 22.6%, with an average of 11.3%. Excluding the measurement errors of individual discrete points shows a high extent of agreement between the predicted curve and the actual subsidence curve. Hence, the requirements of the engineering are fulfilled.

The theoretically predicted additional subsidence slip in the study area A17–A25 is 16.97–33.88 mm and the theoretically predicted additional displacement of horizontal movement slip is 0.64–5.97 mm. Figure 11a shows that, based on the average maximum subsidence of 1075.26 mm in the center of the goaf, the actual subsidence value at A17–A25 increases between 3.51–71.16 mm, with an average of 30.59 mm, which is in line with the expected range. The relative deviation between the theoretically predicted subsidence curve and the actual subsidence curve is 0.2–4.4%, with an average of 1.0%. The actual sinking horizontal movement slip increment cannot be calculated since a benchmark for the horizontal movement is lacking, Figure 11b shows that the relative deviation between the theoretical predicted horizontal shift curve and the actual horizontal shift curve is 2.7–22.6% at A17–A25, with an average of 9.8%. The data document that LDS can effectively predict the subsidence and horizontal movement of the slope due to slip in loess donga.

#### **6. Discussion**

The study of surface subsidence in loess donga of coal mines is of significance for ecological reclamation and mining damage assessment. However, the probability integral model, which is widely used in the prediction of horizontal surface subsidence, is unsuitable because of the influence of the surface of the donga area. To acquire probability integration parameters with a high reliability is the main challenge. However, additional displacement of slip, caused by the instability in loess donga, reduces the prediction accuracy of the probability integration model. Therefore, previous studies have focused on the probability integral parameters of the horizontal surface and on the mining stability of the slope in the donga area. Implementing findings from previous studies, we here propose a new method for predicting coal seam mining subsidence in loess donga from the probability integration model that is based on the distribution and movement principles of rock-soil layers and considers the mining-slip effect of the slope. This predicting method is suitable for subsidence forecasting for mining of surface-level ore seams in most donga areas. Due to the thickness of the sloping ore seam, the mining depth, and the offset caused by the sloping ore seam, the surface mining movement deformation basin is asymmetrical. Under the superposition of the surface in the donga area, the LDS prediction method is no longer applicable. In the future, we will study the impact of sloping seam mining on LDS prediction that will further refine the LDS subsidence prediction method.

#### **7. Conclusions**

We improved the subsidence prediction process of the probability integration model in loess donga areas to solve the problem that surface mining subsidence in such areas is difficult to predict. The following main conclusions are drawn from our new data:


the help of the probability integration software, the subsidence basin formed by the superposition of multiple horizontal stratigraphic regions can easily be obtained.


**Author Contributions:** The present article is addressed by the eight authors mentioned, each of whom were responsible for various aspects of the work. Conceptualization, Y.G. and B.Z.; methodology, Y.G. and B.Z.; software, B.Z.; validation, Y.G. and B.Z.; formal analysis, Y.G., Y.H., and B.Z.; investigation, Y.G., B.Z., X.M., D.Z. (Di Zhai) and J.W.; resources, B.Z.; data curation, Y.G., D.Z. (Defu Zhu) and B.Z.; writing—original draft preparation, Y.G. and B.Z.; writing—review and editing, Y.G., B.Z., X.M., D.Z. (Di Zhai), D.Z. (Defu Zhu), Y.H., Z.S. and J.W.; visualization, Y.G.; supervision, B.Z., X.M., D.Z. (Di Zhai), D.Z. (Defu Zhu), Y.H., Z.S. and J.W.; project administration, B.Z. and Y.G.; funding acquisition, B.Z., Y.G., D.Z. (Di Zhai), Y.H. and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China through (Grant No. 51874230 and 52074208). And Shanxi Applied Basic Research Programs, Science and Technology Foundation for Youths (Grant No. 202103021223073).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

