Abstract
Underground coal mining will cause large-scale surrounding rock movement, resulting in surface subsidence and irreversible deformation of surface morphology, which would lead to geological disasters and ecological environment problems. In this paper, FLAC3D numerical model is built based on the natural slope gradient, slope type, and included angle between the slope and working face, and their influences on the change of surface morphology and soil erosion caused by underground coal mining is studied. Research results show that the change of slope gradient caused by underground mining decreases with the increase of natural slope gradient, while slope length has opposite laws; different slope types have different changes of slope morphology. The order of slope types corresponding to gradient changes is mixed slope < uniform slope < concave slope < convex slope; the length of the concave and uniform slope decreases, and the convex and mixed slope length increases. When the included angle between the slope and working face is 0° ≤ α < 90°, the underground mining will cause the natural slope gradient increase, the change of the slope gradient will increase with the rise of the angle, the slope length will decrease, and the rate of decrease will be reduced with the increase of the angle. Coal mining will cause the increasing of the runoff and erosion modulus of slope, mainly runoff modulus.
1. Introduction
As an important nonrenewable energy source, coal is the main energy supplier of domestic industrial production, which has a direct effect on China’s economic development [1–3]. In recent years, however, the emergence of new energy sources has greatly reduced the demand for coal in China, but it still cannot be completely replaced. Coal reserves are large and buried deep with well mining as the main mining method [4–6] in the country. Underground coal mining can inevitably cause large-scale surrounding rock movement, which would transfer to the surface layer by layer with the advance of coal mining [7, 8], resulting in changes in surface morphology, including slope gradient and slope length, and eventually they are displayed in the way of subsidence, which would lead to changes in physical and chemical properties of the soil, decline in vegetation coverage, soil erosion, and a series of environmental problems [9, 10]. Therefore, it is of great significance to study the change of surface morphology caused by underground coal mining [11, 12]. Researchers in both home and abroad have carried out a lot of research studies on the influence of underground coal mining on surface morphology. Zhang et al. [13, 14] studied the influence of mining depth, mining thickness, natural slope of surface, rock-soil ratio, etc. on the slope morphology of the mining subsidence area by numerical simulation from the perspective of underground mining conditions and probed into effects of the aforementioned factors on soil erosion from a new perspective of combining underground with surface. Liu and Zhao [15] used FLAC3D simulation software to investigate the influence of underground mining on irregular slope morphology and soil erosion. It was found that different slope types had different responses to slope deformation caused by underground mining. Yu and Liu [16] studied the law of surface movement in the wide strip mining of extremely thick loess seam by laying monitoring points on the surface of a mining area and obtained reasonable strip mining parameters by using probability integral and Wilson theory. The influence of underground mining on slope morphology is not only affected by mining conditions, but also by the original surface morphology and its relative position with the working face. Therefore, from the perspective of original surface morphology and its relative position with the working face, MIDAS-FLAC3D numerical simulation software is selected to study the influence of natural slope gradient, slope type, and the angle between the slope and working face on the change of slope morphology, including slope gradient and slope length, caused by underground coal mining in order to provide reference and basis for reducing the influence of underground coal mining on the slope.
2. General Situation of Research Area
We choose Dafosi Coal Mine as the research object. The mine is located at the border between Binzhou City and Changwu County, Shaanxi Province, 12 km southwest of the downtown of Binzhou City. Geographic coordinates: eastern longitude 107°49′40″ and northern latitude 35°05′00″. Dafosi Mine Field covers an area of 86.3 km2. The main coal-bearing strata in the coal mine are Yan’an Formation, and 4# coal is the main coal seam in this group. The dip angle is very small, most of which are between 3° and 5°. The mining conditions are favorable with the average thickness of 11.5 m, mining depth 400–600 m, mining thickness 2.8–12.5 m, the width of the working face 200–220 m, and annual advancing distance 1-2 km. A fully mechanized caving method is adopted, and the method is all caving method, and the roof caving method is used for the roof management.
3. FLAC3D Model Building
The model framework is designed based on the surface slope morphology of Dafosi coal mine in the Binzhou mining area and the characteristics of coal mining. Three-dimensional geological model is shown in Figure 1.

