Deformation prediction analysis of sandstone coal mine goaf based on discrete element method and probability integration method

Accurate prediction of goaf settlement is essential for evaluating surface movement and deformation and for understanding overburden failure mechanisms in sandstone coal mine goafs. Coal seam extraction commonly induces significant overburden failure, stress redistribution, surface settlement, and ground fissures, which pose threats to mining safety and ecological stability, particularly under conditions of thick loose strata and thin bedrock. In this study, a sandstone coal mine in Ningxia, China, is selected as the research site. A three-dimensional discrete element model is established using 3DEC to simulate overburden deformation during coal seam extraction under complex geological conditions. The evolution of plastic zones, displacement fields, stress distribution, and surface settlement at different mining stages is systematically analyzed. Meanwhile, the probability integration method is employed to conduct dynamic prediction of surface settlement, with key parameters determined through inversion of field monitoring data. The results indicate that coal seam mining results in an elliptical settlement basin, with a maximum surface settlement of 3.77 m occurring at the goaf center. Shear failure is dominant along coal pillar edges, while significant stress concentration develops in the central goaf area. The settlement predictions obtained from the probability integration method show strong agreement with numerical simulation results, with prediction errors controlled within 3.6%–6.0%. These findings demonstrate that the combined application of discrete element numerical simulation and the probability integration method is effective for predicting surface settlement and analyzing overburden failure mechanisms under thick loose strata and thin bedrock conditions. The study provides a robust theoretical basis and practical reference for surface stability evaluation, hazard prevention, and ecological restoration in coal mining areas.

1 Introduction

The study of overburden failure characteristics and surface movement and deformation in coal mine goafs is a fundamental issue for ensuring safe coal extraction and effective land reclamation. With the intensification of deep coal mining, the impacts of different mining methods on overburden structures and surface responses become increasingly pronounced (Ren et al., 2021; Zhang et al., 2019; Guo et al., 2020). In room-and-pillar mining, coal pillars of specific widths are retained within the seam to support the overlying strata, maintaining short-term overburden stability. However, under long-term influences such as groundwater seepage, oxidation and mining-induced disturbances, coal pillars gradually degrade and eventually fail, leading to surface settlement and secondary geological hazards (Zhao et al., 2021; Ding et al., 2023; Zai et al., 2020; Li et al., 2020; Zhang et al., 2018). In contrast, longwall mining removes large, continuous panels of coal, causing the overburden structure to lose support and progressively collapse under gravity, forming caving zones, fracture zones and bending settlement zones that together constitute a complex caving structural system (Han, 2020; Wang P et al., 2018; Guo et al., 2019; Zhu et al., 2019). Based on extensive field observations and monitoring data, researchers have proposed various classification models of overburden failure structures, including the granular structure model, block-fracture model and the “lower caving–middle arch–upper bending” spatial structural form, which have significantly advanced understanding of overburden evolution under different mining conditions (Kudłacik et al., 2021; Jiang et al., 2019; Cui and Deng, 2017; Chen et al., 2020; Liu et al., 2021). Moreover, studies have demonstrated that the stress distribution and transmission characteristics of the overlying strata in goafs, under the combined effects of mining-induced disturbance and additional surface loading, exert a substantial influence on surface settlement and the development of secondary hazards (Jones and Hobbs, 2021).

Regarding surface movement and deformation, the theoretical framework has evolved from the early empirical “vertical line theory” to the probability integration model (Li, 2024; Zhang et al., 2021; Wang et al., 2021). This approach simulates surface settlement induced by underground goafs through the construction of an influence function, offering high computational efficiency and broad applicability, and has become one of the primary methods for predicting surface settlement in China (Bai et al., 2025a; Bai et al., 2025b). However, under complex geological conditions—characterized by intricate structures, strong stratigraphic heterogeneity, and intensive mining activities—the accuracy of traditional methods is limited, which has driven the rapid development of next-generation deformation monitoring technologies. In recent years, the integrated application of multi-source monitoring techniques, including InSAR (Interferometric Synthetic Aperture Radar), GNSS (Global Navigation Satellite System), and oblique photogrammetry (Bai et al., 2023), has significantly enhanced the spatiotemporal resolution and accuracy of surface settlement monitoring. These advances enable real-time, high-precision deformation control, particularly in densely populated urban areas and around critical engineering infrastructures (Zhang et al., 2025). Against this background, numerical simulation methods increasingly became indispensable for analyzing surface subsidence induced by coal-mine extraction (Xie et al., 2025), for instance, utilized the FLAC3D numerical simulation software to conduct a systematic investigation of surface subsidence associated with the mining of a working face in a coal mine in North China. Their results demonstrated that the maximum surface subsidence was negatively correlated with mining depth but positively correlated with coal-seam thickness. The resulting subsidence basin exhibited marked asymmetry, with the advance angle smaller than the lag angle, and the evolution of surface-movement velocity displayed a characteristic “S-shaped” trajectory as mining progressed. (Li et al., 2025). performed UAV-based monitoring of surface subsidence in a coal mine in Northwest China and carried out numerical simulations using FLAC3D. The study identified pronounced subsidence and collapse features, with a maximum subsidence reaching 3.78 m. The numerical results further elucidated the stress distribution within the goaf and the dynamic evolution of subsidence along the principal section, thereby confirming the model’s reliability (Wang et al., 2024). Investigated the 2311 working face of the Changcun Coal Mine using the probability-integration method in combination with discrete-element simulation. Their findings indicated that grouting-assisted mining effectively constrained the magnitude of surface subsidence, thus mitigating ground deformation and enhancing coal-resource recovery.

