Prediction and Analysis of Surface Residual Deformation Considering the Impact of Groundwater in Mines

Abstract

With economic development and coal resource exploitation, the area of mined-out zones is expanding continuously. The traditional waste disposal methods no longer meet the current demands, making it urgent to evaluate and reuse the surface stability of these mined-out zones. Surface residual deformation is a process where voids and fissures within the mined-out zones are gradually filled and compacted, affecting the overlying rock structure. Additionally, groundwater significantly impacts the strength of the overlying rock, leading to increased subsidence. Therefore, predicting surface residual deformation while considering the effects of groundwater is crucial for forecasting surface deformation and assessing stability in mined-out zones. This study, taking into account the characteristics of subsidence zones and the impact of groundwater on the compaction of fractured rock masses, uses equivalent mining height and probability integral methods to develop a predictive model for surface residual deformation incorporating groundwater effects. Predictions for the study area show that groundwater exacerbates surface residual deformation, with various deformation values ranging from 33.8% to 51.9%. The surface stability categories are divided into stable and essentially stable regions based on the residual deformation’s impact on the working face. This model fully considers the influence of groundwater on residual deformation in mined-out zones, refining existing mining subsidence theories, addressing deformation issues caused by adverse groundwater factors, and providing a theoretical basis for predicting residual deformation and evaluating stability in mined-out zones, promoting the sustainable development of land and environmental resources in mining areas.

1. Introduction

Coal has long been a fundamental energy source and an important raw material in China. According to the 2023 “Annual Report on the Development of the Coal Industry”, China’s raw coal production reached 4.71 billion tons, with an average annual growth rate of 4.5%. Coal resources continue to play a major role in China’s energy structure. However, with prolonged extraction, coal resources in some regions have gradually depleted. As a result, many coal mines have been gradually closed, and the number of mines has decreased annually. In 2013, there were 12,000 coal mines, but by the end of 2023, the number had dropped to approximately 4300, a decrease of 64.2% over the past decade. As the mining area has expanded, many coal cities have developed extensive voids, with the area of subsidence continually increasing. With the closure of mines and the cessation of pumping, groundwater has gradually accumulated. These abandoned coal mines have led to water accumulation in the voids [1], which exacerbates the instability of the overburden and surface. This situation may reactivate previously stable fractured rock masses, leading to reactivation of the subsidence phenomenon. According to professional statistics, the surface subsidence for every ten thousand tons of coal (with subsidence greater than 10 mm as the standard) is 0.2 to 0.33 hm² [2]. Conservatively calculated with an average of 0.265 hm², and based on the total coal production data from 1949 to 2023, which amount to 104.04 billion tons, it is estimated that approximately 2.75706 million hm² of subsidence had formed by 2023. Assuming a coal production of 4 billion tons, the annual addition to the subsidence area is about 1.06 million hectares. The mining subsidence areas involve 4500 to 5000 km² of urban and rural construction land [3]. The development of mining cities and the scarcity of land resources have forced some buildings and transportation facilities to be located close to subsidence areas, increasing the urgent need for land reuse in these subsiding regions. The reuse of land resources in resource-based cities requires the residual deformation prediction and stability analysis of goaf areas. These efforts reduce the negative environmental impacts of mining activities and help prevent instability risks during land development in mining regions. As a result, they alleviate the pressure of land resource shortages during urban expansion and promote the sustainable development of land and environmental resources in mining areas. Therefore, residual surface deformation prediction and stability analysis of goaf areas under the influence of groundwater become particularly important.

Methods for calculating surface residual deformation in subsidence areas include equivalent mining prediction methods, limit settlement prediction methods, numerical simulation methods, nonlinear prediction methods, and mechanical prediction methods [4,5]. Guo [6] established a random medium model for predicting surface residual settlement in longwall mining subsidence areas. Han [7] proposed the principles and calculation methods for determining the surface subsidence coefficients in the central and boundary areas of the goaf, based on the analysis method for the stability of overlying rock in longwall goaf. Wang [8] divided the compaction characteristics of collapsed rock masses into zones using masonry beam theory and predicted residual settlement using the probability integral method. Guo [9] analyzed the spatial distribution and shape characteristics of overburden cavities, proposed a residual surface settlement prediction model, and validated it using numerical simulation experiments. Li [10] introduced the concept of the surface residual subsidence coefficient and a residual deformation prediction model based on the probability integral method. Yang [11] used the variable mining height probability integral method to predict the deformation, while Li [12] optimized the equivalent variable mining height probability integral method and proposed an improved version of the variable mining height probability integral method.

Currently, many scholars have conducted extensive research on residual deformation in subsidence areas and published various findings. However, there is limited research on the impact of groundwater on subsidence areas, and even less on predicting residual deformation caused by groundwater. How does groundwater affect the surface stability of subsidence areas? To what extent does groundwater influence residual surface deformation in these areas? Research on these questions is still quite limited, particularly regarding surface deformation prediction, stability assessment, and the management of subsidence areas. Therefore, studying the impact of groundwater on residual deformation in subsidence areas is crucial for evaluating surface stability, designing management plans for subsidence areas, and developing effective mitigation measures and disaster prevention strategies. This paper aims to address the issue of groundwater influence on residual deformation in goaf areas. Based on the deformation characteristics of overlying strata in goaf areas, a residual surface deformation prediction model, considering the effect of groundwater, was developed using the principles of equivalent mining height and the probability integral method. Additionally, the surface stability evaluation method was improved to promote the reuse of land resources and sustainable development in resource-based cities.