3.1. Rock Layer Design
Based on borehole information in Dafosi Coal Mine, the underground rock formation is designed. The main coal seam of the mine is 4# coal seam. In order to be conducive to calculation and modeling, the strata are merged according to their mechanical properties. Finally, there are 24 overlying strata (including coal seam and floor) in 4# coal seam. The main types of strata are sandy mudstone, sandstone (coarse, medium, and fine sandstone), conglomerate, and loess. The physical and mechanical properties of each stratum are shown in Table 1.
3.2. Design of Surface Morphology
Slope gradient, slope type, and the included angle between the slope and working face are taken as variables to design the surface morphology, while the underground strata, mining depth, mining thickness, working face width, and annual advancing distance remain unchanged in the model. The horizontal projection of slope is located at model X = (100–300 m).(1)Different Slope Gradients. Five FLAC3D models with natural slope gradients of 0°, 5°, 10°, 15°, 20°, and 30°, slope length of 150.573 m, maximum mining depth of 400 m, and minimum mining depth of 324.713 m were built.(2)Different Slope Types. Four FLAC3D models of straight slope, concave slope, convex slope, and mixed slope were constructed. The average depth is 400 m, and the slope gradient is 15°.(3)Different Included Angles between the Slope and Working Face. Five FLAC3D numerical models with the slope gradient of 15°, slope length of 207.055 m, maximum mining depth of 400 m, and minimum mining depth of 346.41 m were established. The angles between the slope and working face were 0°, 30°, 45°, 60°, and 90°, respectively.
3.3. Numerical Model Design
The average mining depth of 4# Coal Seam is 400 m, mining thickness is 9 m, annual progressive distance is 1000 m, and working face width is 200 m, that is, X = 600 m, Y = 200 m, and Z = 400 m.
3.4. Boundary Conditions and Initial Forms
x-axis direction is the subjective direction of our survey. The two boundary surfaces in X- and Y-direction of the model are set as single boundaries, the bottom boundaries in Z-direction as full boundaries, and the upper boundaries as free boundaries [12, 17]. Initially, no external force is applied to the stratum, which is only affected by its own gravity. The initial in situ stress field is shown in Figure 2.

3.5. Monitoring Point Settings
A series of monitoring points are set up along the working face to real-time monitor the change of slope displacement caused by coal mining. The mining is stopped when the fully mining is achieved, and the real-time monitoring values of points are extracted using “hist write” command flow.
3.6. Mining Methods
The method of step-by-step excavation is adopted in the mining. The distance of each step is 100 m, and the mining stops when the fully mining achieved. Before the simulation, a pilot simulation experiment is carried out. The results show that the fully mining is reached when the excavation advances 500 m; thus, the length, width, and height of the excavation size are 500, 200, and 9 m, respectively. The dark part in Figure 3 is the excavation area.

4. Results and Discussions
4.1. Effect of Natural Slope Gradient on Slope Morphology
The cloud pictures of X- and Z-direction displacement caused by underground mining under fully mining conditions are shown in Figure 4.

(a)

(b)
It can be seen from Figure 5(a) that when fully mining is achieved, the value of slope subsidence increases with the rise of the relative distance between the monitoring point and slope top, and the relationship between subsidence and slope gradient is inversely proportional. However, the relationship will change locally when the relative distance between the monitoring point and the top of the slope is less than 43 m. At the same time, the subsidence decreases with the increase of the relative distance, and the subsidence is proportional to the slope gradient. There are two possible reasons for the local change. One is that the nature slope gradient is different. When the surface morphology is horizontal, no such local change takes place. There is a positive relationship between the subsidence and the relative distance between the monitoring point and the top of the slope. With the increase of the slope gradient, the local change becomes more and more obvious. Second reason is that the position of the slope top is coincidence with the open-off cut. The aforementioned two reasons would at the same time affect the stability of the slope top. The greater the natural slope gradient, the worse the stability of the slope, and the slope has a more sensitive response to the interference of underground mining. From Figure 5(b), it can be concluded that with the increase of natural slope gradient, the horizontal displacement of slope caused by underground mining increases. With the increase of the relative distance between the monitoring point and slope top, the horizontal displacement of the slope increases first and then decreases. The maximum horizontal displacement occurs about 170 m from the slope top, that is, at the central position of goaf. Thus, underground mining has a variable influence on values of the subsidence and horizontal displacement for different slope gradients. It has a greater influence on both the horizontal displacement of the slope with a larger slope gradient and the subsidence of the slope with a smaller slope gradient.