Building on the aforementioned research, this study further adopted an integrated framework that combines numerical simulation with the probability-integration method. A three-dimensional finite-element numerical model was developed according to the geological characteristics and mining conditions of the goaf to simulate the evolution of surface subsidence under complex stratigraphic settings involving thick coal seams, thin bedrock, and loose overburden. Comparative analysis between the numerical simulation results and those derived from the probability-integration method verified the reliability and consistency of the model in predicting the subsidence extent, the spatial distribution of subsidence magnitudes, and the position of the maximum subsidence center. The findings provide a robust theoretical foundation and practical engineering reference for advancing the understanding of subsidence mechanisms in goaf areas, evaluating surface stability, and supporting ecological restoration in mining regions.

2 Project overview

The coal mine is located in Ningdong Town, Lingwu City, Ningxia Hui Autonomous Region. The planar location of the coal mine is shown in Figure 1, The mining lease extends approximately 10 km from north to south and 6 km from east to west, covering an area of about 60.13 km². The total coal resources amount to approximately 1.358 billion tonnes, of which the proved mineable reserves reach 759 million tonnes. Within the lease area, about twenty numbered coal seams are identified, among which six are classified as mineable. The cumulative mineable coal-seam thickness is approximately 26.95 m, corresponding to a mineable coal-bearing coefficient of about 9.18%. The principal mineable seams include the No. 2, No. 8, No. 17, and No. 18 coal seams. Among these, the No. 2 seam—currently under extraction—is positioned at the top of the eighth minor cyclothem within the fifth member of the Yan’an Formation, representing the uppermost mineable coal-bearing horizon in the stratigraphic sequence. Representing the uppermost recoverable coal seams. The operation area adopts a retreat longwall mining method with full-seam extraction in a single slice, and the roof is managed by full caving. Within the mining area, this coal seam is extensively distributed, covering a total area of 48.52 km², all of which is recoverable, yielding a 100% extractable probability. The roof strata primarily consist of thick-bedded coarse-grained sandstone of the Zhiluozhen Formation (Qilizhen Sandstone), while in some areas they are composed of thin-bedded fine-grained sandstone or siltstone. The seam is classified as a typical thick coal seam, exhibiting relatively uniform thickness and a simple structural configuration. The coal type is dominated by non-caking coal, followed by long-flame coal, both of which demonstrate high stability and reliability. As of 2017, the test site of the mining field spans 130.22 km² and is characterized by significant ground fissures and surface collapses. These fissures extend for 1–2 km, reach widths of up to 2 m and penetrate several meters deep. Such deformations constitute severe geological hazards, posing substantial risks to the local ecological environment, public safety, and mining operations.

Figure 1

Figure 1. Plan view of the test site.

3 Simulation results of the settlement for the site by using 3DEC

3.1 Establishment of the numerical model for the coal mine goaf

3.1.1 Structural simplification of coal mine operation areas and discretization of numerical models

The stability of coal mine operation areas requires a comprehensive assessment of the influence of mining-induced effects under current extraction conditions on the overall stability of the mining area. It is essential to simulate the response of overlying strata and the ground surface to the formation of areas, with the simulation domain determined based on the maximum expected displacement or deformation. Considering the coal seam occurrence conditions, borehole data and comprehensive columnar sections and in combination with the engineering conditions of operation areas110201–110208, the simulation domain is selected within operation areas 110201 and 110202, as illustrated in Figure 2. Operation area 110201 has an average coal seam dip of 25°, a strike length of 1,100 m, a dip length of 96 m and a mining height of 5.5 m; Operation area 110202 has an average coal seam dip of 24°30′, a strike length of 1,180 m, a dip length of 260 m and a mining height of 5.5 m. The numerical model covers a volume of 600 m × 500 m × 300 m. Given that the average burial depth of the two operation areas is approximately 600 m, an additional 300 m of overlying strata is included in the model, with a self-weight stress of 7.5 MPa.Boundary conditions play a critical role in numerical simulations, directly influencing the reliability and accuracy of the results; therefore, proper constraints are applied. In the present three-dimensional numerical model, normal displacement constraints in the x-direction are applied to the two lateral x-surfaces, normal displacement constraints in the y-direction are applied to the two lateral y-surfaces and the bottom surface (z = 0) is fixed in all directions. Stress boundary conditions are applied at the top surface. The resulting mechanical model is depicted in Figure 2.

Figure 2

Figure 2. Mechanical model and operation area layout.

The discretization of the numerical model is a fundamental step for achieving accurate and effective simulation, exerting a significant influence on the outcomes of the analysis. The 3DEC simulation adopts a stepwise excavation sequence, consistent with the actual panel advancement process. In accordance with the principles of geometric and physical approximation, the coal mine calculation model is discretized, Mesh size (≈120000 zones), Computational time (14.6 h on a 32-core workstation). Producing the three-dimensional numerical model illustrated in Figure 3. From top to bottom, the stratigraphic sequence consists of coarse sandstone, medium sandstone, fine sandstone, sandy mudstone, siltstone, coal, siltstone, medium sandstone, fine sandstone, medium sandstone, and siltstone.

Figure 3

Figure 3. Schematic diagram of the three-dimensional model.

