Abstract
In many mining areas, there is a lot of pressed coal under buildings, railways, and water bodies, and there are geological conditions of thick unconsolidated layers, so the surface subsidence has certain particularity. The key to solve the problem of pressed coal is the control technology of surface subsidence. The development of strip filling technology provides a kind of economical and effective surface subsidence control technology. With the passage of time, the natural weathering, flow, and fracture of some coal pillars may lead to the instability and failure of some coal pillars. Therefore, the selection of filling parameters plays an important role in the stability of coal pillars. In order to study the influence factors of surface movement and deformation, considering the influence factors of filling rate, key layer thickness, filling body strength, and unconsolidated layer thickness, FLAC 3D was used to simulate the surface movement deformation, and the orthogonal test method was used to analyze the simulation results, and the sensitivity of the main control factors affecting the surface deformation affected by strip filling mining was studied. The results show that the order of importance of the four factors on the four surface movement and deformation parameters is filling rate > key layer thickness > filling body strength > unconsolidated layer thickness. The influence of these four factors on the surface movement and deformation parameters is gradually decreasing, but the influence degree of different factors has drawn a certain gradient.
1. Introduction
The large-scale development and utilization of mineral resources have laid a solid foundation for the rapid and sustainable economic development. However, the development of resources has changed the original geological structure and occurrence environment of rock and soil in the mining area, causing great disturbance to the geological ecological environment, geomorphic pattern, and geological process in the mining area and bringing negative impact on human life and ecological environment. Coal mining causes strata movement and damage, induced rock burst, goaf collapse, underground debris flow, slope collapse and landslide, hydrogeological environment change, and other disasters [1, 2]; at the same time, it also causes water flow and discharge in rock strata and even causes land desertification. In addition, the discharge of coal mine wastewater causes heavy pollution to surface water and vegetation and makes the mining area land poor. It is necessary to implement “green mining” with low mining, high efficiency, and low emission [3, 4]. The geological conditions in the coal mine area are complex, and there are many kinds of geological disasters in the process of coal mining, such as roof disaster, ground collapse, gas explosion and escape, groundwater level drop and pollution, water inrush and sand burst, spontaneous combustion of coal seam and coal gangue dump, and secondary disasters, which are the main sources of environmental pollution in the mining area [5]. Filling mining is a green mining technology, which not only controls the movement and deformation of overburden rock but also makes full use of the solid waste generated in the process of coal mining and washing to protect the ecological environment of the mining area [6–8]. As a kind of mining technology with a high recovery rate, good safety, and more environmental protection, filling mining was born with the trend.
However, with the long-term uncontrolled mining of coal resources, in China, some old mining areas are facing the situation of exhaustion, and most of the coal resources that have not been mined belong to the “three under” (under buildings, railways, and water) pressed coal [9–12]. The key to solve the problem of “three under” pressed coal is the control technology of surface subsidence [13]. According to incomplete statistics, in China, there are more than 14 billion tons of “three under” coal reserves in production mines [14]. Among them, coal under buildings reaches 8.76 billion tons, and 65% of coal under buildings is located under villages. Particularly in North China, Central China, East China, and other mining areas, there are thick unconsolidated strata, such as Yanzhou mining area with the maximum thickness of 230 m, Yongxia mine with the maximum of 310 m in the area, 170 m in Hebi mining area, 270 m in Huaibei mining area, 460 m in Huainan mining area, and 300 m in Kailuan mining area. In the mining area with a thick loose layer, the surface subsidence law has certain particularity (when the loose layer contains water, the subsidence coefficient is relatively large). Through the practical research of mining workers over the years, the research on the overlying strata and surface movement and deformation caused by coal mining under the general geological conditions has been more sufficient, while the surface movement caused by mining under the thick loose layer has been studied. At present, there is still a lack of sufficient understanding of the law of shape.
In recent years, backfill coal mining [4, 15–17] has become a popular and widely applied mining method for the safe extraction of coal resources trapped under buildings, railways, and water bodies. For a long time, the main mining method of coal under villages is strip mining, which alleviates the problem of “three under” pressed coal to a certain extent. Many scholars have done a lot of research on strip mining [18–22]. However, the strip mining still faces some problems: first, the coal recovery rate is low, and nearly 50% of the strip coal pillars are left in the underground, which causes a lot of waste of resources and is not conducive to the sustainable development of the coal industry. Second, with the passage of time, some coal pillars naturally weathered, creep, and rupture, which may lead to the instability and failure of some coal pillars, and affect the adjacent coal pillars, resulting in “domino effect”; finally, the support system of the whole strip coal pillar is destroyed and unstable, resulting in a large area of goaf collapse, resulting in a hurricane, threatening the underground staff and production equipment, and even causing mine earthquake disasters.
In order to control the environmental damage and disasters caused by coal mining subsidence, it is an urgent problem to realize the green mining of coal resources. In recent years, the development of strip filling technology provides an economic and effective solution. As a partial filling mining method, strip filling mining can effectively control the movement of overlying strata under the premise of ensuring the strength of filling body, so as to control the surface subsidence and reduce the damage of mining subsidence to the ecological environment. Paste filling material is a kind of concrete filling material [15, 23], which can effectively resist the erosion of underground water and has no hidden danger of spontaneous combustion. Therefore, it will not be destroyed and destabilized due to the influence of other geological factors and can maintain stability for a long time, thus eliminating the potential safety hazard of large-area collapse and even mine earthquake caused by the instability and damage of left coal pillars in strip mining. Compared with full filling mining, strip filling has the advantages of low filling rate, low filling cost, and high mining efficiency.
In this paper, FLAC3D numerical simulation is used to study the sensitivity of the main controlling factors of surface movement and deformation in strip filling mining. The research results will provide basic data for the production practice of “three under” pressure coal mining area, improve the recovery rate of coal resources, reduce the filling cost, and maximize the economic benefits. It has important practical significance to protect ground structures, reduce ecological damage, and realize green sustainable development of coal mine.
2. Analysis of Main Controlling Factors of Surface Movement and Deformation
In the coal measure strata under the thick loose layer, there are usually thick and relatively hard rock layers, which play a major role in controlling the movement of the overlying strata. It is called the key stratum in the strata movement. When the caving method is used to mine the coal seam, after the coal body is excavated, the roof strata will bend and deform towards the goaf and slide down along the bedding plane under the pin joint action of self-weight stress and overlying strata pressure. When the goaf reaches a certain extent, the direct roof will be fractured, and the caving rock mass will roll to the lower part of the goaf. At this time, the basic roof will lose the direct roof support and gradually bend and deform along the normal direction of the layer, and finally, fracture failure will occur. With the advancement of the working face, this process will continue to repeat until the exposed space is completely filled with waste rock. The overburden movement mechanism of strip filling mining is similar to that of strip caving mining under a thick loose layer. Considering the difference of construction technology and the mining method of the strip caving method, it is necessary to consider the displacement of roof and floor after the coal seam is mined out and before the strength of filling body is stable. Because the overburden and thick loose layer of coal measure strata in strip filling mining are also controlled by the key layer, the filling body and bearing rock layer under the key layer present compression and bending deformation mechanism due to the additional vertical stress. The subsidence of the key layer and the strata below it transfers from bottom to top and reaches the maximum when it reaches the upper interface of the key layer. The surface subsidence of strip filling mining under the whole thick loose layer is the internal subsidence of overburden and the movement deformation of the thickened loose layer.
The main controlling factors of surface movement and deformation of strip filling mining under a strip thick loose layer are as follows: (1)Filling rate
The higher the filling rate is, the smaller the moving space of surrounding rock in goaf is, and the less the displacement of roof and floor is. Therefore, the filling rate should be increased as much as possible to control the movement of surrounding rock. (2)Thickness of loose layer
With the increase of loose layer thickness, the depth of coal seam is increasing, and the in situ stress is also rising. After coal mining, the original stress field is destroyed. The greater the original rock stress is, the greater the energy needed to release to reach the new stress balance state. Therefore, the movement and deformation of surrounding rock become more intense. (3)Thickness of key layer
Under the condition of thick and loose, when the filling method is used to treat goaf, the key layer controls the deformation ability of the upper strata. The larger the thickness of the key layer is, the lower the damage degree is. At the same time, the stronger the control ability of overlying strata is, the smaller the surface subsidence coefficient will be. (4)Strength of filling body
After filling into the goaf, the filling body can play a supporting role on the roof and prevent and restrict the further movement and deformation of the surrounding rock to the goaf in time and space, so that only fracture zone and bending subsidence zone are formed in the upper part of the goaf, so as to reduce the separation between magmatic rock and lower strata and avoid the fracture movement of overlying strata. The greater the strength of filling body is, the greater the effect of bearing overlying strata is, the more obvious the controlling effect of strata is, and the smaller the surface movement deformation is.
Through the analysis of the factors affecting the surface deformation of strip filling mining under a thick loose layer, the main control factors in numerical simulation are filling rate, key layer thickness, filling body strength, and loose layer thickness.
3. Numerical Simulation Analysis of the Factors Affecting the Subsidence of Three Strip Filling Mining
3.1. Selection of the Numerical Model
According to the measured geological conditions in the fifth mining area of a mine, the model parameters are determined, and the basic numerical simulation model is established, as shown in Figure 1. The main physical parameters of each layer used in this simulation are shown in Table 1.