2. Analysis of the Mechanism of Surface Residual Deformation in Subsidence Areas Affected by Groundwater

After the extraction of underground coal seams, the initially balanced rock stress is disrupted, leading to fractures and collapses in the direct roof strata. In the collapsed area, the rocks become fragmented and structurally disordered. From the boundary to the center of the old subsidence area, the rock undergoes a process of transitioning from voids and partial compaction to complete compaction. Although surface displacement and deformation gradually stabilize, the residual voids and fractures in the upper rocks do not fully compact. The intrusion of groundwater compromises the relative stability of these rock layers. The interaction between groundwater and the damaged rock affects the stability of the rock layers, altering their mechanical properties and causing lubrication and softening effects. The lubricating effect of groundwater reduces the internal frictional resistance, decreases the friction angle, and exacerbates shear movement within the rock. As groundwater infiltrates the subsidence area and increases the rock’s moisture content, the rock may transition from a solid to a plastic state, or even become muddy. These significant changes in the mechanical properties lead to instability in the masonry beam structures, the closure of layered rock, fractures in the cantilevered areas, and the compaction of voids and under-compacted regions, resulting in new residual surface deformation [13].

Groundwater has a significant impact on rock strength, notably altering the physical state of the rock. The mechanical properties of rocks under saturated conditions differ markedly from those in dry conditions. Groundwater intrusion significantly reduces rock strength, and compresses and gradually stabilizes deformation. Consequently, groundwater intrusion destabilizes the originally stable state above the subsidence area, potentially worsening pre-existing instability. The research by Luo [14] indicates that water accumulation in subsidence areas is mainly found in fractured zones. During flow, infiltrated water not only removes soluble substances from the rock but also may cause the loss of fine particles, significantly reducing rock strength. The studies by Ma [15] and Yang [16] revealed that, under the same axial load, fractured rocks in a water-saturated state have a smaller residual expansion coefficient than in a dry state, leading to larger subsidence spaces that affect the surface, as shown in Figure 1. The alteration of rock properties due to groundwater affects the stability of the entire collapse zone. This impact needs to be thoroughly considered for more accurate surface deformation predictions and stability analysis.

Considering the differences in the morphology of old subsidence area overburden, these differences determine the distribution of “activation” spaces in various regions. The activation mechanisms in different areas of old subsidence zones vary significantly. The central region is mainly influenced by slow creep under sustained load. In contrast, the edge areas experience more sudden space release due to additional stress disturbances, the relaxation of structural stress, or reduced material strength, potentially leading to instability in key blocks of the masonry beam structures. The distribution of equivalent mining heights in different regions is as follows: In the center of the subsidence area, the equivalent mining height is equal to the residual space and usually remains constant. At the edges, influenced by masonry beam structures, the equivalent mining height varies with the rotation of the masonry beam structure blocks. Based on the characteristics of the collapse zone in subsidence areas, the subsidence area can be divided into two parts: the central compaction zone and the boundary under-compacted zone. These two parts have different overburden states, leading to different residual spaces and varying impacts on residual deformation. Additionally, due to the sequential support role of cantilevers within the subsidence area, there are rock layer fracture angles in the collapse zone. The height of the collapse zone changes, extending toward the center of the subsidence area. Due to the compaction of the overburden, it stabilizes at its highest point. Thus, the collapse zone can be divided into two parts: a height variation area on both sides of the collapse zone and a central stable height area.

3. Prediction Model for Surface Residual Settlement in Mines Considering the Impact of Groundwater

To provide a more detailed description of the deformation characteristics in the goaf, the space affecting residual surface deformation is divided into two parts based on the deformation sources. The first part is the residual space resulting from overburden deformation. The second part is the space influenced by groundwater in the goaf. This division allows for a better understanding of the residual deformation space in the collapse zone, leading to a more accurate prediction of residual surface deformation.

3.1. Residual Movements of Overburden Rock Formations Caused by Excavated Void Deformations

To estimate the maximum residual “activation” capacity of the compaction zone at the limit, it is assumed that the compaction zone is eventually 100% compacted. The maximum settlement value, including the residual settlement in single-layer mining subsidence, will not exceed the mining thickness $m_c$. Typically, the settlement coefficient is used to express the stability of coal seams. Let $q$ be the settlement coefficient at the initial stability; thus, the equivalent mining thickness for completed subsidence is $m_{c}q$.

Therefore, the equivalent mining thickness for residual settlement is calculated as follows:

$$m_1=m_c(1-q)$$

(1)

After the extraction of underground coal seams, the roof strata are left in a suspended state. Due to the inherent stiffness of the rock layers, fractures typically do not occur directly beneath the subsidence zone boundary. Instead, due to the cantilever effect, the fractures shift inward from the boundary. When lower rock layers form a cantilever structure, they provide some support to the upper rock layers, causing the cantilever length of the upper layers to exceed that of the lower layers and extend deeper into the subsidence zone. The compaction characteristics of the collapsed rock exhibit clear zoning from the edge of the old subsidence area to the center. A subsidence zone of length $D_1$ can be divided into an under-compacted zone of length $l_0$ and a central compaction zone. For the boundary of the subsidence area, considering the overburden morphology, the equivalent mining height m is variable. Here, it is simplified to a linear distribution [8].