A total of 748 monitoring points were arranged at 20 m intervals on the model surface (ground surface), forming a comprehensive monitoring network to ensure accurate measurement of surface deformation. To clearly present the numerical simulation results, two data sections were established, as shown in Figure 4, located at y = 112 m and y = 305 m, corresponding to the central positions of the 110201 and 110202 working faces, respectively. These sections enable a clear visualization of the stress and displacement fields within the overlying strata during the fully mechanized top-coal caving extraction of the 110201 and 110202 working faces.

Figure 4

Figure 4. (a) Location schematic of the y = 112 m cross-section (b) Location schematic of the y = 305 m cross-section.

3.1.2 Mechanical parameters of the rock mass medium

The geological conditions of the coal mine are relatively complex, with diverse rock structures and types forming a heterogeneous geological medium. Accordingly, the rock mass medium in the numerical model is appropriately simplified by consolidating rocks with similar mechanical properties to obtain a manageable combination of rock types. Based on the geological characteristics of the mining area and results from rock mechanics tests, the rock mass along the dip of the ore body is simplified into six fundamental types: coal, siltstone, fine-grained sandstone, medium-grained sandstone, coarse-grained sandstone and sandy mudstone. Considering laboratory test results as well as the effects of rock structure and scale, Scale effects were corrected using Hoek–Brown disturbance factor D, Bedding and jointing effects were incorporated by reducing laboratory cohesion and friction angle by 10%–20% according to engineering empirical guidelines. The mechanical parameters of each rock layer are determined, as summarized in Table 1. All mechanical parameters were derived from laboratory triaxial compression tests conducted on representative rock samples from the study area. When direct test data were unavailable, supplementary values were taken from peer-reviewed literature with similar lithology and geological settings.

Table 1

Table 1. Mechanical parameters of coal, loose overburden and surrounding rock.

3.2 Distribution characteristics of plastic zones in the test site

Plasticity is a mechanical state of a material, representing its capacity to undergo stable and permanent deformation under external loads without compromising structural integrity. For most engineering materials, when the applied stress is below the yield limit, the material exhibits elastic deformation, whereby the deformation fully recovers once the external load is removed. When the applied stress reaches the yield limit, the material enters a plastic deformation stage, which is irreversible, although the material’s integrity remains intact. As the load continues to increase, the material ultimately fails upon reaching its strength limit, accompanied by a rapid decrease in stress. Thus, the plastic state refers to the condition in which the material has attained its yield limit but has not exceeded its strength limit, serving as a mechanical indicator for assessing the degree of material damage.

The 3DEC numerical simulation software provides a function for identifying plastic zones, enabling the assessment of the model’s failure state based on the spatial distribution of these zones. Its suitability for brittle rock failure in mining conditions, its validated performance in simulating shear-dominated pillar failure, and its wide acceptance in coal-mine goaf studies. During post-processing in 3DEC, the distribution of plastic zones can be visualized and analyzed. The software defines six mechanical states for the model: “none” indicating that the region has not experienced plastic failure and remains in the elastic stage; “none-p” indicating that the region has previously remained entirely elastic and has not experienced any form of failure. “shear-n” indicating that the region is currently undergoing shear failure; “shear-p” indicating that the region has previously undergone shear failure; “tension-n” indicating that the region is currently undergoing tensile failure; and “tension-p” indicating that the region has previously undergone tensile failure.

Following coal seam extraction, the distribution characteristics of plastic zones in the operation area are illustrated in Figures 5, 6. As shown in Figure 5, the coal pillars in the section undergo shear failure, with plastic zones surrounding the pillars and a core zone located at the center of each pillar; some coal pillars, however, do not exhibit a core zone. In addition, the edges of the coal seam along the boundaries of the goaf also experience shear failure.

Figure 5

Figure 5. Distribution of coal pillar plastic zones (plan view).

Figure 6

Figure 6. Distribution of plastic zones in the overlying strata of the operation area. (a) Operation area 110201 (The corresponding section line for the figure is shown in Figure 4a) (b) Operation area 110202 (The corresponding section line for the figure is shown in Figure 4b).

The distribution characteristics of plastic zones in the overlying strata of the operation area are illustrated in Figure 6. At the floor of the goaf, a discontinuous tensile-shear mixed failure zone develops. Following coal extraction, the floor enters a stress-relieved state, resulting in floor heaving and tensile failure within the floor strata. A comparison of the plastic zone distributions in operation areas 110201 and 110202 indicates that both exhibit similar patterns, failure zones are less extensive in the upslope direction and more pronounced in the downslope direction. However, the overall failure area in operation area 110202 is larger, primarily manifested by more extensive shear failure zones in both the overlying strata and the floor region.

3.3 Distribution characteristics of the displacement field in the test site

The vertical displacement distribution of the goaf roof following coal seam extraction is presented in Figure 7. As shown, an elliptical subsidence basin develops within the goaf, with the coal pillars in the section undergoing settlement. The vertical displacement of coal pillars within the elliptical zone is significantly greater than that in surrounding areas, with the maximum settlement at the center of the goaf reaching 3.77 m. The section coal pillars, together with the boundary coal pillars at the edges of the extraction area, collectively support the gravitational load of the overlying strata. Coal pillars located farther from the goaf center, being closer to the boundary pillars, experience lower stress concentrations as the boundary pillars bear the majority of the roof load, resulting in relatively smaller settlement in these areas. The coal pillars radiate outward from the center of the goaf, forming an elliptical test site.