In order to simplify the numerical model, the length of the model is 700 m and the width is 220 m. The height of the model is 254 m-334 m due to the different combination structure of overlying strata. The upper boundary of the model is the surface, and the lower boundary is the coal seam floor. In order to eliminate the influence of boundary effect, a 200 m boundary is designed outside the mining boundary along the strike direction, and a 60 m boundary condition is set in the inclined direction, that is, the width of the filling working face is 100 m and the mining thickness is 3.0 m.
3.2. Simulation and Monitoring Scheme
According to the purpose of the study, this paper designs 30 simulation schemes, 20 of which are simulated by a single variable method, that is, when one factor is studied, the other three factors are unchanged, and the simulated study is the surface subsidence law of strip filling mining under different loose layer thickness (160 m, 180 m, 200 m, 220 m, and 240 m), filling rate (51%, 55.5%, 60%, 64.5%, and 69%), filling body strength (2.0 MPa, 2.5 MPa, 3.0 MPa, 3.5 MPa, and 4.0 MPa), and different key layer thickness (40 m, 55 m, 70 m, 85 m, and 100 m). The basic simulation scheme is as follows: filling rate 60%, filling body strength 2.0 MPa, key layer thickness 70 m, and loose layer thickness 200 m. In addition, in order to study the influence of different filling rate and loose layer thickness combination on surface subsidence, the filling rate is set at 51% and 40.5%, and 10 comparison schemes for surface subsidence under different filling rate and loose layer combination are designed. The specific simulation schemes are shown in Table 2.
The filling strip is arranged along the advancing direction of the working face, with 60 m as a filling unit. The corresponding filling and retention ratio of different filling rates are 40.5%-27 : 33 (filling 27 m, leaving 33 m), 51%-34 : 26, 55.5%-37 : 23, 60%-40 : 20, 64.5%-43 : 17, and 69%-46 : 14; the calculation of the filling rate is based on the displacement of roof and floor and construction technology after coal mining, and the calculation formula is as follows:
; the actual filling height is 2.7 m, and the mining width is filling units; is an integer.
At the center of the top (surface) of the model perpendicular to the advancing direction of the working face, a monitoring point is set every 10 m. A total of 70 monitoring points are designed to monitor the horizontal displacement () of surface subsidence during mining.
4. Analysis of the Influence of Filling Rate, Key Layer Thickness, Filling Body Strength, and Loose Layer Thickness on Surface Movement and Deformation
4.1. Influence of Filling Rate on Horizontal Deformation of Ground Surface
For strip filling mining schemes with different filling rates, the variation trend of surface horizontal deformation curve and maximum value of horizontal deformation of the main section of subsidence basin are shown in Figures 2 and 3.


