Numerical Study on the Movement Laws of Overlying Strata in Shallow-Buried Stope Based on the Goaf Compaction Effect

&is study presents an integrated approach including the theoretical analysis and numerical modelling to investigate the failure characteristics of the overlying strata in the shallow-buried stope. &e mechanical characteristics of the caving zone are first revealed and then calibrated by using the double-yield model. &e theoretical results show that the mechanical properties of the collapsed rock mass are closely related to its crushing expansion coefficient and uniaxial compressive strength. &e vertical stress of the collapsed rock mass increases slowly with the strain and then increases exponentially after a certain critical strain. &e simulation indicates that the fracture zone volume is 1.7-1.8 times that of the caving zone in the 31108 working face, and the failure volume of the overlying strata is 9-10 times that of the stope. &e simulated height of the caving zone and fracture zone is 9m and 20m, respectively. &e comparison between the numerical and field measurement results demonstrates that the new evaluation method using shear-tensile strain behaviors can accurately predicate the height of the two zones.&e proposed numerical method could be a viable alternative approach to two zones height calculation.


Introduction
Western China is rich in coal resources, accounting for 65% of the country's total reserves, but water resources reserves only account for 3.9% [1,2]. e lack of water resources seriously limits the sustainable development of energy in Western China. Considering the serious contradiction between coal resources development and water resources protection in the western region, Qian et al. [3] presented the concept of green coal mining and water conservation mining technology in the western mining area. Accordingly, underground reservoir technology of coal mine was first proposed by Shenhua Group. In this way, water storage in goaf formed after coal mining has become an important technical way to handle the contradiction between coal mining and water resources protection in Western China. In the process of coal mining, the caving zone, fracture zone, and bending subsidence zone are formed from bottom to top with the movement of the overlying strata [4][5][6]. e water storage capacity of the underground reservoir is directly determined by the failure mode and height of the caving zone and fracture zone [7]. erefore, it is urgent to understand the failure mode and height of the overlying strata during coal seam mining.
To date, many studies have been conducted to investigate the failure law of the overlying strata during underground reservoir construction. For example, Guo et al. [8] investigated the process of overburden failure transmission and further proposed a method for predicting overburden failure height. Xu et al. [9] presented the development law of a water flowing fractured zone with the breaking of the key stratum based on the key stratum theory. Huang et al. [10] proposed a method to characterize the distribution of the fractured zone by using the tensile deformation rate of the overlying strata. Yang et al. [11] performed a field observation to reveal the development law of the water flowing fractured zone of the soft overburden above a thick coal seam by monitoring borehole water injection leakage. Zhang et al. [12] revealed the influence of mining height and the panel layout on the height of the two zones in the soft overlying strata using discrete element numerical simulation. e above research studies provide a better understanding of the movement laws of overlying strata.
Currently, theoretical analysis, field measurement, and numerical modelling are the most used method to predict the development height of the water flowing fractured zone.
ere into, theoretical calculation methods are mostly based on various assumptions with low accuracy, while field measurement improves the accuracy of the measurement results, but the test process is time-consuming and difficult. erefore, numerical simulation has become a simple and effective method for analyzing the failure characteristics of overlying strata [13,14]. In addition, numerical modelling can involve various influence factors in the analysis. Relevant studies have shown that the goaf compaction effect can result in characteristics changes of the overlying strata movement by providing an additional support resistance to the roof strata [15,16]. Nevertheless, existing studies have rarely taken this influence into account. erefore, it is urgent to conduct numerical research on the movement laws of overlying strata considering the compaction effect of goaf.
In this study, a numerical method by using FLAC3D is presented to investigate the range of caving and fracture zones by discerning the plastic strain value of the zones in the model. In this modelling method, the compaction effect of the goaf is achieved by using the double-yield constitution model. e test site is located in the Lijiahao coal mine, Ordos City, Inner Mongolia Autonomous Region, China.
e results presented in this study further improve the accuracy of numerical simulation, which is of great help to calculate the underground reservoir capacity.

Project Overview
e Lijiahao coal mine is located in Dongsheng District, Ordos City, Inner Mongolia Autonomous Region, China, with a designed production capacity of six million tons yearly. Currently, the mining seam is 3-1# coal seam. e average overburden depth and thickness of the coal seam are 201-254 m and 2.5-6.3 m, respectively. e coal seam is a relatively simple structure, and the immediate roof is sandy mudstone, with an average thickness of 3.72 m. e main roof is fine sandstone with an average thickness of 5.00 m, and the floor is mainly sandstone. e test site is located in the 31108 working face. Figure 1 shows the coal seam histogram at this position.