Figure 7

Figure 7. Displacement cloud map of the operation area roof (plan view).

As illustrated in Figure 8, surface settlement forms a basin, with the maximum settlement occurring at the basin center, reaching 3.77 m. The settlement gradually decreases with increasing distance from the basin center along the same direction. Within the operation area, settlement decreases from the bottom upward, with the largest displacement occurring in the strata at the center of the operation area and the smallest near the edges of the goaf. The overlying strata above the coal pillars exhibit relatively minor displacements. To further illustrate the overall movement characteristics of the overburden strata in the vertical direction, the front-view displacement distributions are shown in Figure 9. In operation area 110201, the maximum settlement of the basic roof is only 1.01 m, and the settlement is relatively uniform. In contrast, in operation area 110202, the maximum settlement of the basic roof reaches 3.75 m, with uneven distribution: greater settlement occurs in the upper-left region of the operation area, while smaller settlement occurs in the lower-right region.

Figure 8

Figure 8. Surface displacement contour of the operation area.

Figure 9

Figure 9. Overall movement of the mined overburden (front view). (a) Operation area 110201 (The corresponding section line for the figure is shown in Figure 4a) (b) Operation area 110202 (The corresponding section line for the figure is shown in Figure 4b).

Based on the above analysis, following coal extraction, the overlying strata of the operation area undergo settlement. In operation area 110202, the overlying strata experience relatively greater settlement, whereas the strata above the section coal pillars exhibit smaller settlement, leading to localized uneven settlement of the roof strata, although the affected area is limited. Overall, the settlement of the overlying strata gradually decreases from the bottom upward, with the maximum displacement occurring at the center of the operation area and smaller displacements toward the periphery. Near the surface, the settlement tends to be gentle, and the overall magnitude of settlement remains minimal.

3.4 Distribution characteristics of the stress field in the test site

The vertical stress distribution throughout the operation area following coal extraction is presented in Figures 10, 11. As shown in Figure 10, coal mining disrupts the original in-situ stress equilibrium of the strata, inducing a secondary stress redistribution in the surrounding rock of the operation area. After coal extraction, the section coal pillars and the coal walls on both sides of the goaf support the full load of the overlying strata, resulting in relatively high vertical stress in the coal pillars at the center of the operation area. In contrast, the stresses in the roof and floor decrease significantly compared to the pre-mining state, placing these regions within a stress-relieved zone.

Figure 10

Figure 10. Vertical stress distribution (front view). (a) Operation area 110201 (The corresponding section line for the figure is shown in Figure 4a) (b) Operation area 110202 (The corresponding section line for the figure is shown in Figure 4b).

Figure 11

Figure 11. Stress distribution of coal pillars (front view).

4 Calculation results of the settlement for the site by using the probability integration methods

4.1 Fundamental principles of the probability integration methods

The probability integration method considers that surface strata movements induced by mining are random events. From a statistical perspective, under any mining condition, the entire extraction process can be divided into numerous small mining units, and the overall impact on the strata and surface is equal to the cumulative effect of all individual units. Experiments indicate that under ideal conditions the profile of a settlement basin follows a normal distribution, consistent with the distribution of probability density. The probability integration method uses the normal distribution function as an influence function and employs integration to represent the profile of surface settlement basins.

The surface settlement rate in the mining-affected area decreases over time following a negative exponential curve, primarily influenced by factors such as mining depth, the lithology of the overlying strata and the mining rate. Additionally, it represents the cumulative impact of different mining blocks on the surface at various times, under differing geological and mining conditions and employing different mining methods. The mathematical model used to predict surface movement and deformation is expressed as Equation 1:

$Y\left(\mathrm{x},\mathrm{y},\mathrm{z}\right)={\displaystyle \sum _{i=1}^j}{C}_i{\displaystyle \sum _{i=1}^k}{\displaystyle \underset{qk}{\overset{qk+1}{\int }}}{F}_1\left(Rk,z\right)F_2\left(q\right)dq(1)$

In the equation: (x, y, z) denotes the coordinates of the calculation point; C_i represents the time influence coefficient; j indicates the number of calculation blocks; k indicates the number of inflection points for any mined block; q_k is the angle between the x-axis of the Cartesian coordinate system used for calculation and the line connecting the calculation point with the inflection point h; F₁ and F₂ are computational functions; R_k = R_k (q, x, y) represents the polar coordinate radius.

$R_k=\frac{\left(x_k-x\right)\left(y_{k+1}-y_k\right)-\left(y_k-y\right)\left(x_{k+1}-x_k\right)}{\left(y_{k+1}-y_k\right)\cos q-\left(x_{k+1}-x_k\right)\sin q}(2)$

Based on the model presented in Equation 2, the predicted settlement at point p (x, y) within the mining area is derived and expressed as shown in Figure 12:

$w\left(x,y,t\right)=\frac{w_{\max }}{2\pi }{C}_i·{\displaystyle \sum _{i=1}^k}{\int }_{q_k}^{q_{k+1}}\left[1-e^{-\pi \frac{R_k^2}{r^2}}\right]dq(3)$

Figure 12

Figure 12. Schematic diagram of polar coordinate closed integration.

In the Equation 3: w_max represents the maximum settlement resulting from full mining-induced deformation; C_i = (1−e^cti) denotes the time-dependent settlement coefficient, which is determined by mining depth, the lithology of the overlying strata and the mining rate.