(a)

(b)
4.1.1. Effect on Slope Gradient
By calculating the data extracted from monitoring points of different slope gradients, the variation of the slope gradient and length can be obtained, as shown in Table 2.
The variation of the slope gradient and length of different natural slope gradient caused by underground coal mining is shown in Figure 6. No matter the size of natural slope gradient, mining would cause the increase of slope gradient, but the increase rate is different. With the increase of natural slope gradient, the increase of the gradient caused by the mining gradually decreases; the gradient increases from 0° to 30° and the increment decreases from 1.417° to 0.583°. It shows that underground mining has little influence on the slope gradient with larger natural slope gradient. The reason is that when the natural gradient of the slope is larger, the lateral stress is higher, which results in a higher tensile deformation, then the slope deformation shows that the increase of slope gradient decreases.

(a)

(b)
4.1.2. Effect on Slope Length
When the natural slope gradient is greater than 5°, the underground mining would cause the increase of slope length. The change of slope length is proportional to the natural slope gradient. The slope gradient increases from 5° to 30°, and the slope length increases from 0.164 to 0.524 m. The latter is 3.2 times of the former. When the natural gradient of slope is less than 5°, underground mining will cause the decrease of slope length. With the gradual increase of slope gradient, the reduction of slope length decreases. It can be seen from Figure 6(b) that the natural slope gradient is in the range of 5°–10°, the change of slope length caused by the mining has a zero value point and the slope length does not change at the point. Thus, underground mining has a greater influence on the length of the slope at a smaller slope gradient, which is caused by that the response of the lateral stress of the slope to the mining is gradually strengthened with the increase of the slope gradient, while the vertical stress is gradually decreased.
4.2. Influence of Slope Type on Slope Morphology
Slope type will have a significant influence on the change of surface morphology caused by underground mining. We have found from the related literatures [18] that the actual slope type is simplified as uniform slope by researchers in the process of study, in which there must be certain errors in the research results, that is, the estimation topographic factors for concave slope is higher than the actual value, and the estimation of topographic factors for both the convex slope and the mixed slope is lower than the actual value. In order to reduce errors caused by the irregular slope, Castro and Zobeck [19] divide the slope into several segments and assume that the slope of each segment is uniform. Undoubtedly, the calculation errors could be reduced by the method, but it is still difficult to accurately calculate the slope length and gradient of the irregular slope due to the complexity and variability of natural slope shape. In view of the problems in the abovementioned method, we propose a method of replacing irregular slopes with uniform slopes under the same conditions (topography, soil, etc.) in the same region to obtain the replacement length and gradient of irregular slopes. Because the irregular slope has the same erosion energy as that of the equivalent uniform rainfall, so the proposed method is feasible in practice. In this paper, the same substitution method is adopted to replace the irregular slope with the uniform slope, thus the length and gradient of the irregular slope can be substituted.
The sketch of irregular slope type is shown in Figure 7, and the maximum subsidence corresponding to different slope types is shown in Table 3. The cloud pictures of X- and Z-direction displacement caused by underground mining under fully mining conditions are shown in Figure 8 (taking the mixed slope as an example).


(a)

(b)
The fully mining can be achieved when the workface of underground coal mining advances forward by 500 m. The maximum subsidence of slopes with different slope types is shown in Table 3. The order of slope types corresponding to the maximum subsidence is convex slope > uniform slope > mixed slope > concave slope. The maximum subsidence occurs at the foot of slope.
The subsidence and horizontal displacement of different slope types are shown in Figure 9. Regardless of the slope type, the value of slope subsidence increases with the relative distance between the monitoring point and the top of the slope. The order of slope types corresponding to the subsidence value is concave slope < mixed slope < uniform slope < convex slope, and the difference of slope subsidence value between different slope types is small. But when the relative distance between the monitoring point and the top of the slope is less than 65 m, the relationship will change locally. At the same time, the value of slope subsidence varies greatly. The order of slope types corresponding to the value is convex slope < uniform slope < concave slope < mixed slope. Regardless of the slope type, the change of horizontal displacement caused by underground mining increases first and then decreases with the increase of the relative distance. When the relative distance is about 65 m, the horizontal displacement reaches its maximum value. Horizontal displacement varies with slope type, showing a low contrary to subsidence value. The order of slope types corresponding to the horizontal displacement is convex slope > uniform slope > mixed slope > concave slope. When the relative distance between the monitoring point and the slope top is less than 65 m, concave slope will change locally. Horizontal displacement of the concave slope is obviously larger than that of the mixed slope, but this kind of locality change does not affect overall trends.