Mechanical Properties of the Collapsed and Compacted
Rock Mass. A large number of field observation data on strata movement has demonstrated that overlying strata can be divided into caving, fracture, and bending subsidence zone from bottom to top after the coal seam is excavated. e collapsed rock mass is transformed from a loose body to a supporting body under the compaction action induced by the gravity and overlying strata movement [17][18][19]. e fractures development and stress distribution of the overlying strata is significantly affected by the strain hardening behavior of the collapsed rock mass. Salamon [20] proposed the stress-strain relationship of rock mass in the caving zone by taking the caving block as a granular material. e expression is shown in the following equation: where σ cap , ε, ε max , and E 0 are the vertical stress of the collapsed rock mass, the strain of collapsed rock mass under the action of vertical stress, the maximum vertical strain, and the initial elastic modulus of the collapsed rock mass, respectively. ε max depends on the crushing expansion coefficient K p of the collapsed rock mass, and the expression is as follows [21]: e relationship between the initial modulus E 0 , crushing expansion coefficient K p , and compressive strength of the collapsed rock mass σ c is expressed as follows [22]: (3) Figure 2 illustrates the effect of the crushing expansion coefficient and uniaxial compressive strength on the initial modulus of a collapsed rock mass. e initial modulus of the collapsed rock mass decreases exponentially with the increase in the crushing expansion coefficient but decreases gently with the uniaxial compressive strength decreasing.   e initial modulus of the collapsed rock mass decreases greatly when the crushing expansion coefficient is less than 1.6; conversely, it decreases gently when the crushing expansion coefficient is larger than 1.6.
From equations (1)-(3), the expression of the mechanical properties of the collapsed rock is deduced, as shown in the following equation.
e relationship of the stress-strain of the collapsed rock mass with different crushing expansion coefficients and the compressive strength is illustrated in Figure 3. With the strain increasing, the vertical stress of the collapsed rock mass increases slowly; and then, it increases rapidly at an exponential function when the strain reaches a certain critical value. When the crushing expansion coefficient is equal, the compressive stiffness of the collapsed rock mass increases with the increase of its uniaxial compressive strength. In other words, the strain value of the low-strength rock mass is larger than that of the high-strength rock mass with the same stress value [23,24].

Determination of Mechanical Parameters of the Caving
Zone. In this study, the double-yield constitutive model in FLAC3D is employed to depict the stress recovery behavior of collapsed rock mass, which considers the permanent volumetric strain under isotropic pressure [25,26]. Two kinds of parameters of cap pressure and material properties are needed in the double-yield constitutive model. And the material properties include density, bulk modulus, shear modulus, internal friction angle, and dilatancy angle. e cap pressure can be calculated by a theoretical equation, while the material parameters can be determined by the "trial-and-error" method, as shown in Figure 4.
For the 31108 working face, the bulking coefficient is 1.3, and the uniaxial compression strength is 41 MPa. Based on equations (2)-(3), the maximum strain and the initial modulus of the gob materials can be estimated as 0.23 and 66 MPa, respectively. Furthermore, the cap pressures for the double-yield model are obtained and plotted in Figure 5.
A trial-and-error method was employed to determine the gob materials parameters by matching; the stress-strain curve is obtained by numerical modelling to that found by equation (1) [16]. For this purpose, a single-element submodel with dimensions 1 m × 1 m × 1 m was generated. A constant velocity was applied to the top of the model in the negative z-direction to generate vertical loading on the model. e velocity magnitude was set at 10-5 m/s. e displacement of the four vertical planes of the model was restricted in the normal direction, and a zero vertical displacement condition was set at the base of the model. e input parameters were calibrated by an iterative change in the bulk modulus, shear modulus, the angle of dilation, and the angle of friction of the gob materials. e final material parameters are listed in Table 1.

Development Law Analysis of Water Flowing Fractured Zone in Overlying Strata
As we all know, there are many bedding, discontinuities, and joints in coal or rock mass, which affect their failure characteristics. Although the finite difference method using FLAC3D is incapable of truly capturing their characteristics, coal or rock mass properties that are upscaled from intact   cores properties using the Hoek-Brown criterion can be performed in the model. In other words, the coal and rock properties in our model are effective properties that accounted for the rock discontinuities. erefore, in this section, we use FLAC3D software to analyze the height of the caving zone and fracture zone in the 31108 working face [2,27]. e double-yield model was used for the gob modelling, and the strain-softening Mohr-Coulomb model was used for the coal and rock mass modelling. e mechanical properties for the double-yield model and strain-softening model are listed in Tables 1 and 2, respectively. Note that the mechanical properties of the coal and rock mass are calculated by the RockLab software, which is based on the generalized Hoek-Brown strength criterion [27][28][29][30][31][32].