Thus, the approximate equivalent mining height for the left under-compacted zone is as follows:

$$m_z=m_c-{\displaystyle \frac{m_{c}q}{l_0}}s$$

(2)

Similarly, the approximate equivalent mining height for the right under-compacted zone is as follows:

$$m_y={\displaystyle \frac{m_{c}q}{l_0}}(s-D_1)+m_c$$

(3)

3.2. Predictive Model for the Central Part of Old Mined-Out Areas Considering Groundwater Impact

For faulted rock boundaries, the area can be conceptualized as a semi-elliptical arch shape due to the sequential support effect of the cantilevered layers. In practical applications, a simplified trapezoidal model is used [17]. In this model, the bottom of the trapezoid represents the mining dimensions, with a strike length of $D_1$ and a width of $D_3$, as illustrated in Figure 1. Since the caving zone is located below the layer, a certain height of the trapezoid can be used to represent the caving zone. Thus, the height of the caving zone is equivalent to the height of the trapezoid. Consequently, the height of the caving rock layers exhibits a linear variation along each side of the trapezoid.

Based on the model of the old mining area, adjustments for the equivalent mining height must consider the impact of groundwater on the fractured rock mass. After the influence of groundwater, the swelling coefficient of the rock mass decreases. Therefore, the prediction model should account for the effect of groundwater on the fractured rock. The equivalent mining height needs to be adjusted for the groundwater impact in different areas. In the trapezoidal equilibrium model, the surface deformations in the central part are caused by the compaction of voids between fractured rocks within the caving zone. Combining this with the swelling coefficient allows for the derivation of the equivalent thickness of the compaction zone. The equivalent deformation of the trapezoid’s central part due to groundwater can be calculated using the residual swelling coefficients of the fractured rock before and after saturation, as shown in the following formula:

$${m_2=h}_k\left(K_g-K_s\right),$$

(4)

where $m_2$ is the equivalent mining thickness of the trapezoid’s central part, $K_g$ is the residual swelling coefficient in the natural state, $K_s$ is the residual swelling coefficient under water saturation, and $\left(K_g-K_s\right)$ represents the difference in the swelling coefficients of the rock mass due to groundwater impact. To predict the height of the caving zone, the calculation formula from “Regulations on Retaining Coal Pillars and Coal Mining in Buildings, Water Bodies, Railways, and Major Shafts” can be used to estimate the height of the caving zone.

Due to the variability of the overburden in the mining area, the swelling coefficient variation for the caving zone should be considered using a weighted average of the swelling coefficients of the fractured zone. The height changes of the caving rock mass before and after water saturation can be calculated as follows, representing the comprehensive variation in the caving zone:

$$m_1=\sum h_{ki}(K_{gi}-K_{si})$$

(5)

where $h_{ki}$ is the thickness of the i-th rock layer within the caving zone, $K_{gi}$ is the residual swelling coefficient of the i-th layer in its natural state, and $K_{si}$ is the residual swelling coefficient of the i-th layer after water saturation.

Simultaneously, the cantilever beam layers remain relatively intact and are less affected by groundwater. The primary impact is further compressed in the caving zone due to groundwater. The equivalent thickness of the trapezoid’s left side under groundwater influence can be approximated as follows:

$$m_{zw}={\displaystyle \frac{\sum h_{ki}\left(K_{gi}-K_{si}\right)}{\sum h_{ki}}}s·\mathrm{t}\mathrm{a}\mathrm{n}\beta$$

(6)

Similarly, the approximate thickness of the under-compressed zone on the right side under groundwater influence is as follows:

$$m_{yw}={\displaystyle \frac{\sum h_{ki}\left(K_{gi}-K_{si}\right)}{\sum h_{ki}}}\left(D_1-s\right)\mathrm{t}\mathrm{a}\mathrm{n}\beta ,$$

(7)

where β is the angle of the rock layer’s fault.

3.3. Analysis of Dimensions for Deformation Caused by Mining-Induced Movement

The overburden of old mining areas, due to the masonry beam structure, results in fractured rock masses with varying degrees of compaction. This leads to the formation of under-compacted and compacted zones. From the perspective of overburden structures, the boundary between compacted and under-compacted zones can be defined as the point where the masonry beam structure blocks rotate to a horizontal position or where the rate of change of the displacement curve of the masonry beam structure is zero. From the viewpoint of mining pressure, the boundary between the under-compacted and compacted zones approximates the position where the original rock stress is restored behind the caving rock mass at the mining boundary. Assuming that the position where the original rock stress is restored on the collapsed rock mass is at a distance $l_0$ from the boundary of the goaf, the size of the under-compacted zone extends $l_0$ from the cut, the stop line, and both upper and lower drifts towards the interior of the goaf. In the strike direction, the range of the compacted zone is $(D_1-2l_0$).

In practice, the value of $l_0$ varies across different mining areas and depends on factors such as the overburden structure, mining thickness, mining depth, mining method, and roof management methods. A reliable method to obtain $l_0$ is to measure the rock layer movement curves or the mining pressure behind the mining area. For old mining areas, obtaining $l_0$ through measurements is challenging. When using the probabilistic integral method for subsidence calculations, the computed value is the actual mining value minus the inflection point offset $l_0$. The inflection point offset reflects the voids formed by the cantilever beam structure of the coal seam roof, and thus, $l_0$ can be used to replace the length of the boundary void region.

