Overburden and surface subsidence with slicing paste filling mining in thick coal seams

To overcome the difficulties of overburden failure and surface subsidence induced by the slicing mining of extra-thick coal seams, slicing filler paste is typically utilized. In this paper, a continuous curved beam mechanical model of paste filling mining was established by theoretical analysis against the background of the 3305 working face of Yangcheng Coal Mine, the main controlling factors of surface subsidence were analyzed using an orthogonal experiment method (OEM) and an analytic network process (ANP) coupled comprehensive assignment, and the vertical displacement distribution of the overlying rock under different filling intensity was simulated using numerical simulation software. The following main findings were obtained. First, the elastic modulus of the backfill is the primary regulating factor, as demonstrated by the results. Second, the greater the filling body’s elastic modulus, the more efficiently it carries the overburden load and minimizes the maximum surface subsidence. Third, the distance from the open-cut mine likewise reflects the degree of surface subsidence, with the amount of subsidence increasing as the distance from the mine increases.


Introduction
In recent years, coal resources have continued to dominate China's energy structure, comprising more than 50% of the total share of all energy sources. However, the environmental issues caused by roof caving incidents and surface subsidence as a result of disturbances in coal mining remain severe. According to the relevant statistics (Wu and Li, 1995), it is typical in East China to mine coal under buildings, railways, and water bodies, which not only wastes resources but also compromises the integrity of the earth if the mining method is inappropriate. The backfill mining technique, which may prevent surface subsidence and maximize coal resources, is currently receiving widespread attention from domestic experts.
Numerous professionals and academics are now conducting research on roof overburden damage and surface subsidence, mostly through indoor experiments, theoretical analysis, numerical simulation, and field observations. Deng et al. (2015) suggested a novel mining technique of upward slicing filling mining for the geological circumstances of extra-thick coal seams using indoor experiments, investigated and researched this filling effect using similar material modeling, and conducted engineering practice experiments. Yang et al. (2015) utilized similarity theory to create comparable physical models to simulate the mining operation of underground ore deposits and monitor displacement changes. Liu et al. (2017a) constructed the theory of continuous curving beams for backfill mining using a large number of engineering applications and theoretical analysis as a starting point, revealing the mechanism for backfill control on the roof. Gao et al. (2017) established a defect-coupled constitutive model based on the theory of damage mechanics and showed that the piecewise nonlinear model more accurately depicts the effect of complex backfill defects on the stress-strain curve. Using geological data from Xiaotun Mine and numerical modeling, Zhang et al. (2011) assessed the influence of the filling ratio on surface subsidence. The results indicated that when the filling ratio increases, the control impact on surface subsidence improves. Kostecki and Spearing (2015) utilized finite-difference software to model the rheological characteristics of coal pillars in backfill mining and studied the link between the mechanical properties of the backfill material and the coal pillars under backfill mining settings. Numerous researchers have undertaken studies on the techniques and methods of backfilling, as well as the materials used in conjunction with the site, to avoid and manage roofcaused surface subsidence problems. Zhao et al. (2019) established an enhanced analytical method for evaluating the safety of vertically exposed backfill faces and maintaining the stability of the exposed backfill face. Thirukumaran et al. (2016) stated that the influence of the filling body on the surrounding rock is dependent on the mechanism, which consists of local support, surface support, and overall support. Gong et al. (2017) and Guo et al. (2018) investigated the physical parameters of filling with coal gangue as aggregate and the rationale and applicability of coal gangue mixture as backfill material.
During thick coal seam paste filling mining, the filling body in the mined-out region may handle the majority of the overburden load, which is one of the characteristics that separates it from caving. There is currently no consensus regarding important theories and application solutions for the site's study status of the stress development law of the surrounding rock, overburdened rock, and surface movement deformation characteristics. The correlation between the strength of the filling material and the quantity of surface subsidence requires more investigation. To solve the problems of roof fall accidents and surface subsidence caused by thick coal seam mining, this paper analyzes the influencing factors of backfill, the effect of roof control, and predicts the amount of surface subsidence. Combining statistical principles and evaluation methods, the mining plan was finally chosen as using cement as the cementitious material, using tailing sand as the aggregate, mixed with fly ash and industrial slag as the filling material for backfill mining with the research background of the 3305 working face of Yangcheng Coal Mine.