As shown in Figure 2, the surface horizontal deformation curve of the main section of the subsidence basin in each scheme is very similar to that of the curvature curve. It is also about the symmetrical distribution of the center point of the subsidence basin, and the middle part is roughly shaped like a bowl with a convex bottom, similar to the effect formed by the superposition of two sine function curves. When and , the horizontal deformation value is positive. The positive value of horizontal deformation appears near and . The positive value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. When , the horizontal deformation value is negative, and the negative value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. The horizontal deformation value of the interval increases first and then decreases, and the horizontal deformation value of the middle part fluctuates slightly. The overall horizontal deformation value of the interval is greater than that of the two ends of the model.
It can be seen from Figure 3 that with the increase of filling rate, the maximum value of surface horizontal deformation of different mining schemes decreases in stages, and the reduction range first decreases and then increases.
4.2. Influence of Loose Layer Thickness on Horizontal Deformation of Ground Surface
For strip filling mining scheme with different thickness of a loose layer, the variation trend of surface horizontal deformation curve and maximum value of horizontal deformation of the main section of subsidence basin are shown in Figures 4 and 5.


As can be seen from Figure 4, the surface horizontal deformation curve of the main section of the subsidence basin in each scheme has a high similarity with the change shape of the curvature curve, which is also about the symmetrical distribution of the center point of the subsidence basin, and the middle part is similar to the superposition of two sine function curves. When m and m, the horizontal deformation value is positive. The positive value of horizontal deformation appears near m and m. The positive value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. When , the horizontal deformation value is negative, and the negative value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. The negative value of horizontal deformation is greater than the absolute value of the positive value, and the horizontal deformation value of the interval increases first and then decreases, and the horizontal deformation value of the upper part of the goaf changes a little and fluctuates slightly.
It can be seen from Figure 5 that the change trend is basically consistent with that of the maximum curvature.
4.3. Influence of Backfill Strength on Horizontal Deformation of Ground Surface
For strip filling mining scheme with different filling body strength, the variation trend of surface horizontal deformation curve and maximum value of horizontal deformation of the main section of the subsidence basin are shown in Figures 6 and 7.