(a)

(b)
4.2.1. Effect on Slope Gradient
The method of uniform slope substitution is adopted for the calculation of both irregular slope gradient and slope length under the same conditions, that is, the slope length and gradient of irregular slope are replaced by uniform slope length and gradient. The substitution principle is that the potential energy of runoff generated on irregular slope during the same rainfall process is equal to that generated on the alternative uniform slope. That is,
Therefore,where and indicate the potential energy of runoff of the irregular slope and the uniform slope in the same rainfall process, respectively; and indicate rainfall intensity and soil infiltration rate at any point on a single wide surface, respectively; indicate rain density; indicate gravity acceleration; indicate curve equation of longitudinal profile of irregular slope; and indicate linear equation of longitudinal profile of replacing uniform slope.
Taking the mixed slope as an example (shown in Figure 10), in the same rainfall process the runoff potential energy generated on the mixed slope “ab” and the uniform slope “bc” are as follows:where and indicate potential energy of runoff of irregular slope “ab” and uniform slope “bc” in the same rainfall process, respectively; and indicate curve equation of “ab” and linear equation of “bc”, respectively; the rest are the same as formulas (1) and (2).

From the definite integral principle, it can be seen that is the area (Aab) of the surface surrounded by the curve “ab” and the x-and y-axes, and is the area (Abc) of the surface surrounded by the straight line “bc” and the x- and y-axes. In order to make the runoff potential energy of the mixed slope “ab” and the uniform slope “bc” be equal in the same rainfall process, it is necessary to make Aab = Abc. Therefore, the slope length and gradient of the equivalent uniform slope can be obtained by calculating the area of the curved surface. In this model, the integral formulas for calculating the area of concave, convex, and mixed slopes (the area of uniform slopes can be directly obtained) are shown in Table 4.
The replacement of slope length and gradient of four types of slopes before and after mining (concave, convex, mixed, and uniform slope) can be obtained, and the calculation results are shown in Table 5.
It can be seen from Table 5 that no matter what the slope type is, the underground mining will inevitably lead to the increase of natural slope gradient, but changes in the gradient vary with the slope type. The order of slope types corresponding to the increase of slope gradient is mixed slope < uniform slope < concave slope < convex slope. The greatest increase in the gradient of slope for the convex slope is 69 times that of the mixed slope with the smallest increase in the gradient of slope.
4.2.2. Effect on Slope Length
Table 5 shows that the influence of underground mining on slope length of different natural slope types is different. Among them, the length of the concave slope and uniform slope decreases by 1.109 and 0.513 m, respectively, and the length of the concave slope decreases by 2.16 times that of the uniform slope. Obviously, underground mining will cause an increase of slope length of both the convex slope and mixed slope, with the increase of 5.746 and 1.497 m, respectively. The increase of the length of the convex slope is 3.84 times that of the mixed slope.
4.3. Effect of Included Angle between Slope and Propulsion Direction on Slope Morphology
The maximum subsidence corresponding to different included angles are shown in Table 6. The cloud pictures of X- and Z-direction displacement caused by underground mining under fully mining conditions are shown in Figure 11.

(a)