Summary of the study area
Yangcheng Coal Mine is located in Jining City, Shandong Province, in the Liangshan coalfield in the southwestern block of the Luxi fault block and belongs to the North China sedimentary coalfield. Its stratigraphic distribution is Ordovician, Carboniferous, Permian, and Quaternary in order from the bottom to the top. There are two major coalbearing strata in the study area, namely, the Upper Carboniferous Taiyuan Formation (C3t) and the Lower Permian Shanxi Formation (P1s). The 3# coal seams are the recoverable coal seams in the whole area, with the thickness ranging from 4 m to 9.5 m, and the average thickness being 7.5 m. The mine has an annual production capacity of 2.1 million t/a and a recoverable coal reserve of 133 million t in a near-horizontal seam. The 3305 working face is located on the west side of the third mining area, and the upper part is the 3303 working face goaf. The average coal seam thickness of the working face is 7.2 m, and the seam floor elevation is −870 m to −990 m. Due to the damage caused to the buildings and farmland of the nearby villages above the ground surface during the fully mechanized caving mining of the working face, the surface subsidence caused by mining was mitigated by using the paste filling mining method when extracting the coal seam of working face 3305; because the coal seam was thick, the whole coal seam could not be filled at one time, so the thick coal seam was extracted by using the slicing paste filling method. The working face layout and filling site are shown in Figure 1.

Mechanical modeling of roof overburden
Analysis of the mechanical model of continuous curved beams During the mining operation, as the slurry is continuously filled into the goaf, the backfill material gradually solidifies and forms a stable support carrier, thereby limiting the sinking of the overlying rock layer, causing the rock layer to bend and sink, and preventing the irregular collapse of the roof (Shi et al., 2021). The pressure on the roof above the thick coal seam can be considered the uniform load (q 0 ), while the roof's rocking beam is considered the elastic foundation beam. In addition, we consider the subterranean ceiling, the filling body, and the coal seam floor to be the elastic base. Figure 2   depicts a structural mechanical model of thick coal slicing paste filling mining based on the aforementioned assumptions.
Due to the supporting effect of the backfill body, the roof is more stable and less prone to fracture during the mining of thick coal seams with slicing paste filling. After mining damage, the roof is no longer supported on four sides and resembles a "simply supported beam" (Sandholtz et al., 2020). This structure's mechanical model matches that of the upper and lower slicing mining surface. To assist further qualitative calculations, the load q is approximated as a uniform load q 0 , and a schematic diagram of the force is presented in Figure 3.
According to the Winkler foundation hypothesis, the pressure on any point on the foundation is positively correlated with the settlement at that point (Jena et al., 2020). Take the upper slicing face as an example, then, where p 1 (x) is the support of the main roof by the foundation in the backfill area; p 2 (x) is the support of the main roof by the foundation in the coal seam area; w 1 (x) is the deflection of the main roof rock beam in the backfill area; w 2 (x) is the deflection of the main roof rock beam in the coal seam area; k 1 , k 2 is the elastic foundation coefficient determined by the elastic modulus; L is the strike length of the backfill area; and q 0 is the uniform load.
When mining the thick coal seam with slicing and filling, for the main roof deflection differential equation:   Frontiers in Earth Science frontiersin.org 03 where E is the roof elastic modulus; I is the moment of inertia of the roof section; and u 1 is the final subsidence amount of the immediate roof contacting the backfill body. Set and substitute the last formula: By solving the aforementioned equation we get Since the mining method and the backfill body of the upper and the lower slicing face are the same, the amount of subsidence in the aforementioned Eq. 4 is the same. Combined with the boundary conditions, substituting x=0, x ± π/2 into Eq. 4, it can be determined that A, B, C, and D are as shown in Eq. 5: Finally, the deflection curve equations of the upper slicing face and coal seam area are obtained as follows: In the same way, the deflection curve equation of the lower slicing face and the coal seam area can be obtained as follows (Jena et al., 2020): Main roof load and foundation coefficient

Main roof load
The main roof load during backfill mining can be regarded as the sum of the total pressure of the overlying strata, which is as follows: Substituting the corresponding parameters of the overlying strata on the 3305 working face of Yangcheng Coal Mine 3# coal paste filling mining into the aforementioned formula, the main roof load q can be obtained as 22.5 MPa.