As can be seen from Figure 6, the surface curvature curve of the main section of the subsidence basin in each scheme is approximately symmetrical with respect to the center point of the subsidence basin. When the strength of filling body is 2.0 MPa, the upper part of the goaf is approximately a bowl-shaped convex at the bottom, which is similar to the effect formed by the superposition of two sine function curves; when the strength of filling body is between 2.5 MPa and 4.0 MPa, the curvature curve is approximately irregular sine function distribution, and the greater the strength of filling body is, the smoother the curve is. From the boundary of the goaf to the edge of the subsidence basin, the convexity on the curve is larger, and the change range near the boundary of the goaf is larger than that near the edge of the subsidence basin. When m and , the horizontal deformation value is positive. The positive value of horizontal deformation appears near m and m. The positive value of horizontal deformation is, respectively, mm/m, mm/m, mm/m, mm/m, and mm/m. When , the horizontal deformation value is negative, and the negative value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. The main reason is that the deformation value of the mined-out area tends to be smaller in the upper part of the mining area than that in the negative pole area.
It can be seen from Figure 7 that with the increase of filling body strength, the maximum value of horizontal surface deformation of different mining schemes decreases linearly and fluctuates when the strength of filling body is 3.0 MPa. The maximum value of horizontal deformation increases with the decrease of filling body strength. When the strength of filling body exceeds 3.5 MPa, the decrease range of maximum horizontal deformation changes from large to small.
4.4. Influence of Key Layer Thickness on Horizontal Surface Deformation
For strip filling mining scheme with different thickness of key stratum, the variation trend of surface horizontal deformation curve and maximum value of horizontal deformation of the main section of the subsidence basin are shown in Figures 8 and 9.