The value of c is defined as follows: when the mining depth is shallow and the overlying strata are loose and soft: c = 2.5–3.0; when the mining depth is shallow and the overlying strata are relatively hard: c = 2.0∼2.5; when the mining depth is large and the overlying strata are soft: c = 1.5∼2.0; and when the mining depth is large and the overlying strata are hard: c = 1.0∼1.5. Under repeated mining-induced deformation, the value of c is generally less than 1.

The expressions for the components induced by predicted surface dynamic settlement, including inclination (i_x, i_y), curvature (K_x, K_y), horizontal movement (u_x, u_y), and horizontal deformation (ε_x, ε_y), are as Equations 4–10:

$i_x=\frac{w_{\max }}{r}{\displaystyle \sum _{i=1}^j}{C}_i{\displaystyle \sum _{i=1}^k}{\displaystyle \underset{q_k}{\overset{q_{k+1}}{\int }}}{f}_1\left(\frac{R_k}{r}\right)\cos qdq(4)$

$i_y=\frac{w_{\max }}{r}{\displaystyle \sum _{i=1}^j}{C}_i{\displaystyle \sum _{i=1}^k}{\displaystyle \underset{qk}{\overset{q_{k+1}}{\int }}}{f}_1\left(\frac{R_k}{r}\right)\sin qdq(5)$

$K_x=\frac{w_{\max }}{r^2}{\displaystyle \sum _{i=1}^j}{C}_i{\displaystyle \sum _{i=1}^k}\left(P_k-Q_k\right)(6)$

$K_y=\frac{w_{\max }}{r^2}{\displaystyle \sum _{i=1}^j}{C}_i{\displaystyle \sum _{i=1}^k}\left(P_k+Q_k\right)(7)$

$u_x=w_{\max }{\displaystyle \sum _{i=1}^j}{b}_i{C}_i{\displaystyle \sum _{i=1}^k}{\displaystyle \underset{q_k}{\overset{q_{k+1}}{\int }}}{f}_1\left(\frac{R_k}{r}\right)\cos qdq(8)$

$u_y=w_{\max }{\displaystyle \sum _{i=1}^j}{b}_i{C}_i{\displaystyle \sum _{i=1}^k}{\displaystyle \underset{q_k}{\overset{q_{k+1}}{\int }}}{f}_1\left(\frac{R_k}{r}\right)\sin qdq(9)$

${\epsilon }_x=\frac{w_{\max }}{r^2}{\displaystyle \sum _{i=1}^j}{C}_i{b}_i{\displaystyle \sum _{i=1}^k}\left(P_k-Q_k\right)(10)$

Equations 11–15 define and explain the symbols appearing in Equations 4–10.

$f_1\left(s\right)=s{e}^{-\pi s^2}-{\displaystyle \underset{0}{\overset{s}{\int }}}{e}^{-\pi t^2}dt(11)$

$f_2\left(s\right)=\pi s^2{e}^{-\pi s^2}(12)$

$f_3\left(s\right)=1-e^{-\pi s^2}-\pi s^2{e}^{-\pi s^2}(13)$

$P_k={\displaystyle \underset{q_k}{\overset{q_{k+1}}{\int }}}{f}_2\left(\frac{R_k}{r}\right)dq(14)$

$Q_k={\displaystyle \underset{q_k}{\overset{q_{k+1}}{\int }}}{f}_3\left(\frac{R_k}{r}\right)\cos 2qdq(15)$

4.2 Parameter determination

When predicting surface movement and deformation induced by mining using the probability integration method, it is first necessary to determine five predictive parameters: the settlement coefficientη, the main influence angle tangent tanβ, the horizontal movement coefficient b, the inflection point offset distance d, and the mining influence propagation angle θ, all of which are obtained through analysis of measured data. The most critical parameter is the maximum settlement. Under conditions of full mining-induced deformation (ensuring full mining is achieved in both the strike and dip directions), the maximum settlement along each survey line is determined from observational data collected at movable monitoring stations.

4.2.1 Settlement coefficient η

Since full mining-induced deformation is achieved in both the strike and dip directions of the operation area, the maximum settlement observed along the survey lines of the principal strike section and the principal dip section represents the maximum settlement under conditions of full mining-induced deformation.

The settlement coefficient is determined based on the relationship between the maximum settlement at the measurement points and the mining height, expressed as η = (W_max/m)/cosα, with the settlement coefficient in Table 2.

Table 2

Table 2. Maximum surface settlement value and settlement coefficient.

From the perspective of protecting the ground surface, the surface settlement coefficient for coal mining is assigned a value of 0.55.

4.2.2 The main influence angle tangent tanβ

The magnitude of the tangent of the main influence angle tanβ is closely related to the lithology of the overlying strata and serves as a parameter reflecting the relationship between the overlying strata and the mining depth. Softer overlying strata correspond to larger tanβ values, whereas harder strata correspond to smaller values. The angle β formed between the line connecting the boundary points of the main influence zone and the mining boundary and the horizontal is defined as the main influence angle.

The main influence angles along the surface Line A and Line B of the 110202 fully mechanized operation area are approximately equal. Therefore, the moving angle is adopted as the main influence angle for both the Line A and Line B, with β = 75.4° and the corresponding tangent of the main influence angle tanβ = 3.84.