Foundation coefficient
Through the aforementioned analysis, it can be seen that the force and deformation of the roof are closely related to the foundation coefficient k (Chen and Park, 2020). When the coal seam is filled and mined, the elastic foundation is composed of the roof and the backfill body, which is regarded as an elastic body as a whole, and the constitutive relationship can be obtained as follows: that is, the total compression of the immediate roof and the backfill body is as follows: where: E 1 is the immediate roof elastic modulus; E 2 is the backfill body elastic modulus (Liu et al., 2017b); y 1 is the immediate roof compression; y 2 is the backfill body compression; h 1 is the immediate roof height (Mangal, 2021); and h 2 is the backfill body height.
According to the Winkler foundation assumption: where σ is the overburden load; k is the foundation coefficient; and y is the total compression of the roof and backfill body.
The foundation coefficient of the 3# coal seam lower slicing filling mining is similar to that of the upper slicing filling mining. Combining Eqs. 9-11, there is the following relationship between the foundation coefficient of the 0≤x≤L section and the foundation coefficient of the x≤0 section in the upper and lower slicing backfill mining: Frontiers in Earth Science frontiersin.org where k 1 is the foundation coefficient of upper slicing filling mining on 3# coal in the 0≤x≤L interval; k 2 is the foundation coefficient of upper slicing filling mining on 3# coal in the x≤0 interval; k 3 is the foundation coefficient of lower slicing filling mining on 3# coal in the 0≤x≤L interval; k 4 is the foundation coefficient of lower slicing filling mining on 3# coal in x≤0 interval; E 1 is the immediate roof elastic modulus; E 2 is the backfill body elastic modulus; E 3 is the coal seam elastic modulus; E 4 is the main roof elastic modulus; h 1 is the immediate roof height; h 2 is the height of the upper slicing backfill body; h 3 is the coal seam thickness; and h 4 is the backfill body height.
The relevant parameters of the upper face are as follows: . Substituting into the relevant parameters, we can obtain the foundation coefficient k 1 =0.14 GN/m 3 , k 2 =0.54 GN/m 3 .
The parameters of the immediate roof and the backfill body during the lower slicing mining are the same as those of the upper slicing mining. The overall height h 4 of the 3# coal backfill body is 6.48 m. Substituting the relevant parameters can obtain the foundation coefficient k 3 =0.073 GN/m 3 , k 4 =0.33 GN/m 3 .