As can be seen from Figure 8, the surface curvature curve of the main section of the subsidence basin in each scheme is approximately symmetrical with respect to the center point of the subsidence basin. The scheme with the thickness of the key layer of 70 m, 85 m, and 100 m is roughly irregular sine function distribution, the upper part of the goaf is roughly irregular bowl shape, while the scheme with the key layer thickness of 40 m and 55 m has two sinusoids above the goaf. The wave effect is formed by superposition of function curves. When m and , the horizontal deformation value is positive. The positive value of horizontal deformation appears near m and m. The positive and negative extrema of horizontal deformation are mm/m, mm/m, mm/m, mm/m, and mm/m. When , the horizontal deformation value is negative, and the negative value of horizontal deformation is mm/m, mm/m, mm/m, mm/m, and mm/m. The negative value of horizontal deformation is greater than the absolute value of the positive value, and the horizontal deformation value of the interval increases first and then decreases, and the horizontal deformation value of the middle part changes a little and fluctuates slightly.
It can be seen from Figure 9 that when the thickness of the key layer is within a certain range and the stability of the filling body is ensured, the maximum surface curvature decreases approximately linearly with the increase of the thickness of the key layer. When the thickness of the key layer changes from 85 m to 100 m, the maximum absolute value of surface tilt changes from decreasing trend to increasing trend.
5. Sensitivity Analysis of Filling Rate, Key Layer Thickness, Filling Body Strength, and Loose Layer Thickness to Surface Movement and Deformation
5.1. Orthogonal Test Analysis
From the point of view of mathematical statistics, according to the principle of orthogonality, the orthogonal test method is a method to select high quality, representative, and typical test points from a large number of complex test points and then arrange the experiment with “orthogonal table.” The selected points in the experiment have the characteristics of “uniform dispersion, uniformity, and comparability.” The principle of orthogonal experiment is as follows: first, before the orthogonal test, with the help of the prepared orthogonal table, the experimental scheme is arranged in a planned and purposeful way, and then, the experiment is carried out; secondly, after the uncomplicated table operation, the test results are correctly analyzed based on the calculation value. The advantages of orthogonal test are as follows: on the one hand, it can distinguish the primary and secondary effects of each factor in the experiment through less test times; on the other hand, it can distinguish the size of the effect of each factor on the index.
5.2. Overview of Range Analysis
According to the orthogonal table experiment, we can get the test results of a certain index (single index or multi-index). Through “range analysis” or “variance analysis,” we can get the best experimental scheme and the order of primary and secondary factors.
The steps of “range analysis” test results (take four factors and three levels as an example, see Table 3) are as follows: (1)The results of the first level of each factor in each test were summed up, namely, ; then, the two-level results of each factor in each test, namely, ; and the three-level results of each factor in each test, namely, (2)The average values of the test results of each factor and level were calculated, respectively, , , and , and filled in the orthogonal table(3)The difference (also known as range) of the average value of each level of each factor is calculated, and the difference between the maximum value and the minimum value of the average value of each level is calculated
5.3. Orthogonal Experimental Design and Analysis of Experimental Results
The purpose of this paper is to determine the sensitivity of filling rate, filling body strength, key layer thickness, and loose layer thickness to surface deformation. Due to the large number of numerical simulation models, influencing factors, and the influence levels of each factor, it is more complicated to calculate one by one. By using the orthogonal experimental design method and orthogonal table to determine the experimental simulation scheme, the workload of calculation can be effectively reduced. According to the principle of orthogonal experiment design, for the above four factors, each factor takes three representative levels, so L9(34) orthogonal table is selected to arrange the experiment, and the factor levels are shown in Table 4.
5.4. Orthogonal Test on Factors Influencing Maximum Horizontal Deformation
Nine numerical simulation schemes are designed for the orthogonal test of the factors affecting the maximum curvature. The schemes and results of the orthogonal test are shown in Table 5, and the range analysis results are shown in Table 6.
As can be seen from Table 6, the range of filling rate is the largest, which is 0.38; the range of key layer thickness and filling rate is 0.20 and 0.17, respectively; and the range of loose layer thickness is the smallest, which is 0.06. The range analysis results of the orthogonal test show that the filling rate has the greatest influence on the surface subsidence control effect, which is the main control factor of strip filling mining control surface subsidence; the key layer thickness and filling body strength are ranked second, and the loose layer thickness is the smallest relative to the other three factors.
6. Conclusion
In summary, this paper draws the following conclusions: (1)The sensitivity of filling rate, key layer thickness, filling body strength, and loose layer thickness to surface movement and deformation is in the order of filling rate > key layer thickness > filling body strength > loose layer thickness(2)With the increase of filling rate, the horizontal deformation of the main section decreases gradually; when the filling rate increases from 64.5% to 69%, the reduction range of the maximum value of the surface horizontal deformation first decreases and then increases(3)With the increase of loose layer thickness, the lower horizontal deformation value of surface movement deformation first decreases and then increases. When the loose layer thickness is in the range of 160 m-180 m, the maximum value of surface movement deformation parameters gradually decreases with the increase of loose layer thickness. When the thickness of a loose layer is in the range of 180 m-240 m, the maximum value of surface movement deformation parameters gradually increases, and the increase range is small in the range of 200 m-220 m(4)With the increase of backfill strength, the maximum value of horizontal deformation decreases linearly with the increase of backfill strength and fluctuates when the strength of backfill body is 3.0 MPa. The maximum value of horizontal deformation increases with the decrease of backfill strength. When the strength of backfill body exceeds 3.5 MPa, the decrease range of maximum horizontal deformation decreases from large to small
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest related to the publication of this paper.
Acknowledgments
This research was funded by the Postdoctoral Innovation Project of Shandong Province (Grant No. 202003090), the Postdoctoral Science Foundation of China (Grant No. 2020M682267), and the Engineering Laboratory of Deep Mine Rockburst Disaster Assessment Open Project (LMYK2020002).