From the above analysis, the length of the intermediate compacted zone is the mining length minus twice the inflection point offset, i.e., $l_1=D_1-s_3-s_4$, where l₁ is the length of the intermediate compacted zone, and $s_3$ and $s_4$ are inflection point offsets used to replace the strike boundary region $l_0$. For the dip direction, the width is $D_3$, with l₂ representing the width of the intermediate compacted zone, and $s_1$ and $s_2$ are inflection point offsets used to replace the dip boundary region $l_0^{\prime }$. For the trapezoidal equilibrium model of the fractured rock mass in the caving zone, the dimensions can be determined by the angle of rock layer faulting. Typically, the fault angle β is between 55° and 60° [17], resulting in a width from the mining boundary of $h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta$. It is important to note that the calculation method for the inflection point offset differs from that for the trapezoidal slope, so the width is generally different.

3.4. Predictive Models for Old Mining Areas Considering Groundwater Impact

3.4.1. Predictive Model for Residual Deformation of Main Cross-Section Considering Groundwater Impact

Based on the classification above and the influence of groundwater on rock masses, the mining face is divided into two parts: the boundary under-compacted zone and the central compacted zone. Each part is activated and calculated separately. The caving zone is also divided into two parts. This leads to the development of a predictive model for surface residual movement deformation in mining areas, considering groundwater impact. In the coal seam mining process, as shown in Figure 2, the actual mining length of the coal seam is af. Due to the influence of the roof layer above the coal seam, the surface subsidence location shifts, resulting in lateral inflection point offsets $s_3$ and $s_4$. The length of the compacted zone is be, and the lengths of the caving zone slopes are ac and df.

In the probability integral method, if the horizontal coordinate of a mining unit is s and the horizontal coordinate of any point A on the surface is x, the residual subsidence mathematical model for each part of the old mining area caused by equivalent mining thickness is calculated as follows:

$$w(x)=\mathrm{c}\mathrm{o}\mathrm{s}\alpha \int m{\displaystyle \frac{1}{r}}{e}^{-\pi {\displaystyle \frac{{(x-s)}^2}{r^2}}}ds,$$

(8)

where $w\left(x\right)$ represents the surface subsidence value for the mining unit; m denotes the equivalent mining height of each area; and r is the primary influence radius.

Considering the boundary under-compacted zone and the compacted zone, and the impact of groundwater, the equivalent mining heights for different regions are adjusted. The cumulative calculation for each region yields the following:

$$w(x)={\int }_{ab}{m}_3\mathrm{c}\mathrm{o}\mathrm{s}\alpha w(x)ds+{\int }_{bc} m_2\mathrm{c}\mathrm{o}\mathrm{s}\alpha w(x)ds+{\int }_{cd} m_1\mathrm{c}\mathrm{o}\mathrm{s}\alpha w\left(x\right)ds+{\int }_{de} m_2^{\prime }\cos \alpha w\left(x\right)ds+{\int }_{ef} m_3^{\prime }\mathrm{c}\mathrm{o}\mathrm{s}\alpha w\left(x\right)ds,$$

(9)

where $m_3$, $m_2$, $m_1$, $m_2^{\prime }$, and $m_3^{\prime }$ are the equivalent mining heights for the regions ab, bc, cd, de, and ef, respectively.

The corresponding formulas for calculating the inclination, curvature, horizontal movement, and horizontal deformation of each impact area are as follows:

$$i(x)={\displaystyle \frac{dw(x)}{dx}}=-\mathrm{c}\mathrm{o}\mathrm{s}\alpha \int m{\displaystyle \frac{2\pi }{r^3}}(x-s)e^{-\pi {\displaystyle \frac{{(x-s)}^2}{r^2}}}ds$$

(10)

$$k(x)={\displaystyle \frac{di\left(x\right)}{dx}}=-\mathrm{c}\mathrm{o}\mathrm{s}\alpha \int m{\displaystyle \frac{2\pi }{r^3}}\left(1-{\displaystyle \frac{2\pi {(x-s)}^2}{r^2}}\right)e^{-\pi {\displaystyle \frac{{(x-s)}^2}{r^2}}}ds$$

(11)

$$u\left(x\right)=bri\left(x\right)=-\mathrm{c}\mathrm{o}\mathrm{s}\alpha \int m{\displaystyle \frac{2\pi b}{r^2}}(x-s)e^{-\pi {\displaystyle \frac{{(x-s)}^2}{r^2}}}ds$$

(12)

$$\epsilon (x)=brk(x)=-\mathrm{c}\mathrm{o}\mathrm{s}\alpha \int m{\displaystyle \frac{2\pi b}{r^2}}\left(1-{\displaystyle \frac{2\pi {(x-s)}^2}{r^2}}\right)e^{-\pi {\displaystyle \frac{{(x-s)}^2}{r^2}}}ds$$

(13)

According to Formulas (8) and (10)–(13), cumulative calculations for each impact area and factor are performed to compute the residual movement deformation values at any point based on the different residual deformations in the regions and factors.

3.4.2. Prediction Model for Residual Deformation at Any Point Considering the Impact of Groundwater