4.2.3 Inflection point offset distance d

When the mined space undergoes full mining-induced deformation in both directions, the point where the settlement reaches 0.5 Wmax (at which both inclination i and horizontal movement u attain their maximum values) is defined as the inflection point of the settlement basin. The distance from this inflection point to the nearest coal wall is referred to as the inflection point offset distance. The magnitude of the offset distance is influenced by mining depth, the lithology of the overlying strata, and the hardness of the coal seam; greater mining depth and harder overlying strata or coal seam result in a larger offset distance, and vice versa. Based on data from the B observation line, the inflection point corresponding to 0.5 Wmax yields an inflection point offset distance of d = 43.8 m, as illustrated in Figure 13.

Figure 13

Figure 13. Inflection point offset distance along B survey line.

4.2.4 Horizontal movement coefficient b

The horizontal movement coefficient represents the ratio between the maximum horizontal movement and the maximum settlement. Based on observational data from the two survey lines, the calculated horizontal movement coefficients are summarized in Table 3. The final horizontal movement coefficient, determined from a lateral line on the principal dip section is 0.44.

Table 3

Table 3. Table of horizontal movement coefficient determination.

4.3 Surface settlement and influencing factors analysis

The surface movement resulting from the mining of coal seams beneath thick loose strata and thin bedrock exhibits certain regularities. However, the surface movement and deformation induced by mining settlement is highly influenced by the distribution of the loose surface layers, and under general conditions, its regularity is relatively low. The predicted settlement results using the probability integration method are presented in Figure 14. In Figure 14, the 110202 working face has dimensions of 600 m × 260 m. For the 110202 operation area, under conventional mining conditions, the maximum surface settlement caused by coal seam extraction is 3.601 m. During the actual mining process, the settlement varies with changes in mining thickness. The theoretically predicted inflection point offset distance is 43.8 m.

Figure 14

Figure 14. Predicted settlement results of the 110202 operation area using the probability integration method.

4.4 Analysis of surface horizontal movement and deformation

Figures 15, 16 illustrate the predicted distribution patterns of surface horizontal deformation for the 110202 operation area. All coordinate values, subsidence measurements, and horizontal displacement data from the monitoring points are imported and processed. Through coordinate transformation, these parameters are converted into functional forms, after which the mining-induced subsidence prediction parameters are solved. The known parameters and the estimated parameters are then substituted into Equations 8, 9 of the probability integration method to calculate the surface deformation within the mining-induced subsidence area. Using the “ezmesh” and “ezcontour” functions in Matlab, the subsidence surface and subsidence contour maps are generated. This approach enables the prediction or inversion of surface deformation, and the corresponding contour diagrams are subsequently produced. As shown, the maximum horizontal movements along the advance direction of the operation area are −1834 mm and +2014 mm, while the maximum horizontal movements perpendicular to the advance direction reach −2042 mm and +2255 mm.

Figure 15

Figure 15. Surface horizontal movement along the advancing direction of operation area 110202 (m).

Figure 16

Figure 16. Surface horizontal movement perpendicular to the advancing direction of operation area 110202 (m).

In this study, the “±” signs for horizontal movement values along the operation area advance direction indicate direction only, not magnitude: “+” denotes the advance direction of the operation area, while “−” denotes the opposite direction. Likewise, for horizontal movement values perpendicular to the advance direction, the “±” signs represent direction only: “+” indicates movement from the return air roadway toward the center of the goaf, and “−” indicates movement from the belt roadway toward the center of the goaf.

4.5 Comparison of numerical simulation and probability integration method results

Cross-sectional profiles along the A and B observation lines are extracted from Figure 14 to plot the predicted settlement curves in these directions. The comparison between numerically simulated settlement and predicted settlement along the A and B observation lines is presented in Figures 17, 18, respectively, providing a basis for evaluating the fit between the selected theoretical prediction model and parameters and the numerical simulation results. As shown in Figure 17, the numerically simulated surface settlement profiles generally align with the theoretically predicted curves, although some deviations exist due to the following reasons:

1. Mining induces shear failure, resulting in step-like fractures along the boundaries of the goaf. The geological conditions and mining techniques of the 110202 operation area cause discontinuities in the surface movement and deformation.

2. Near the shutdown line of the 110202 operation area, the surface is overlain by a relatively thick loose sand layer. This layer absorbs mining-induced energy, effectively reducing the magnitude of settlement and extending the duration of residual settlement.

Figure 17

Figure 17. Comparison of numerical simulation and theoretical predicted surface settlement curves along observation Line A.

Figure 18

Figure 18. Comparison of numerical simulation and theoretical predicted surface settlement curves along observation Line B.

The comparison of Figures 17, 18 shows that when no step-like settlement occurs at the surface, the numerically simulated curves closely match the predicted curves. Conversely, when step-like settlement is present, notable discrepancies arise between the numerical simulation and the predicted results. Additionally, the buffering effect of the loose sand layer significantly influences the predicted surface settlement.

5 Conclusion

This study addresses critical scientific issues surrounding overlying strata failure and surface settlement in sandstone coal mine goafs by integrating 3DEC numerical simulation with the probability integration method. The research establishes a comprehensive settlement prediction framework that combines field monitoring, mechanical analysis, and numerical modeling. The main contributions, findings, and limitations of this work are summarized as follows:

1. A high-precision surface settlement prediction system is developed by coupling probability integration theory with parameter inversion based on field measurements. The calibrated parameters—including a settlement coefficient of 0.55, a tangent of main influence angle of 3.84, a horizontal movement coefficient of 0.44, and an inflection point offset of 0.30H₀—are validated against monitoring data from the 110202 operation area. The predictive accuracy, with deviations maintained within 3.6%–6.0% when compared with numerical simulations, demonstrates the robustness and applicability of the probability integration method for mining areas characterized by thick loose deposits and thin bedrock. The novelty lies in the systematic parameter calibration under complex stratigraphic conditions.

2. A three-dimensional 3DEC numerical model is constructed to elucidate the mechanical response of overlying strata during mining. The model captures the formation of an elliptical settlement basin and a maximum subsidence of 3.77 m, consistent with measured values, confirming its predictive reliability. Moreover, the simulation reveals the internal failure mechanism, including shear-dominated yielding at coal pillar edges, elastic retention within pillar cores, and significant vertical stress concentration (up to 37.67 MPa) at the goaf center. These findings provide mechanical insights that extend the understanding of roof instability and stress redistribution in deep mining settings.

3. The integrated probability-numerical prediction framework improves the accuracy of settlement basin delineation under conditions of thick loose overburden and thin bedrock, which are typical geological settings in many northern Chinese mining areas. The refined deformation patterns, especially the identification of maximum settlement centers and boundary gradients, provide more reliable references for evaluating the stability of buildings, roads, and surface infrastructure above goaf zones. The validated prediction parameters (settlement coefficient, main influence angle, inflection point offset) can be directly applied in surface deformation early-warning systems for similar mining conditions. This supports the development of proactive monitoring schemes and improved risk-prevention mechanisms.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

GR: Conceptualization, Investigation, Writing – original draft. CX: Data curation, Formal Analysis, Writing – review and editing. WY: Methodology, Validation, Writing – review and editing. LG: Supervision, Writing – review and editing, Formal Analysis. LJ: Writing – review and editing, Validation, Resources.

Funding

The author(s) declared that financial support was received for this work and/or its publication. Funded by Economic and Technical Research Institute of State Grid Shaanxi Electric Power Company (SGTYHT/24-JS-004).

Conflict of interest

Authors GR, CX, and LJ were employed by Economic and Technical Research Institute of State Grid Shaanxi Electric Power Company. Authors WY and LG were employed by Shaanxi Local Electric Power (Group) Co., Ltd.

The authors declare that this study received funding from State Grid Shaanxi Electric Power Company. The funder had the following involvement in the study: Conduct of this study.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bai, B., Zhou, R., Yang, G., Zou, W., and Yuan, W. (2023). The constitutive behavior and dissociation effect of hydrate-bearing sediment within a granular thermodynamic framework. Ocean. Eng. 268, 113408. doi:10.1016/j.oceaneng.2022.113408

CrossRef Full Text | Google Scholar

Bai, B., Wu, H. Y., Zhou, R., Wu, N., and Zhang, B. X. (2025a). A granular thermodynamic framework-based coupled multiphase-substance flow model considering temperature driving effect. J. Rock Mech. Geotech. Eng. 17, 5816–5828. doi:10.1016/j.jrmge.2024.11.017

CrossRef Full Text | Google Scholar

Bai, B., Zhang, B. X., Chen, H. J., and Chen, P. P. (2025b). A novel thermodynamic constitutive model of coarse-grained soils considering the particle breakage. Transp. Geotech. 50, 101462. doi:10.1016/j.trgeo.2024.101462

CrossRef Full Text | Google Scholar

Chen, J., Zhao, J., Zhang, S., Zhang, Y., Yang, F., and Li, M. (2020). An experimental and analytical research on the evolution of mining cracks in deep floor rock mass. Pure Appl. Geophys. 177 (11), 5325–5348. doi:10.1007/s00024-020-02550-9

CrossRef Full Text | Google Scholar

Cui, X. M., and Deng, K. Z. (2017). Research review of predicting theory and method for coal mining subsidence. Coal Sci. Technol. 45 (1), 160–169.

Google Scholar

Ding, R. L., Guo, Q., Bai, X., Sun, Z. J., and Zhang, B. C. (2023). Research on stability risk assessment of expressway crossing coal mine goaf. Mod. Inf. Technol. 7 (13), 136–140+144.

Google Scholar

Guo, Q. B., Wang, L., Lu, X., and Jiang, C. L. (2019). Stability evaluation of the construction site above a goaf based on cloud model and fuzzy hierarchy analysis. Metal. Mine (11), 43–48.

Google Scholar

Guo, S., Guo, G. L., Li, H. Z., and Yang, X. S. (2020). Goaf-collapse sites stability evaluation based on principal component hierarchical clustering model. Chin. J. Geol. Hazard Control. 31 (06), 116–121.

Google Scholar

Han, K. M. (2020). Theoretical study on stability evaluation of goaf overburden under surface load. Beijing, China: China Coal Industry Publishing House.

Google Scholar

Jiang, Y., Misa, R., Li, P. Y., Yuan, X., Sroka, A., and Jiang, Y. (2019). Summary and development of mining subsidence theory. Metal. Mine 10, 1–7.

Google Scholar

Jones, L., and Hobbs, P. (2021). The application of terrestrial LiDAR for geohazard mapping, monitoring and modelling in the British geological survey. Remote Sens. 13 (3), 395. doi:10.3390/rs13030395