(b)
It can be seen from Table 6 that when fully mining is achieved and the angle between the slope and the working face is 0° ≤ α ≤ 45°, the maximum subsidence value of surface slope caused by underground mining is increased with the increase of the angle. However, when the angle is 0°, the maximum subsidence value is much larger than that of other slope gradients. When the angle is 45° ≤ α < 90°, the maximum subsidence is decreased with the increase of the angle. The included angle of 90° is a special case, which does not conform to the aforementioned rules. Thus, the working face should avoid being parallel to the surface slope when it is laid out. If conditions permit, the vertical setting of the working face with the surface slope can reduce the value of slope subsidence caused by underground mining.
4.3.1. Effect on Slope Gradient
Under the conditions of underground mining, the influence of the angle between the slope and working face on the slope morphology is shown in Table 7, and the change of slope gradient and length with the included angle is shown in Figure 12.

(a)

(b)
When the included angle between the slope and working face is 0° ≤ α < 90°, the underground mining will cause the natural slope gradient to increase, and the change of slope gradient would rise with the increase of the angle. The slope gradient would decrease by only 5% when the angle is 90°, which can be neglected, and the relationship abovementioned between them is not established. The reason is, when the angle is 90°, the slope inclination is perpendicular to the working face and the vertical and horizontal displacement of the slope are basically same at each place along the inclination from the top to the foot of the slope, therefore, the slope shape changes little and the slope shows integral sinking. When the angle between the slope and working face is 0° ≤ α < 90°, the fitting equation between the change of slope gradient and the angle is y = 1.3428 sin x2 − 2.9812 sin x + 2.3402 (0° ≤ x < 90°, R2 = 0.998).
4.3.2. Effect on Slope Length
When the included angle between the slope and working face is 0° ≤ α < 90°, the change of slope length is decreased with the rise of the angle, and the decrease rate is also increased with the rise of the angle. Meanwhile, the fitting equation between the change of slope length and the angle is y = −5.0115 sin x2 − 5.2697 sin x + 11.775 (0° ≤ x < 90°, R2 = 0.999). When the angle between the slope and the working face is 90°, the slope length is decreased, but the reduction rate is 0.08 × 10−3, which can be neglected.
5. Influence of Underground Mining on Soil Erosion
Soil erosion is mainly affected by the slope gradient and length when other conditions (topography, soil, etc.) are relatively unchanged. The relationship between the slope gradient and erosion were studied extensively by numerous researchers [20, 21]. They believed that the relationship between slope gradient and erosion is a power exponential function and the quantity of soil erosion is increased with the increase of slope gradient. Similarly, the relationship between the slope length and soil erosion also shows a power exponential function, and the relationship is M = aLb. When b < 0, the amount of erosion is negatively correlated with slope length, and when b > 0, the amount of erosion is positively correlated with slope length [22]. Binzhou mining area belongs to the Loess Plateau area, which is dominated by hydraulic erosion. Therefore, rainfall factor is also an important influential factor on soil erosion. The research by Xie et al. shows that the rainfall standard of soil erosion on the Loess Plateau is 12 mm [23]. According to the research results of Wu et al. [22], on the relationship between topography and soil erosion in Loess Gully area, we chose 16.6 mm as rainfall condition to calculate the erosion and runoff modulus. The quantitative relationship between the slope gradient and length and the erosion and runoff modulus is shown in Table 8.
5.1. Effects of Different Natural Slope Gradient on Soil Erosion
Bring the values of the slope gradient and length in Table 2 into the relations of the erosion and runoff modulus in Table 8. The pre- and post-mining erosion and runoff modulus of different natural slope gradients are calculated under erosive rainfall conditions (see Table 9).
Underground mining will cause the increase of the erosion and runoff modulus on slope. With the increase of natural slope gradient, the erosion and runoff modulus is increased. However, the increase of the erosion and runoff modulus shows a decreasing trend. When the natural slope gradient is increased from 0° to 30°, the increment of erosion modulus decreases from 2.854 to 0.337 t km−2 which decreases by 85%, and the increment of runoff modulus decreases from 374.214 to 42.643 t km−2 which decreases by 88%, which indicates that the underground mining has great influence on the slope gradient with smaller natural slope gradient. The reason may be that the displacement value of any point in any direction caused by the mining can be decomposed into mining component and slip component. Similarly, the sand and runoff production of the slope could be decomposed into mining component and slip component. The increment of both sand production and runoff production caused by the mining would decrease with the increase of natural slope gradient, while the slip component will increase with the increase of natural slope gradient. However, we just studied the mining component instead of the slip component because the latter would lead to inconsistent results.
5.2. Influence of Slope Type on Soil Erosion
Regardless of the slope type, the underground mining would cause the increase of the mining component of the erosion and runoff modulus, and the type of slope is different which could lead to the quite different increment in the erosion and runoff modulus. The order of slope morphology corresponding to the increment of the erosion and runoff modulus is mixed slope < uniform slope < concave slope < convex slope. It can be seen from Table 10 that the influence of underground mining on erosion and runoff modulus of the mixed slope is less than that of other slope types, which may be due to the fact that the influence of mining on the mixed slope will offset each other in the concave and convex parts, resulting in the erosion and runoff modulus being significantly less than that of other slope types. The erosion and runoff modulus of the convex slope is the largest among the four types of slope. The erosion and runoff modulus of convex slope is 238 and 120 times that of the mixed slope, respectively. The reason for the obvious differences may be that different slope types bear different stress at different parts of the slope surface. Thus, in order to reduce the risk of soil erosion caused by underground mining, earth filling, and excavation should be carried out according to the actual situation before underground mining.
5.3. Effect of Included Angle between the Slope and Working Face on Soil Erosion
As shown in Table 11, when the angle between the slope and working face is 0° ≤ α < 90°, the mining components of the erosion runoff modulus caused by underground mining show an increasing trend with the increase of the angle, which also indicates a positive proportion relationship between the increase rate and the angle. When the angle is 90°, the mining component of sand and runoff production decreases and the reduction rate is 0.38% and 0.37%, respectively. Reducing the angle between the slope and working face can reduce the impact of mining on soil erosion. Therefore, when the working face is laid out, the included angle between the slope and working face should be reduced as much as possible.
6. Conclusions
(1)No matter what kind of surface morphology, underground coal mining would lead to the increase of slope gradient. The change of slope gradient is influenced by natural slope gradient, slope type, and the included angle between the slope and working face. When the slope gradient is between 0° and 30°, the underground mining disturbance to the slope is inversely proportional to the natural slope gradient. The sensitivity of underground mining response to different slope types are different, the order of slope types corresponding to the increase of slope gradient is mixed slope < uniform slope < concave slope < convex slope. The influence of underground mining on slope morphology is highly affected by the included angle between the slope and working face. When the angle is 0° ≤ α < 90°, the slope gradient rises as the angle increases, and its rate of increase rise rapidly.(2)Underground mining has less influence on the slope gradient with the larger natural slope gradient and less influence on the slope length with the smaller natural slope gradient. The gradient and length of the slope show a trend of increase and decrease, that is, one is rising, and the other is falling.(3)The influential regularity of underground mining on slope length is weak. When the natural slope gradient is less than 5°, the length of slope decreases, and the change of slope length is inversely proportional to the natural slope gradient. When the slope gradient is between 5° and 30°, the length of slope increases due to underground mining, and the increment of slope length is proportional to the natural slope gradient. The influence of underground mining on slope length varies with different types of slopes. The length of the uniform slope and concave slope decreases, while the length of the convex slope and mixed slope increases. The maximum change in the length of the convex slope is 5.746 m. The influence of the included angle between the slope and working face on the slope length should not be neglected. When the included angle is 0° ≤ α < 90°, the variation of the slope length decreases with the increase of the angle, and the slope is perpendicular to the working face with the slight decrease in the length of slope.(4)Underground coal mining will increase the erosion and runoff modulus of slope, mainly increasing the runoff modulus of slope. The increment of the erosion and runoff modulus is affected by the natural slope gradient, the slope type, and the included angle between the slope and working face. The soil erosion, runoff, and sand yield modulus are inversely proportional to the slope gradient, but underground mining would increase the risk of landslide with larger slope gradient. The order of slope morphology corresponding to the increment of the erosion and runoff modulus is mixed slope < uniform slope < concave slope < convex slope. When the included angle between the slope and working face is 0° ≤ α < 90°, the mining components of the sand and runoff production show an increasing trend with the increase of the included angle.
Data Availability
All relevant data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (grant number No. 414023080) and supported by the Special Research Project of Shaanxi Province (grant number No. 17JK05070) and the Project of Geological Research Institute for Coal Green Mining (No. MTy2019-01).