In practical applications, surface subsidence often occurs, not only at the main section, but also at any point within the surface subsidence basin, affecting movement and deformation in various directions. According to the aforementioned classification and considering the impact of groundwater on the rock mass, the mining face is divided into two plane regions: the boundary under-compacted zone and the central compacted zone. Additionally, the subsidence zone is divided into two plane regions for activation calculations. This approach allows for the development of a residual movement and deformation prediction model for the mining area surface, taking groundwater impact into account.

In the coal mining process, as illustrated in Figure 3, the action of the overlying rock layer’s roof causes shifts in the surface subsidence location, leading to strike offset distances $s_3$ and $s_4$, and dip offset distances $s_1$ and $s_2$. The under-compacted zone is labeled as D_A, and the compacted zone is represented as D_B + D_C. The red area in the figure indicates the initially stable subsidence space, the yellow area denotes the residual deformation space, and the blue area represents the residual deformation space affected by groundwater.

According to the probability integral method, the subsidence model for any point (x, y) on the surface caused by the unit mining (s, t) in different regions is calculated as follows:

$$w\left(x,y\right)=\underset{D}{{\displaystyle \iint }} m\mathrm{c}\mathrm{o}\mathrm{s}\alpha {\displaystyle \frac{1}{r^2}}{\mathrm{e}}^{-\pi {\displaystyle \frac{(x-s{)}^2+(y-t{)}^2}{r^2}}}dsdt,$$

(14)

where D denotes the influence region.

The subsidence caused by different regions and factors in old mined-out areas is as follows:

$$\begin{gathered} w\left(x,y\right)=\underset{D_A}{{\displaystyle \iint }} m_{D_A}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi \left[\right(x-s{)}^2+(y-t{)}^2]}{r^2}}}dsdt+\underset{D_B}{{\displaystyle \iint }} m_{D_B}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi \left[\right(x-s{)}^2+(y-t{)}^2]}{r^2}}}dsdt \\ +\underset{D_c}{{\displaystyle \iint }} m_{D_c}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi \left[\right(x-s{)}^2+(y-t{)}^2]}{r^2}}}dsdt \end{gathered}$$

(15)

Since the boundaries between the under-pressure and over-pressure zones are defined based on inflection point offset distances, while the trapezoidal central region and slopes are defined based on the rock layer rupture angle and collapse belt height, these boundaries do not necessarily coincide. Therefore, they are delineated and calculated separately. According to this method, the mining area is divided into the under-pressure zone D_A and the over-pressure zones D_B + D_C. The equivalent mining height for the over-pressure zone D_A is calculated using Formula (1). The under-pressure zone D_A is further divided into four regions: left, right, upper, and lower. Each region is integrated separately. Here, the boundaries of these regions are determined by connecting the corner points of the trapezoid. The equivalent mining heights for the left and right regions are calculated using Formulas (2) and (3), while the equivalent mining height formulas for the upper and lower regions are as follows:

$$m_s={\displaystyle \frac{m_{c}q}{s_1}}\left(t-D_3\right)+m_c$$

(16)

$$m_x=m_c-{\displaystyle \frac{m_{c}q}{s_2}}t$$

(17)

Considering the presence of the rock layer rupture angle, and based on the aforementioned trapezoidal collapse belt equilibrium model, the impact range of groundwater is divided into two parts: the trapezoidal central region D_C and the trapezoidal slopes D_A + D_B. The equivalent mining height for the trapezoidal central region D_C is calculated using Formula (5). For the trapezoidal slopes D_A + D_B, the equivalent mining height is divided into four regions. The boundary lines of these regions are defined by connecting the corner points of the trapezoid. The equivalent mining height for the left and right regions is calculated using Formulas (6) and (7), while the equivalent mining height formulas for the upper and lower regions are as follows:

$$m_{sw}={\displaystyle \frac{\sum h_{ki}\left(K_{gi}-K_{si}\right)}{\sum h_{ki}}}\left(D_3-t\right)\tan \beta$$

(18)

$$m_{xw}={\displaystyle \frac{\sum h_{ki}\left(K_{gi}-K_{si}\right)}{\sum h_{ki}}}t·\mathrm{t}\mathrm{a}\mathrm{n}\beta$$

(19)

Due to the differing regions and factors causing residual deformation subsidence, and since each region’s equivalent mining height varies, individual integration for each source region is required to compute the residual deformation displacement value for any point. The subsidence, based on Formulas (1)–(3), (5), (6) and (15)–(19), is calculated as follows:

$$\begin{array}{cc}\hfill w\left(x,y\right)={{\displaystyle \int }}_0^{s_3} & ({{\displaystyle \int }}_{{\displaystyle \frac{s_2}{s_3}}\cdot s}^{-{\displaystyle \frac{s_1}{s_3}}\cdot s+D_3} m_z\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \\ & +{{\displaystyle \int }}_{D_1-s_4}^{D_1} ({{\displaystyle \int }}_{-{\displaystyle \frac{s_2}{s_4}}\cdot s+{\displaystyle \frac{s_2}{s_4}}\cdot D_1}^{{\displaystyle \frac{s_1}{s_4}}\cdot s+\left(D_3-{\displaystyle \frac{s_1}{s_4}}\cdot D_1\right)} m_y\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \\ & +{{\displaystyle \int }}_{D_3-s_1}^{D_3} ({{\displaystyle \int }}_{{\displaystyle \frac{s_3}{s_1}}\cdot \left(D_3-t\right)}^{{\displaystyle \frac{s_4}{s_1}}\cdot t+D_1-{\displaystyle \frac{s_4}{s_1}}\cdot D_3} m_s\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}ds)dt\hfill \\ & +{{\displaystyle \int }}_0^{s_2} ({{\displaystyle \int }}_{{\displaystyle \frac{s_3}{s_2}}\cdot t}^{D_1-{\displaystyle \frac{s_4}{s_2}}\cdot t} m_x\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}ds)dt\hfill \\ & +{{\displaystyle \int }}_{s_3}^{D_1-s_4} ({{\displaystyle \int }}_{s_2}^{D_3-s_1} m_1\mathrm{c}\mathrm{o}\mathrm{s}\alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \\ & +{{\displaystyle \int }}_0^{h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta } ({{\displaystyle \int }}_s^{-s+(D_3)} m_{zw}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \\ & +{{\displaystyle \int }}_{D_1-h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta }^{D_1} ({{\displaystyle \int }}_{-s+D_1}^{s-D_1+D_3} m_{yw}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \\ & +{{\displaystyle \int }}_{D_3-h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta }^{D_3} ({{\displaystyle \int }}_{-t+(D_3)}^{t-(-D_1+D_3)} m_{sw}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}ds)dt\hfill \\ & +{{\displaystyle \int }}_0^{h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta } ({{\displaystyle \int }}_t^{D_1-t} m_{xw}\cos \alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}ds)dt\hfill \\ & +{{\displaystyle \int }}_{h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta }^{D_1-h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta } ({{\displaystyle \int }}_{h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta }^{D_3-h_k/\mathrm{t}\mathrm{a}\mathrm{n}\beta } m_2\mathrm{c}\mathrm{o}\mathrm{s}\alpha {\displaystyle \frac{1}{r^2}}{e}^{-{\displaystyle \frac{\pi [(x-s{)}^2+(y-t{)}^2]}{r^2}}}dt)ds\hfill \end{array}$$

(20)

Based on Formula (14), the inclination $i(x,y,\phi )$ along the φ direction is calculated as follows:

$$\begin{gathered} i\left(x,y,\phi \right)={\displaystyle \frac{\partial w(x,y)}{\partial x}}\mathrm{c}\mathrm{o}\mathrm{s}\phi +{\displaystyle \frac{\partial w(x,y)}{\partial y}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \\ =m\underset{D}{{\displaystyle \iint }} (-\cos \left(\alpha \right){\displaystyle \frac{2\pi }{r^4}}{\mathrm{e}}^{-\pi {\displaystyle \frac{(x-s{)}^2+(y-t{)}^2}{r^2}}}\left[\cos \left(\phi \right)\left(x-s\right)+\sin \left(\phi \right)\left(y-t\right)\right])dsdt \end{gathered}$$

(21)

The curvature $k(x,y,\phi )$ along the φ direction is calculated as follows:

$$\begin{array}{cc}\hfill k\left(x,y,\phi \right)={\displaystyle \frac{\partial i\left(x,y,\phi \right)}{\partial x}}& \mathrm{c}\mathrm{o}\mathrm{s}\phi +{\displaystyle \frac{\partial i\left(x,y,\phi \right)}{\partial y}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \hfill \\ & =m\underset{D}{{\displaystyle \iint }} \cos \left(\alpha \right){\displaystyle \frac{2\pi }{r^4}}{\mathrm{e}}^{-\pi {\displaystyle \frac{(x-s{)}^2+(y-t{)}^2}{r^2}}}[\mathrm{c}\mathrm{o}\mathrm{s}\phi (-\cos \phi +{\displaystyle \frac{2\pi (x-s{)}^2}{r^2}}\cos \phi \hfill \\ & +{\displaystyle \frac{2\pi \left(x-s\right)\left(y-t\right)}{r^2}}\sin \phi )+\mathrm{s}\mathrm{i}\mathrm{n}\phi (-\sin \phi +{\displaystyle \frac{2\pi (y-t{)}^2}{r^2}}\sin \phi \hfill \\ & +{\displaystyle \frac{2\pi \left(x-s\right)\left(y-t\right)}{r^2}}\cos \phi )]dsdt\hfill \end{array}$$

(22)

The horizontal displacement $u\left(x,y,\phi \right)$ along the φ direction is as follows:

$$u\left(x,y,\phi \right)=bri\left(x,y,\phi \right)$$

(23)

The horizontal deformation $\epsilon \left(x,y,\phi \right)$ along the φ direction is as follows:

$$\epsilon (x,y,\phi )=brk(x,y,\phi )$$

(24)

Using Formulas (15) and (21) through (24), cumulative calculations are performed for different influencing factors in each region. This allows for the calculation of residual displacement deformation values at any point based on the residual deformation in different regions and caused by different factors.

In conclusion, by using the parameters from the probability integral method and combining them with the empirical values of the fracture angle, the model can analyze the residual deformation space of the overburden and the space influenced by groundwater. This allows for residual deformation prediction considering the impact of groundwater. Subsequently, a surface stability analysis can be performed in accordance with the relevant standards to evaluate surface stability.

3.5. Sensitivity Analysis of Residual Surface Movement Deformation Prediction Parameters Considering the Impact of Groundwater

The factors influencing residual deformation under the effect of groundwater are divided into two parts: one is the residual deformation space in the goaf without groundwater influence, and the other is the space influenced by groundwater. Based on the model, the key parameters affecting the maximum residual subsidence include the residual subsidence coefficient, coal seam dip angle, mining thickness, height of the collapse zone, and the difference in the swelling factor between dry and wet collapse zones. The parameters are set as follows: the residual subsidence coefficient $q^{\prime }$ ranges from 0.1 to 0.3, the mining thickness m ranges from 1 to 3 m, the height of the collapse zone $h_k$ ranges from 2.6 to 13.5 m, the difference in swelling factor between dry and wet collapse zones $(K_{gi}-K_{si})$ ranges from 0.01 to 0.04, and the coal seam dip angle α ranges from 0° to 10°. The Sobol sensitivity analysis method was employed to analyze the parameter sensitivity of maximum subsidence. The first-order Sobol indices from the sensitivity analysis are presented in

Table 1. First-order Sobol indices for each parameter.

Parameter	First-Order Sobol Index	Confidence Interval
$q^{\prime }$	0.3367	0.0245
$m$	0.337	0.0238
$h_k$	0.1563	0.0164
$K_{gi}-K_{si}$	0.1229	0.0165
$\alpha$	0.0002	0.0006

.

The first-order Sobol indices for $q^{\prime }$ and m are both 33.7%, indicating that they have a very high individual contribution to the results and represent the dominant sensitivity factors. Together, they account for the majority of the sensitivity contribution. The contribution of $h_k$ is approximately 15.6%, suggesting that while its influence is secondary to $q^{\prime }$ and m, it is still a significant parameter. The contribution of $(K_{gi}-K_{si})$ is 12.3%, indicating a notable impact, though smaller than that of $h_k$. The influence of α is minimal, with a contribution of only 0.0002%, making it an insignificant factor in the results. The confidence intervals for all parameters are relatively small, indicating that the results are reliable.

4. Case Study

4.1. Study Area Overview

The study area is located in the Xuzhou mining district in Jiangsu Province, China, as shown in Figure 4. The mining area is part of the Huang-Huai alluvial plain, characterized by flat terrain, with surface elevations ranging from +32 to +36 m. The area slopes from southwest to northeast, with a gradient of approximately 0.001. In the southeastern part of the mining district, small hills are formed by Cambrian and Ordovician strata. To the northeast, the district borders the Weishan Lake, and the Beijing–Hangzhou Grand Canal and the Shun Ti River traverse the northern edge of the district diagonally. To the south of the mining area is the Taoyuan River. The region has a temperate climate with abundant annual rainfall. This study focused on working face 2114 to estimate the residual deformation pattern influenced by groundwater. The mining method used was subsidence mining. The working face is approximately rectangular, with a mining length of 422 m, a width of 98 m, an average depth of 139.5 m, and an average thickness of 1.2 m. The coal seam has a relatively small dip angle of 14°.

4.2. Estimation of Groundwater Impact on Surface Residual Deformation Prediction Analysis

Based on the mining geological conditions of the study area, this section analyzes the residual deformation patterns of the mined-out area under the influence of groundwater. Considering the geological conditions of the mining district [18] and typical geological profiles [19,20,21], the subsidence coefficient in the central region of the old mined-out area is 0.9. The parameters for the predictive model, taking into account the influence of groundwater, are outlined in

Table 2. Predictive model parameters.

Parameter	Goaf Zone Height ${\mathit{h}}_{\mathit{k}}$/m	Main Effect on Tangent Value	Inflection Point Offset Distance ${\mathit{s}}_1$/m	Inflection Point Offset Distance ${\mathit{s}}_2$/m	Inflection Point Offset Distance ${\mathit{s}}_3$/m	Inflection Point Offset Distance ${\mathit{s}}_4$/m	Horizontal Movement Coefficient	Difference in the Fragmentation and Swelling Coefficients between Dry and Wet Caving Zones
Value	4.87	1.8	7.6	6.35	6.975	6.975	0.22	0.0345

.

Using Python 3.10 programming, we calculated the site’s deformation based on the predictive model. For ease of display and analysis, we set the mining plane along the strike as the horizontal axis and the dip as the vertical axis. The residual subsidence without groundwater influence and the residual subsidence, residual tilt, residual curvature, residual horizontal deformation, and residual horizontal movement considering groundwater influence are shown in Figure 5. The black rectangle indicates the boundary of the working face.

In Figure 5a,b, it can be observed that the subsidence basin is relatively small due to the low coal seam inclination. The residual deformation considering groundwater impacts shows significant changes compared to when groundwater is not considered, with increased deformation. The calculations revealed that groundwater influences generally lead to increased amplitudes in various curves. Without groundwater influence, the maximum residual subsidence is located near the inflection point inside the mining area’s cutting face and stop line, not at the center of the mining area. For old mining areas, the cavity is larger near the cutting face and stop line due to cantilever beam support, and deformation progressively reflects on the surface as compaction occurs and groundwater effects are considered. With groundwater impact, the bottom of the subsidence basin is nearly horizontal, due to changes in the height of the collapse zone edges and reduced active space compared to the center.

Figure 5b–j indicate that groundwater influence causes the re-compaction and deformation of fractured rock in the mining area, increasing subsidence in the center and decreasing it toward the sides. The main increase in residual subsidence is in the central area of the mining area, with the basin’s center remaining roughly horizontal. The maximum residual subsidence increases from 138.522 mm to 280.356 mm, with the largest values near the inflection points on the sides. The patterns of residual tilt along the strike and dip are similar, with varying degrees of change by region. The absolute tilt values are highest outside the mining plane, decreasing inward and approaching zero, with a maximum absolute residual tilt of 3.879 mm/m. Residual horizontal movement follows a similar pattern, with the maximum value of 66.134 mm occurring outside the mining area under groundwater influence, decreasing inward and approaching zero. The maximum absolute value of residual curvature is 0.084 mm/m², with positive values outside the mining boundary and negative values inside. The pattern of residual horizontal deformation is similar to residual curvature, with maximum negative values of 1.429 mm/m and positive values outside the mining boundary, and negative values inside.

Based on the principles outlined above, further analysis was conducted on the study area. Surface residual deformations were calculated for both scenarios: with and without considering groundwater effects. The results are summarized in Figure 6, which shows the proportions of each part. The maximum residual subsidence increased from 134.74 mm to 280.36 mm, representing a 108.07% increase. The maximum residual tilt grew from 2.094 mm/m to 3.879 mm/m, an increase of 85.24%. The maximum residual curvature rose from 0.056 mm/m² to 0.084 mm/m², marking a 50% increase. The maximum residual horizontal deformation expanded from 0.946 mm/m to 1.429 mm/m, reflecting a 51.06% increase. The maximum residual horizontal movement increased from 35.706 mm to 66.134 mm, showing an 85.22% rise. It is evident that the impact of groundwater shows a significant increase. This is related to the subsidence coefficient: a larger subsidence coefficient results in a smaller residual subsidence coefficient. Consequently, the relative increase in residual subsidence due to groundwater impact is more pronounced.

Figure 6 indicates that, when considering groundwater effects, the extent of surface deformation in the subsidence zone varies significantly. The proportions of each residual deformation zone under groundwater influence range from 33.3% to 51.9%. Therefore, groundwater effects must be predicted and considered when evaluating the stability of subsidence zones.

For assessing the stability of the surface in subsidence zones, reference should be made to the stability classification criteria outlined in the relevant standard [22], as shown in

Table 3. Assessing the stability of mining-induced foundations based on surface (residual) deformation values.

Stability Level	Horizontal Deformation ε (mm/m)	Curvature K (10−3/m)	Tilt i (mm/m)	Criteria
Stable	ε ≤ 2.0	K ≤ 0.2	i ≤ 3.0	All conditions met
Basically Stable	2.0 < ε ≤ 4.0	0.2 < K ≤ 0.4	3.0 < i ≤ 6.0	Any one condition met
Unstable	ε > 4.0	K > 0.4	i > 6.0	Any one condition met

.

Based on Figure 5, it can be concluded that the strike direction has been fully mined, while the dip direction remains partially mined. As a result, horizontal deformation and curvature along the strike exhibit smaller deformation in the central area of the basin, with larger deformation near the cutting face and stop line. However, for horizontal deformation and curvature along the dip direction, since the dip has not been fully mined, the central basin area still shows relatively large deformation. Tilt deformation is generally smaller in the central basin area along both the strike and dip directions. The calculations indicate that when both the strike and dip directions are fully mined, the values of horizontal deformation, curvature, and tilt in the central basin area are relatively small.

The study area was subdivided based on Table 2, and the calculation results are shown in Figure 7.

In Figure 7, it can be observed that the surface residual deformation caused by this working face is generally in a stable state, with only parts near the mining boundary being in a basically stable condition.

5. Conclusions

This study analyzed the equivalent mining height in the collapse zone of a mined-out area under the influence of groundwater. The analysis was based on the distribution characteristics of fractured rock, the principles of the equivalent mining height and the probability integral method, and the variation patterns of the residual swelling coefficients of fractured rock. A groundwater-influenced residual deformation prediction method, based on the trapezoidal distribution of the collapse zone, was established. This method improves the surface deformation prediction model for mined-out areas affected by groundwater. The specific conclusions are as follows:

(1) The predictive model fully accounts for the changes in the residual swelling coefficients of rocks in different regions of the collapse zone, as well as the deformation characteristics of the overburden in the mined-out area. It corrects surface deformation caused by the varying impacts of groundwater and the original subsidence space in the mined-out area. By thoroughly considering the original deformation space and the influence of groundwater, the prediction results are more scientifically sound and reasonable, providing a valuable reference for analyzing the surface stability of old mined-out areas after flooding. The research findings contribute to the reuse of land resources and the sustainable development of resource-based cities;

(2) Applying the proposed model to the study area yielded results consistent with the characteristic slight “W” shape of the residual subsidence basin on the surface of the old mined-out area. Under the influence of groundwater, the residual deformation was greater than the original residual deformation, with the proportions of various residual deformation values in the study area ranging from 33.8% to 51.9%. Overall, the residual deformation was relatively minor. The presence and proportion of residual deformation caused by groundwater influence in the mined-out area should be fully considered. The residual deformation of the working face affects the surface stability, resulting in stable and basically stable areas on the surface;

(3) By analyzing the contour lines of residual horizontal deformation ϵ, curvature k, and tilt i in the study area, and evaluating the stability of the study area according to the stability criteria for foundations in mined-out areas, it was determined that the residual deformation caused by this working face is relatively small. A small portion of the area is in a basically stable state, while the majority is in a stable state. If the surface is affected by multiple mining faces, a cumulative analysis should be conducted.