CrossRef Full Text | Google Scholar

Kudłacik, I., Kapłon, J., Lizurek, G., Crespi, M., and Kurpiński, G. (2021). High-rate GPS positioning for tracing anthropogenic seismic activity: the 29 January 2019 mining tremor in Legnica- głogów copper district, Poland. Measurement 168, 108396. doi:10.1016/j.measurement.2020.108396

CrossRef Full Text | Google Scholar

Li, H. W. (2024). Study on fully-mechanized solid-filling coal mining technology. Energy and Energy Conservation 10, 208–211. doi:10.16643/j.cnki.14-1360/td.2024.10.022

CrossRef Full Text | Google Scholar

Li, X. B., Qiu, J. D., Zhao, Y. Z., Chen, Z. H., and Li, D. Y. (2020). Nstantaneous and long-term deformation characteristics of deep room-pillar system induced by pillar recovery. Trans. Nonferrous Metall. Soc. China. 30 (10), 2775–2791. doi:10.1016/s1003-6326(20)65420-6

CrossRef Full Text | Google Scholar

Li, X. D., Zhang, L. J., and Zhao, X. C. (2025). Research on integrated dynamic monitoring and simulation prediction for surface and underground subsidence in coal mining. China Energy Environ Prot. 47 (07), 79–84+94.

Google Scholar

Liu, J., Wu, N., Si, G., and Zhao, M. (2021). Experimental study on mechanical properties and failure behaviour of the pre-cracked coal-rock combination. Bull. Eng. Geol. Environ. 80 (3), 2307–2321. doi:10.1007/s10064-020-02049-6

CrossRef Full Text | Google Scholar

Ren, L. W., Ning, H., Zhou, Y. F., Dun, Z. L., Guo, W. B., and Tian, Z. B. (2021). Research status and prospect on deformation control of high-speed railway subgrade in goaf site. J. China Coal Soc. 46 (08), 2534–2547.

Google Scholar

Wang, Y., Liu, D., Dong, J., Zhang, L., Guo, J., Liao, M., et al. (2021). On the applicability of satellite SAR interferometry to landslide hazards detection in hilly areas: a case study of shuicheng, Guizhou in southwest China. Landslides 18 (7), 2609–2619. doi:10.1007/s10346-021-01648-y

CrossRef Full Text | Google Scholar

Wang, P. F., Miao, H. Z., Chai, J. F., Wang, L., and Wu, H. Y. (2024). Probability integral prediction and numerical simulation of surface subsidence in bed separation grouting and filling mining in changcun coal mine. ShanXi Coal. 44 (03), 132–139+148.

Google Scholar

Wang, P. F., Feng, G. R., Zhao, J. L., Chugh, Y., and Wang, Z. Q. (2018). Effect of longwall gob on distribution of mining-induced stress. Chin. J. Geotech. Eng. 40 (07), 1237–1246.

Google Scholar

Xie, Z. C., Fu, W. G., Guo, J. F., Huang, Y. H., and Lian, G. L. (2025). Surface subsidence characteristics in process of coal mining based on numerical simulation. Energy Energy Conservation 11, 73–76.

Google Scholar

Zai, H., Wang, W., and Shen, W. (2020). Safety and damage assessment method of transmission line tower in goaf based on artificial intelligence. I&CPS Asia, 1474–1479. doi:10.1109/icpsasia48933.2020.9208404

CrossRef Full Text | Google Scholar

Zhang, Y., Xu, Y., Wang, K., Chen, P., Wang, X., Zheng, Q., et al. (2018). The fracturing characteristics of rock mass of coal mining and its effect on overlying unconsolidated aquifer in Shanxi, China. Arab. J. Geosci. 11 (21), 666. doi:10.1007/s12517-018-4034-0

CrossRef Full Text | Google Scholar

Zhang, J., Liu, Q. Z., Yang, T., Zhou, F. W., and Wang, B. (2019). Analysis of “box-girder-bridge” structural for upward fully mechanized mining in room-and-pillar mining goaf of shallow close distance seam. J. Min. and Saf. Eng. 36 (05), 867–872.

Google Scholar

Zhang, Y., Wu, H., Li, M., Kang, Y., and Lu, Z. (2021). Investigating ground subsidence and the causes over the whole Jiangsu Province, China using Sentinel-1 SAR data. Remote Sens. 13 (2), 179. doi:10.3390/rs13020179

CrossRef Full Text | Google Scholar

Zhang, M., Tian, L. H., Xie, J. L., Hu, D. K., Zhou, F. L., Zhang, Y. F., et al. (2025). Experimental study on the swelling performance of GMZ24 bentonite buffer material under Fe (II) reducing environment. Nucl. Eng. Des. 442, 114280. doi:10.1016/j.nucengdes.2025.114280

CrossRef Full Text | Google Scholar

Zhao, B., Zhao, Y., and Wang, J. (2021). New stability forecasting model for goaf slope based on the AHP-TOPSIS theory. Arabian J. Geosci. 14 (1), 17. doi:10.1007/s12517-020-06292-9

CrossRef Full Text | Google Scholar

Zhu, W. B., Xu, J. L., Chen, L., Li, Z., and Liu, W. T. (2019). Mechanism of disaster induced by dynamic instability of coal pillar group in room-and-pillar mining of shallow and close coal seams. J. China Coal Soc. 44 (02), 358–366.

Google Scholar