Simulation Scheme and Process. Previous studies have
shown that tensile failure mainly occurs in the caving zone, while shear failure mainly occurs in the fracture zone [33]. Accordingly, during the numerical simulation process, we distinguish the caving and fracture zone by the tensile and shear strain value by a self-developed FISH language. If the plastic tensile strain value of the zone is larger than 0.05, it is judged as a caving zone, while it is judged as a rupture zone if the shear strain value is larger than 0.05. e modelling process is as follows. After the initial balance of the numerical model, preliminary excavation is conducted; then, the range of the caving and fracture zones is determined by the plastic tensile and shear strain values; meanwhile, the caving zone is defined as a double-yield model and calculated to equilibrium. e above modelling process is conducted circularly until the working face is mined out. Figure 7 shows a threedimensional evolution diagram of the failure zone during the working face retreating, and Figure 8 shows the volume statistical analysis diagram of failed zones in two zones.

Analysis of Simulation Results.
It can be inferred in Figures 7 and 8 that when the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume of failed zones is 6.6 times, 7.59 times, 9.11 times, and 9.89 times that of the stope. Meanwhile, with the working face retreating, the volume of the failed zones increases quickly first and then tends to be stable. And the ultimate volume ratio of the failed zones and the stope maintains at 9-10. When the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume ratio of the fracture zone to the caving zone is 1.14, 1.74, 1.72, and 1.80, respectively. It can be observed clearly that the volume ratio changes tend to be stable with the working face retreating and finally maintains 1.7-1.8 when the working face retreats to the "square" stage (the panel retreat length equates to the dip length). Figure 9 is the height evolution law of the two zones at various retreating distances. It can be observed that the height evolution law of the caving and fracture zones kept a similar tendency, that is, it first increases rapidly and then tends to be stable with the working face retreating. When the working face retreats 25 m, the heights of the caving and fracture zones are 6 m and 10 m, respectively. And when the working face retreats 50 m (square stage), the height of the caving and fracture zones are 8 m and 16 m, increased by 33% and 60%, respectively. And furthermore, the increase

eoretical Calculation and Field Measurement.
In the longwall mining process, the height of the two zones reflects the movement characteristics of the overlying strata. erefore, it is helpful to master the height of the two zones for the understanding of the movement law of overlying strata.
Based on the field measured data, the calculation method of the heights of two zones under different geological conditions is deduced by Bai et al. [34]. e obtained formula is as follows: where σ 1 is the uniaxial compressive strength; H c is the height of the caving zone; H f is the height of the fracture zone; c 1 , c 2 , c 3 , and c 4 are the rock strength coefficients, as given in Table 3. e mining height of the 31108 working face is 2.9 m. According to the uniaxial compressive test of the roof strata, the average strength of strata in the roof is 15 MPa, which belongs to the weak rock layer. Substituting the parameters into equation (5), the heights of the caving and fracture zones are calculated to 8.5 m and 20.7 m, respectively.
Studies have shown that field observation of borehole water injection leakage can obtain the height of two zones accurately [11]. In the field, a test of borehole water injection leakage was performed in a nearby working face, which has similar geological conditions to 31108. e observation results show that the height of the caving and fracture zone is 11 m and 21.5 m, respectively.

Result from Analysis.
e heights of the caving zone and fracture zone obtained by numerical simulation, theoretical calculation, and field measurement are given in Table 4. It       Shock and Vibration can be observed clearly that the results obtained by simulation are less than that by the field measurement. is may be attributed to the great simplification in terms of geological conditions in the numerical model. However, within the allowable error range, the measured results and simulation results are almost the same. It is well known that the failure mode and height of the strata overlying the stope are affected by various factors in the field, such as mining height, retreating rate, and buried depth.
is study is only based on a specific geological condition of Lijiahao coal mine, and more studies are needed to prove the reliability of the numerical calculation method. However, the modelling procedure and evaluation method presented in this study are helpful to the other coal mines.

Conclusions
(1) Considering the goaf compaction effect, a new numerical simulation method by using FLAC3D is proposed to judge the fracture zone and caving zone of the overlying strata. e method can accurately predicate the height of two zones. (2) With the strain increasing, the vertical stress of the collapsed rock mass increases slowly; and then, it increases rapidly at an exponential function when the strain reaches a certain critical value. e critical strain increases with the increase in the crushing expansion coefficient.
(3) e simulation results show that the volume of the fracture zone in the 31108 working face is 1.7-1.8 times that of the caving zone, while the rock failure is 9-10 times that of the stope. e simulated height of the caving zone and fracture zone is 9 m and 20 m, respectively.

Data Availability
e data used to support the findings of this study are included within the article.