Mechanical calculation of roof subsidence
Regarding the main roof subsidence characteristics of foundation in the filling area, from the aforementioned analysis, it can be seen that the main roof subsidence of the filling area of 3# coal upper and lower slicing paste filling mining is expressed as follows: The buried depth of the 3# coal seam is 900 m. The upper and the lower slicing mining areas have the same roof subsidence during the slicing paste filling mining; both are 0.36 m. We find the foundation coefficient k 1 =0.14 GN/m 3 , k 2 =0.57 GN/m 3 in the upper slicing filling mining, and the foundation coefficient k 3 =0.073 GN/m 3 , k 4 =0.33 GN/m 3 , with the main roof load q=22.5 MPa in the lower slicing filling mining. The mechanical parameters of each influencing factor are shown in Table 1. We then substitute the relevant parameters in Table 1 into Eq. 13 and analyze the main roof subsidence's characteristics of the upper and lower slicing mining filling areas (using MATLAB) and draw the fitting curve (using Origin) (Heinze et al., 2021). The schematic diagram of the influence curve of various factors on the main roof subsidence during upper slicing mining is shown in Figure 4, and the schematic diagram of the influence curve of various factors on the main roof subsidence during lower slicing mining is shown in Figure 5 ( Deng et al., 2021).
In the upper and lower slicing mining, the influence of many sources on the principal roof subsidence has a clear, consistent tendency. As the distance from the coal wall (i.e., the distance of the open-off cut from the coal wall) grows, the backfill pressure first increases dramatically, peaks after a given distance, and then decreases gradually until it reaches the original rock stress (Yang et al., 2014). The following are the mechanical properties of the influencing factors' parameters: 1) The modulus of elasticity and the thickness of the main roof has a certain influence on the pressure of the backfill body during coal seam filling and mining, which is manifested as follows: as the elasticity modulus and the thickness of the main roof increase, the maximum value of the pressure on the backfill body is increasingly further from the wall. Compared with the main roof elastic modulus, the influence of the main roof thickness is relatively significant.
2) The buried depth of the coal seam has a more obvious impact on the pressure of the backfill. The further the buried depth of the coal seam is from the surface, the greater the pressure of the backfill and the higher its peak value. 3) During coal seam filling and mining, the amount of subsidence when the roof touches the backfill also has a significant influence on the pressure on it. The greater the subsidence value, the higher the peak pressure of the backfill body. 4) During coal seam filling and mining, regarding the backfill body, the larger the elastic modulus, the larger the peak pressure on it, and as the peak pressure point gradually approaches the wall, the smaller the range of pressure on the backfill body. 5) Regarding the elastic modulus of the coal, as the elastic modulus of the coal gradually increases, the peak pressure of the backfill gradually decreases, but the pressure on the backfill is less affected.
Comprehensive risk assessment and analysis of main roof subsidence based on the adopted method The alternative ranges of the main control factors are selected by the foundation coefficients and hydrogeological data, and the ranking of the importance of the weights is determined by using the comprehensive assignment method coupled with the analytic network process (ANP) and orthogonal experiment method (OEM) (Mosaad and Basheer, 2020). With the assessment methods based on expert experience and documentary data, the advantages of subjective analysis and objective analysis can be combined to improve prediction accuracy and credibility.
Frontiers in Earth Science frontiersin.org 05      Analytic network process (ANP) An overview of the analytic network process (ANP) Thomas L. Saaty enhanced the analytic network technique based on the analytic hierarchy process (AHP) (Huang et al., 2020), which is an evaluation approach that combines qualitative and quantitative analysis. ANP can take into account the internal connections between nodes of various factor groups in a more thorough way than the classic AHP, which merely stresses the flaws of reciprocal impact across criteria levels (Kundu et al., 2021).
The weights are calculated as follows: 1) Construction of the judgment matrix. The control layer of the ANP has criterion elements B 1 , B 2 , B 3 ..., and the network layer under the control layer has element groups B i1 , B i2 , B i3 ..., where i=1,2,....,n (Mirarabrazi and Navrodi, 2020). Take the control layer elements as the criteria and the network layer elements as the secondary criteria. The elements in the element group are compared according to their influence on B ij , and the construction judgment matrix is obtained:

FIGURE 7
Distribution diagram of the weights for the main roof subsidence factors.

FIGURE 6
Three-cluster ANP network model.
Frontiers in Earth Science frontiersin.org    2) The ANP super cluster-weighted matrix and limit matrix (Camara et al., 2020). Combining the ranking vectors of the degree of mutual influence of the elements of the network layer, the supermatrix under the control element is obtained: Each element of the aforementioned matrix is a non-negative matrix, and the columns of each supermatrix's submatrices are normalized, but the entire matrix's columns are not. The relevance of each element group is compared using the control level B as the criterion and the arbitrary element group B i as the sub-criterion to generate the characteristic matrix: Use the formula W W · P to obtain the weighted supermatrix and use the formula W ∞ lim n → ∞ W n to obtain the limit matrix so as to obtain the importance ranking of each index.

Weight calculation by the analytic network process (ANP) model
According to the analysis of the mechanical model, we obtain six factors that affect the main roof subsidence: main roof elastic modulus; main roof thickness; buried depth of the coal seam; coal elastic modulus; backfill elastic modulus; and the subsidence of the roof contacting the backfill. Based on the attributes of the factors, they can be classified into three categorizations: roof structure; backfill characteristics; and coal lithology. A threecluster ANP network model was developed as shown in Figure 6.
Using group decision-making expert scoring, using 1-9 scales as the measurement criteria, and using (0,1) as the scoring range for the importance of the criteria or elements at the same level, pairwise comparison scores are obtained. The supermatrix composed of the judgment matrix is shown in Table 2, Table 3, Table 4, and Table 5.
The maximum eigenroots and the corresponding eigenvectors are calculated according to the limit matrix, and then the eigenvectors are normalized to obtain the final weight of the roof subsidence amount evaluation index, as shown in Figure 7.
It can be seen from Figure 7 that the backfill elastic modulus and the subsidence of the roof contacting the backfill are the main controlling factors. They play a decisive role in the failure state of the overlying rock during sliced filling mining.

Verification experiment
Overview of the orthogonal experiment method (OEM) Orthogonal experimental design and analysis methods are scientific calculation methods based on mathematical statistics, probability theory, and practical experience, using standardized orthogonal table design schemes and analyzing the results of precision and sensitivity. In the analysis of actual problems, the use of orthogonal tables to optimize the test plan can reduce the number of experiments while ensuring a certain level of accuracy, efficiently handling multi-factor optimization problems, and finding the main factors that have a significant impact on the test indicators among many factors. The design process of the OEM (Xuan and Leung, 2011) is as follows: 1) determine the experimental factors and the number of levels; 2) select the applicable orthogonal table; 3) list the experimental plan and the experimental results; 4) perform range analysis and variance analysis on the experimental results; and 5) determine the optimal or better combination of factor levels.  Table 6.
We use the data in Table 6 to obtain the k i value and the range value R of each factor through orthogonal calculation, and the results are shown in Table 7.
The importance ranking was performed by finding the comprehensive average of the five levels (k 1 , k 2 , k 3 , k 4 , and k 5 ) of the six main control factors and conducting the range analysis to find the range value R corresponding to each factor. The order of importance positively reflects the importance of the factors, i.e., the first one is the most important and has the greatest weight, and so on. As can be seen from Table 7, the first two are the backfill elastic modulus and the subsidence of the roof contacting the backfill.

Determine the main control factors
The weight of contribution to the aim, as determined by the ANP expert scoring panel, is as follows: backfill elastic modulus > roof subsidence contacting the backfill > coal elastic modulus > buried depth of the coal seam > main roof thickness > main roof elastic modulus. The results were confirmed by the construction of orthogonal experiments based on objective data, and the ranking of the relevance of the main control factors derived by analysis of variance was compatible with the empirical analysis-based ANP group decision scoring approach. It can be determined that the elastic modulus of the backfill body is the main control factor for the subsidence of the backfill mining, and simulations were conducted for various orders of the backfill elastic modulus. In addition, the filling material that can reach the ideal modulus was chosen to handle the surface subsidence and destruction problem.
Numerical simulation of movement and deformation of overlying strata in the slicing paste filling mining of thick coal seams Geomechanical parameters of the model The numerical calculation model is established based on the engineering background of the geological mining conditions of the 3305 filling face of Yangcheng Coal Mine (Miranda et al., 2011). The simulated coal seam thickness is 7.2 m, the coal seam advancement distance is 750 m, and the working face width is 120 m. Due to the deep location of the coal seam, only the positions between 120 m above and below the coal seam (distance from the z-direction 240 m) are simulated, and the overburden gravity of 21.25 MPa is applied on the top. Considering the boundary effect factor, the simulated size is determined to be 1,500 m*350 m*463 m. The working face adopts the operation of aligning itself with the inclination, mining toward the direction, and filling while mining. The advancing amount is 1,000 m with 250-m coal pillars reserved at both ends, and 50-m protection coal pillars are reserved at both ends of the inclination. The simulation excavation is set with 20 excavation steps, where each step is 50 m, and a total of 109,224 grids are generated by the finite difference method of FLAC3D. From the 60 observation points, we selected five representative points (A, B, C, D, and E) with different orientations as sample data collection points. Taking the orientation of the model as the X axis, the inclination as the Y axis, and the height as the Z axis, an XOY space rectangular coordinate system is established. The coordinates of A, B, C, D, and E are as follows: A (300; 100;400); B (500;150;360); C (700;200;320); D (900;280;280); and E (1,100;300;240).The geomechanical parameters of coal and rock are shown in Table 8, and the numerical model is shown in Figure 8.
Overlying strata and surface deformation characteristics of thick coal seam slicing paste filling mining under different backfill strengths Paste filling materials can be calculated by the Salamon formula double-yield plasticity model Wang and Du, 2020;Yang et al., 2020). From the conclusion of the previous section, it can be seen that the strength of the backfill is a key factor affecting the Frontiers in Earth Science frontiersin.org filling effect of the paste (Liu et al., 2017c). In order to verify the effect of backfill strength on overburden and surface control, numerical simulation was conducted to analyze the overburden and surface subsidence by comparing three cases of 0.3 GPa, 0.5 GPa, and 1.0 GPa backfill strength at a 95% fill rate. We compared three different backfill strengths of 0.3 GPa, 0.5 GPa, and 1.0 GPa. We took the vertical displacement distribution of overburden strata after all the advancing of the 3# coal 3305 working face as an example for comparison, as shown in Figure 9, and the surface subsidence comparison graph is shown in Figure 10.
Five points with different orientations (A, B, C, D, and E) were selected from the 60 observation points in FLAC3D numerical simulation to observe the amount of roof subsidence corresponding to the filling of the slicing mining with different elastic moduli of 0.3 GPa, 0.5 GPa, and 1.0 GPa at a 95% filling rate. With the advancement of the 3305 working face, the overlying rock layer gradually sinks on the original foundation, and the amount of subsidence increases with the advancement of mining. The upper slicing face reaches its peak at an excavation step distance of about 800 m, and the lower slicing face reaches its peak at an excavation step distance of about 1,500 m; until the direction of mining is infinite, the amount of subsidence tends to be stable. We observed the change of the effect of surface subsidence and overburden damage at five points (A, B, C, D, and E) using backfill with different elastic moduli. The subsidence of the roof with mining damage was ranked 0.3 GPa>0.5 GPa>1.0 GPa. In the maximum settlement curves of the upper and lower strata, 1.0 GPa<0.5 GPa<0.3 GPa, and the peak of the main roof subsidence of the upper slicing face is greater than that of the lower slicing face.

Analysis of mine pressure monitoring data
According to the mine monitoring pressure data before and after backfill mining at the 3305 working face, a comparative analysis was conducted. With the advancement of mining, the pressure on the protective coal pillars showed different changes, and the mine monitoring pressure graphs before and after filling were obtained as shown in Figure 11.
From Figure 11A, it can be seen that in the unfilled area, the mine pressure is obvious and the pillar is under greater force. During the mining process, the mine pressure started to decrease at 10 m of excavation because the coal body was extracted one after another. When the excavation reaches 40 m, the periodic incoming pressure appears, and comparing this part of mine pressure with the previous numerical simulation, it can be found that it basically matches. According to Figure 11B, it can be seen that, after filling, the mine incoming pressure shows a large and stable decrease, and the increase tends to be stable in the later period. The peak and average values of the pillar pressure are significantly reduced, which reduces the risk of rock impact and roof collapse and deformation. Therefore, it can be concluded that slicing paste filling mining with cement, tailing sand, fly ash, and industrial slag as the filling materials for thick coal seams can effectively support the roof and reduce the amount of subsidence.

Conclusion
1) The mechanical model of the roof elastic foundation beam model was constructed using mechanical knowledge, and the deflection curve equations in the filling area and coal body area were established. Through the force analysis of the roof, the range of factors that affect the amount of roof sinking was obtained. The relationship curves between the roof elastic modulus, the roof thickness, the buried depth of the coal seam, the elastic modulus of backfill, the coal elastic modulus, and the subsidence were drawn. 2) From the expression of the main roof subsidence in the filling area of 3# coal upper and lower slicing paste filling mining and the geological data of the 3305 working face of Yangcheng Coal Mine, six influencing factors affecting the amount of main roof subsidence were inductively analyzed, and the weights were derived using the subjective group decision method of the ANP and verified by the OEM based on objective data, and the importance ranking coinciding with the weight ranking was derived. The results show that the backfill elastic modulus is the main controlling factor. 3) FLAC3D was used to model backfill elastic moduli of 0.3 GPa, 0.5 GPa, and 1 GPa during excavation while maintaining a constant filling rate. The simulation revealed the effect of the filling body's elastic modulus on the overburden and surface sinking. The results demonstrate that increases the elastic modulus of the backfill reduces the degree of overburden damage and the quantity of surface subsidence, indicating that increasing the elastic modulus of the filling body can effectively enhance the effect while maintaining the filling rate.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary material. intellectual contribution to the work and approved it for publication. QZ was responsible for writing the paper, deriving the model, and translating and embellishing it. CW directed the structure of the thesis. LP was responsible for the literature proofreading.