# Dynamic Prediction Method of 3D Spatial Information of Coal Mining Subsidence Water Area Integrated with Landsat Remote Sensing and Knothe Time Function

Authors: Hui Liu; Yu Li
Journal: Geofluids
 (2022)

---

## Abstract

Subsidence water disaster has become a major problem to the ecological environment because of the submerged villages, farmland destruction, and ecological destructive changes in the mining subsidence area of high-water level. Taking the 1031 working face of Wugou Coal Mine in Huaibei, Anhui, China, as the research subject, (1) a three-dimensional (3D) spatial information dynamic prediction method was proposed for high-water-level coal mining subsidence areas by combining the Knothe time function based on the probability integration method (PIM) and the principle of water balance. (2) The dynamic evolution law of the water accumulation area in the high-water-level coal mining subsidence area was studied. (3) The applicability of the dynamic prediction model of the water accumulation range in the high-water-level coal mining subsidence area was verified. The results showed that the dynamic prediction of the 3D spatial information of the high diving area is highly accurate and can be suitable for the dynamic prediction by comparison with the results of remote sensing monitoring and field measurement. As a result, technical reference and theoretical basis for the comprehensive assessment and remediation of the ecological environment of the high diving mining subsidence area were proposed.


---

## Body

## 1. Introduction

Coal is China’s main energy source, accounting for approximately 70% of primary energy production and consumption [1,2]. It has an important position in the economy [3] and is an important factor that drives world economic development. With the rapid development of contemporary socioeconomic and industrial production, coal consumption as a conventional fuel is continuously growing [4]. According to the World Energy Yearbook, China’s coal resource output reached 3846 billion tonnes in 2019, accounting for 47.3% of the global total output, remarkably higher than the world’s average output level [5]. However, coal mining not only caused serious safety accidents [6–8] but also caused a range of profound environmental problems, such as decreased soil quality, ecosystem degradation, biodiversity loss, and landscape destruction [9]. Fractured rocks in the subsurface are ubiquitous, and the dynamics of mass transfer in the fractured rocks play an important role in understanding the problem in engineering geology and environmental geology [10]. Amongst them, underground coal mining directly causes large surface subsidence [11–13]. Subsidence and accumulated water disasters have become the primary ecological environment problems faced by coal mine subsidence areas with high phreatic levels. Village flooding, farmland damage, and disruptive changes in the ecosystem have brought great challenges to the construction of ecological civilization in the mining area [14,15], and these are the urgent problems in the coal mining areas that need to be addressed.The negative impact of coal mining in the areas with high phreatic water levels in eastern China is particularly serious. This area has the characteristics of high phreatic level and thick loose layer coverage [16]. The mining area is an area with coal resources. Diving is gravity water with a free surface above the first continuous aquifer below the surface. Its upper part has no continuous and complete waterproof top plate, and the diving water surface is a free water surface, called a diving surface. Any point on the diving surface is called the diving level according to the absolute height of the reference level. A high diving level, as the name suggests, is a higher level of diving underground. High-water-level mining areas are mainly distributed in the eastern part of China. High underground water level causes the surface to sink, and water is accumulated after coal mining, thus forming a unique ecological environment of the mining area with a combination of water and land.Understanding in situ stress is the foundation of studying deep rock mechanics [17]. Generally, deformation, fracture propagation, and energy release are highly associated with mining-induced stress evolution [18]. When underground coal is extracted, the surrounding rock loses the original stress balance [19,20], causing stress redistribution. With the increasing maximum stress, the average stress applied and the duration of rock suffering from a high-stress increase. Therefore, reaching the strength limit becomes easier, and the fatigue life is shorter [21]. When the concentrated stress of the roadway roof exceeds the strength limit, roof layer deformation, fracture, and collapse [22] cause pressure subsidence and water consolidation subsidence. These conditions eventually cause the movement and deformation on the surface subsidence basin and produce a large number of cracks in high diving areas, where groundwater formation easily connects to the subsidence area. The water in the subsidence area mainly comes from atmospheric precipitation, surface runoff, surface water evaporation, and groundwater seepage. When the water input and output of various factors reach equilibrium, the range of water in the subsidence area remains unchanged. Figure 1 shows the image of the surface subsidence caused by coal mining in a typical high-water-level mining area of Lianghuai coal mines in China. The images on the left and right are the Huainan Guqiao Coal Mine and Huaibei Coal Mine, respectively. The subsidence water disaster has become a primary ecological environment problem because of the high diving mining subsidence areas, submerged villages, farmland damage, and disruptive changes to the ecosystem, causing great challenges to the sustainable development of the mining area. In the past 25 years, the water accumulation area in the coal mining subsidence area of Anhui Province has increased by approximately six times, from 18.95 km2 to 118.09 km2, with an annual growth of 3.97 km2. In addition, some water areas account for 30% to 50% of the total coal mining subsidence area [23]. The large subsidence area gradually evolves from an original land ecosystem into a water land composite ecosystem, causing irreversible damage to the original local ecosystem balance. On the other hand, surface subsidence destroys large areas of cultivated land and prominent human and land contradictions, thereby seriously affecting the local economic development, social structure, and stability.Figure 1
High-water-level coal mining subsidence area.Surface subsidence caused by coal mining directly damages the ecological environment in the mining area. Especially in the areas with high diving sites, coal resource mining not only causes surface subsidence but also accumulates a large range of water in the subsidence area. As a result, the ecological environment is seriously damaged. In the view of the surface subsidence water in high diving mining subsidence area and a series of environmental problems, field mapping, mining subsidence estimation, and remote sensing monitoring have been widely used to obtain the geospatial information data of the coal mining subsidence area.(1)
Field surveying and mapping: the surface mobile observation station is set, a large number of measured data are collected, theoretical analysis is performed, the surface mobile data of different mining areas are obtained, and the surface deformation laws under different geological mining conditions are analyzed. The monitoring of water mainly uses the entire station instrument, theodolite, GPS, and other traditional measurement instruments. In accordance with the water characteristics, the characteristic point coordinates are measured, and then water statistics are obtained according to the internal industry treatment.However, the surface form and water situation in the subsidence basin are dynamically affected by many factors. This method is suitable for small water areas with limited applicability.(2)
Mining subsidence is expected to happen because of the construction of underground mines. In recent years, with the further deepening of mining under different conditions, scholars have developed and established many new mining subsidence prediction methods, such as the Knothe surface subsidence prediction theory based on the traditional influence function method [24] and the mining subsidence prediction model based on spatial statistical methods. From the mechanical perspective, the linear elastic analysis principle and the surface element principle are introduced. Different geological conditions and mining methods have led to the development of the prediction theories and models of inclined coal seam surface subsidence prediction, surface settlement prediction of abandoned column old mining area, surface subsidence prediction of the shallow mining area, and surface subsidence prediction of chamber column mining.However, in a special geological situation, the surface subsidence characteristics are different from the conventional mining subsidence characteristics. Thus, the conventional subsidence prediction model cannot fully apply to the surface subsidence prediction and needs further verification of the accuracy of the surface subsidence prediction model of coal seam mining.(3)
Remote sensing monitoring: in terms of water evolution, Suping et al. considered the Huainan mining area and used transmission electron microscopy (TM) images to extract the information of water at different times using the method of principal component analysis and analyzed the changes in the water area [25]. With the development of remote sensing technology, domestic high-resolution satellite images have been used for monitoring the range of stagnant water and scope information extraction.However, the surface form and water in the mining subsidence area are dynamically affected by many factors. The observation data obtained by the existing methods lack timeliness, thus making the timely prediction of the 3D spatial information of water in the surface subsidence area difficult and affecting the accuracy of the planning and design of the coal mining subsidence area.Studies on the field mapping of the high mining subsidence, mining subsidence dynamic forecast, and remote sensing monitoring are abundant. However, the surface subsidence water and many factors, such as comprehensive mining subsidence, surface water range evolution, and water circulation, require further investigation. In addition, high mining subsidence and water range dynamic evolution law should be analyzed. (1) The dynamic evolution law of water in the subsidence area of the high diving level is still unclear. Restricted by the means of field mapping, the existing studies mainly targeted the range change of surface water before and after subsidence. However, studies on the production, increase, stability, and dynamic laws of surface water at different mining stages are limited. (2) The accuracy of the mining subsidence expected model should be verified. Further validation is expected for surface subsidence with special geological conditions. (3) The dynamic prediction of surface water range in high diving coal mining subsidence area is limited. Although the development of remote sensing technology has been relatively mature, most existing studies focus on a single mining area or a working surface, few large areas, and long-time effect of water dynamic evolution characteristics. In addition, the impact of the dynamic change of water analysis is less, the subsidence area evolution law analysis is not comprehensive, and the observation data could not be taken in a timely manner. Thus, satisfying the requirements of land reclamation and ecological restoration is difficult.In this study, a dynamic prediction model of the water accumulation area in the high-water-level coal mining subsidence area was established. The equivalent relationship between the volume of the subsidence basin and the water accumulation in the subsidence basin was established using the probability integration method based on the Knothe time function and the principle of water balance in the coal mining subsidence area, combined with the main influencing factors of the evolution of the water accumulation area to determine the mining capacity of the working face. The dynamic evolution of the surface water range was predicted, and information, such as the range of water accumulation, the range of subsidence basins, the level of water accumulation, and the amount of water accumulation during the evolution of water accumulation is obtained. The general law of the evolution of the water accumulation area in the coal mining subsidence area was studied. The time of water accumulation in subsidence basins generally lags behind the mining of the working face. The evolution of surface water accumulation can be divided into four stages: unformed period, synchronous growth period, residual growth period, and relatively stable period. The accuracy of the model was analyzed.Therefore, taking the working surface of Wugou Mine 1031 as the research area, the three-dimensional (3D) spatial information of the coal mining subsidence water accumulation area of 1031 was dynamically predicted, and the dynamic evolution law of the water accumulation range was analyzed. Initially, the study area overview, as shown in section2, was introduced. Then, the overall idea of the dynamic prediction method was proposed, as shown in section 3. The complete dynamic prediction method of the 3D spatial information of high diving mining subsidence is introduced in section 4. Finally, in section 5, model evaluation and validation are described by the engineering example of 1031. This study provides a theoretical basis for land use planning, land reclamation, and the establishment of a land and water complex ecosystem in the subsidence area of high diving sites.


## 2. Overview of the Study Area

Wugou Coal Mine is the 10th pair of mines developed and constructed by North Anhui Coal and Power Group. It is located in the Wugou Town, Suixi County, which is 50 km away from the Huaibei urban area, as shown in Figure 2. The underground diving level in the mining area is 1.5 m, which is a typical high diving level mining area. The coal-containing formation in this mine is carbonite-Permian system, mainly comprising four layers of coal seam (72, 81, 82, and 10 coal seams, with the average thicknesses of 1.80, 2.50, 2.53, and 3.86 m, respectively). The average recoverable total thickness is 10.69 m because the coal seam is seriously damaged by fault cutting. A considerable part of the shallow area is in the wind oxidation zone or waterproof coal column, resulting in the recoverable area of the coal seam. The total resource reserves of the entire well field are 137,304,600 tonnes, 125,910 million tonnes, and 40,339,500 tonnes, of which 10 coal seams’ recoverable reserves are 25,956,200 tonnes, accounting for 64.3%.Figure 2
Location of the study area and current status of the water system.The regional climate is mild. It is a monsoon warm temperate semihumid climate, mild spring and autumn rain, hot and rainy in summer, and cold and windy in winter. The annual average temperature is 14.1°C, the annual average sunshine is 2345.3 h, the average annual rainfall is 834 mm, the groundwater is buried at a depth range of 1.5–4 m, and the average annual buried depth of the water level is 2.48 m [25]. The rainfall is mostly concentrated in July and August, the annual evaporation is 1400 mm, the annual frost-free period is 208–220 days, and the freezing period is generally from early December to mid-February of the following year.Figure2 shows the location of the working surface and surface remote sensing images with a total length of 1136 and width of 180 m. The coal seam inclinations are from 3° to 16° and 8°. The coal seam average thickness is 3.8 m, and the average buried depth is 364 m. Figure 3 shows the loose layer at approximately 94 m, as shown in Figure 3. The working surface is flat with a surface elevation of 26.50 m to 27.56 m, with no large flow river on the surface. The water on the surface forms after the working surface mining because of the buried depth of the groundwater in the area, high rainfall, and the large subsidence coefficient of the surface.Figure 3
Generalized columns of working faces for 1031.


## 3. Overall Idea of the Prediction Approach

### 3.1. Obtaining the Geological Mining Parameters

According to the working surface specifications of the mining area, the following parameters were established, namely mining promotion distancev, the tangent of major influence angle tan β, average mining thickness of the working surface m, coal bed pitch α, subsidence factor q, and major influence radius r. The 3D coordinate system of the mining subsidence space was established. The inflection point of the subsidence curve of the main section on the surface includes the coordinate origin o. The x-axis points along the surface toward the drain area are parallel to the coal seam direction. The y-axis over the coordinate origin leads straight down from the vertical z-axis of the x-axis. The Knothe time function was used to estimate the mining subsidence model.


### 3.2. Calculating the Subsidence Basin Volume

The dynamic subsidence was interpolated by the Kriging interpolation method, the contour of the subsidence basin of the coal mining subsidence area was generated successively, and the subsidence basin volume corresponding to different contour sections was calculated.


### 3.3. Calculating the Water Accumulation Volume in the Subsidence Basin

The hydrological and water resource data of the coal mining subsidence area were collected, and the water balance iterative equation was established according to the water balance principle to obtain the water volume at different times in the subsidence basin of the coal mining subsidence area.


### 3.4. Calculating the 3D Spatial Information

The water volume in different sinking basins was calculated according to the iterative equation of water balance. In addition, combined with the sinking basin volume corresponding to different contour sections at this time, the relationship between the subsidence basin volume and the volume of water in the basin was established.This model was used to predict the 3D dynamic spatial information at different moments with the advancement of the working surface, including water depth, water area, water range, water volume, and maximum storage capacity.The technical flow chart is shown in Figure4.Figure 4
Technical flow chart.


## 3.1. Obtaining the Geological Mining Parameters

According to the working surface specifications of the mining area, the following parameters were established, namely mining promotion distancev, the tangent of major influence angle tan β, average mining thickness of the working surface m, coal bed pitch α, subsidence factor q, and major influence radius r. The 3D coordinate system of the mining subsidence space was established. The inflection point of the subsidence curve of the main section on the surface includes the coordinate origin o. The x-axis points along the surface toward the drain area are parallel to the coal seam direction. The y-axis over the coordinate origin leads straight down from the vertical z-axis of the x-axis. The Knothe time function was used to estimate the mining subsidence model.


## 3.2. Calculating the Subsidence Basin Volume

The dynamic subsidence was interpolated by the Kriging interpolation method, the contour of the subsidence basin of the coal mining subsidence area was generated successively, and the subsidence basin volume corresponding to different contour sections was calculated.


## 3.3. Calculating the Water Accumulation Volume in the Subsidence Basin

The hydrological and water resource data of the coal mining subsidence area were collected, and the water balance iterative equation was established according to the water balance principle to obtain the water volume at different times in the subsidence basin of the coal mining subsidence area.


## 3.4. Calculating the 3D Spatial Information

The water volume in different sinking basins was calculated according to the iterative equation of water balance. In addition, combined with the sinking basin volume corresponding to different contour sections at this time, the relationship between the subsidence basin volume and the volume of water in the basin was established.This model was used to predict the 3D dynamic spatial information at different moments with the advancement of the working surface, including water depth, water area, water range, water volume, and maximum storage capacity.The technical flow chart is shown in Figure4.Figure 4
Technical flow chart.


## 4. Dynamic Prediction Method of 3D Spatial Information

### 4.1. Dynamic Sinking Is Expected

#### 4.1.1. Random Medium Theory

The theoretical basis of the probability integration method is the stochastic medium theory. Thus, this method is also known as the stochastic medium method. Littvinishen, the founder of random media theory, is a doctor in mathematics, who is familiar with probability theory. It is impossible to identify whether the rock mass medium is elastic, plastic, continuous, loose, or whole. The mine rock mass medium is generally called random medium. The basic assumptions and corollary of random medium theory indicate that the rock medium comprises many sufficiently small rock block particles. These particles are completely disconnected and can move relatively. Their movement is a stochastic process, as shown in Figure5. The laws of rock strata and surface movement caused by mining are macroscopically similar to those described in the granular medium model acting as a random medium. The mutual verification of theory and practice was carried out, and the effect function based on the random medium theory is very similar to the one-Gaussian function of the Knott–Bradlake theory based on the measured data.Figure 5
Random medium model. (a) Theoretical model. (b) Stochastic model of granular medium movement.
(a)(b)


#### 4.1.2. Probabilistic Integral Method

The randomastic medium theory was initially introduced into the rock layer movement study by the Polish scholar Li Tvinishen in the 1950s and was then developed into the probability integration method by Chinese scholars Liu Baochen and Liao Guohua. After more than 20 years of research by mining and subsidence workers in China, the probability integration method has become one of the more mature and most widely used prediction methods in China. Probability integral is named after the probability integral (or its derivatives) in the movement and deformation prediction formula.This method is based on the movement law of the surface and is generally similar to the random medium model. It is a comprehensive response to geographical and mining conditions to mining settlements, such as mining depth, mining thickness, rock mass properties, and coal seam inclination, regardless of the effect of special structures, such as faults, folds, and collapse columns.The probabilistic integration model is one of the most widely used models in the field of mining subsidence currently as it can predict the degree of surface movement deformation through the working surface geological mining condition parameters and the expected parameters. The probability integral method model is represented by equation (1).(1)Wx,y=WmaxCxCy,Wmax=mqcosα,Cx=122π∫0πx⋅tanβ/H0e−λ2dλ−2π∫0πx−l⋅tanβ/H0e−λ2dλ,Cy=122π∫0πy⋅tanβ/H1e−λ2dλ−2π∫0πy−L⋅tanβ/H0e−λ2dλ,l=L0−S1−S2,L=L1−S3−S4×cosα+sinα⋅cotθ0,where x and y are the coordinates of any point in the working surface coordinate system within the influence range of underground mining, and Wx,y is the final sinking amount of the point (x,y) in the working surface coordinate system. Wmax is the surface maximum subsidence, when the goaf is critical or supercritical in size [26]. m is the mining thickness. q is the sinking coefficient. H1 and H2 are the mining depths of downhill and uphill direction, respectively. H0 is the average mining depth. θ0 indicates the propagation angle influence.tanβ is the positive tangent value of the main influence angle. L0 and L1 are the working surface direction and tendency length, respectively. S1, S2, S3, and S4 are the inflection point deviation distances corresponding to both ends of the uphill and downhill directions. α is the seam inclination, and the practical significance of the parameters can be used. The parameters m, H0, H1, H2, L0, L1, and α are called geological mining condition parameters and are generally provided by the production unit, and they are the known quantities: q, θ, tanβ, S1, S2, S3, and S4 are the expected parameters of the probability integration method, and they can be obtained from the measured data. With the known geological mining condition parameters of the working surface, the sinking amount of any point in the working surface can be expected by equation (1) as long as the rate integration prediction parameter is accounted into the model.


#### 4.1.3. Knothe Time Function

To more accurately reverse and predict the dynamic process of surface subsidence, mine staff and scholars attempt to simulate the process of the dynamic subsidence of surface points using mathematical models, amongst which the most widely used model is the Knothe function of time. According to the working surface information of the mining area, the dynamic subsidence in the basin was obtained by the probability integration method based on the Knothe time function.In 1952, the Polish scholar Knothe established the Knothe time function by analyzing the relationship between the sinking speed at a certain surface point and the sinking amount at the final subsidence at that time. This function can partly reflect the relationship between the surface subsidence and time [24]. The Knothe model represents the relationship between the surface subsidence velocity and subsidence volume, as expressed by equation (2).(2)dWtdt=vt=CW0−Wt,where t is the surface subsidence moment, W0 is the final surface subsidence amount, W(t) is the subsidence amount of the surface at time t and V(t) is the subsidence speed of the surface point at time t, C is the time factor influence coefficient or lithology factor and is generally considered to be related to the mechanical nature of the overlying rock layer [27].The surface sink is 0 att = 0, namely, W(0) = 0. Thus, equation (3) was integrated with respect to t as follows:(3)Wt=W01−e−Ct,where t is the independent variable, representing the time calculated from the moment the surface begins to sink, W0 is the final amount of surface subsidence and can be obtained or estimated by the probabilistic integration model, C value is the time factor influence coefficient and is related to the geological mining conditions of the working surface, and equation (3) is the Knothe time function. Substituting the final surface subsidence W0 and the time factor influence coefficient C into equation (3) can lead to the surface subsidence at any time. In addition, the horizontal movement, slope, and parameters like curvature can be calculated.Whent tends to zero, W(t) = 0; when t⟶∞, W(t) = W0. Thus, W0 controls the upper and lower boundary limits of the Knothe time function and is obtained by changing the time factor. C is the Knothe time function impact factor.


#### 4.1.4. Dynamic Sinking Is Expected

The probability integral method is based on the Knothe time function, and geological mining condition parameters and equation (3) are used to calculate the dynamic subsidence prediction value of the subsidence basin, which is represented as follows:(4)wxi,yi,T=mqcosαr2⋅1−e−cT⋅∬De−πxi−ς2+yi−τ2/r2dςdτ,where wxi,yi,T is the expected value of the dynamic subsidence at any point T on the surface of the subsidence basin, m is the average thickness of the working surface per m, α is the inclined angle of the coal seam, q is the sinking coefficient, r is the main influence radius and is measured per meter, c is the time coefficient, e is the natural constant and is approximately 2.718, dξ and dξ are the binary integral variables, ξ and τ are the integral variables in the s and y directions, respectively, and D is the range of the sinking basin at time T.


### 4.2. Volume Calculation of the Subsidence Basin

The Kriging interpolation method is an unbiased optimal estimation of regionalization variables in finite regions based on the variant function theory and structural analysis [27] and was established back in the 1950s with the prototype algorithm and is expressed as follows:(5)z0∧=∑i=0nλizi,where z0 is the value of the estimated point and zi is the actual measurement of the interpolation element at point i. λi= the unknown weight of the actual measurement at the position. The weights are not only related to the distance between the measurement points and the predicted location but also based on the overall spatial arrangement of the measurement points. In the Kriging method, weight λi is a fit model of the optimal set of coefficients that can satisfy the smallest difference between the estimates and the actual values at the unknown point, satisfying the unbiased estimation conditions, depending on the spatial relationship amongst the measured points, the distance of the predicted position, and the measurements around the predicted location [27]. Therefore, the Kriging interpolation method considers not only the known relationships of the points to be estimated but also the spatial correlations of the variables. At present, this interpolation method can theoretically estimate the model error point-by-point without using the error test model, and it is not theoretically possible by other interpolation methods [28].According to the expected value of the subsidence basin, the interpolation method initially uses formula (6) to calculate the difference between the bottom of the sinking basin and the lowest horizontal section in the basin d′, as follows:(6)d′=Wm−10+nd,where n is the number of horizontal sections, except the edge of the sinking basin, Wm is the maximum subsidence value of the surface at any time, and the contour value of 10, 10 + d, 10 + 2d, …, 10 + nd in mm corresponds to the equal section area of different sinking basins S, S1, S2, S3, …, Sn. The volume of the sinking basin corresponding to the different contour sections was calculated using the area of the isoline section of each sinking basin.Secondly, as shown in Figure6(a), V is the volume of the entire sinking basin, V1 is the volume between S and S1, and V2 is the volume between S1 and S2. Vn is the sinking volume between Sn−1 and Sn, and V′ is the sinking volume between Sn and Wm. Thus, the volume of the sinking basin with different contour cross-sections can be obtained, VmaxT is the maximum volume of the basin that can be obtained by formula (7).(7)VmaxT=V1+V2+V3+⋯+Vn+V'=d2S+dS1+dS2+⋯+dSn+d2+d′3Sn,Vbk=V′+Vn+Vn−1+Vn−2+⋯+Vk=d2+d′3Sn+dSn−1+⋯+d2Sk,1≤k≤n,where VmaxT is the maximum volume of the basin, Vbk is the sinking basin volume corresponding to the cross-section of the contour 10 + (k − 1)d, V1 is the volume between the contour 10 and the contour 10 + d is the cross-section, V2 is the volume between the contour 10 + d and the contour 10 + 2d cross-section, Vn is the volume between the contour 10 + (n − 1)d and the contour 10 + nd cross-section, V′ is the volume between the contour of 10 + nd and the maximum sinking point Wm.Figure 6
(a) Volume calculation of subsidence basin. (b) Schematic diagram of the working surface.
(a)(b)As the face advances, the newly mined rock begins to move. When the scope of the goaf is large enough after the impact of mining reaches the surface, the ground affected by the mining will sink down from the original elevation, thereby the surface above the goaf will form a much larger depression than the area of the goaf. Such depressions are often called surface mobile basins or ground subsidences, as shown in Figure6(b).Finally, the horizontal cross-sectional area of each contour of the sinking basin is calculated using formula (8), as follows:(8)Sn=12∑i=1Nxi10+nd⋅yi+110+nd−xi+110+nd⋅yi10+nd,where xi10+nd and yi10+nd are the contour point coordinates of the sink basin boundary of the sink value 10 + nd (n is a positive integer in unit mm).


### 4.3. Calculation of Water Accumulation in the Subsidence Basin

The supplementary and the excretion items are mainly considered when calculating the water volume changes in the basin, and the factors with less influence on the evolution of the water range are not considered. Supplementary items include water surface precipitation, surface runoff, and groundwater recharge. Excretion items include water surface evaporation and permeable water volume.(1)
Water face precipitation refers to the amount of water that directly falls in the water surface area during the atmospheric precipitation process, and this water is directly supplied to the atmospheric precipitation area, and the water surface precipitation can be determined according to the product of the water area and precipitation, which is as follows:(9)P=Swater⋅Pday⋅t1000,whereP is the water surface precipitation/m3, Swater is the water surface area/m2, Pday is the daily precipitation/mm, and t is the duration/d.(2)
Surface runoff refers to the water flowing into the basin during the atmospheric precipitation. The flow area of the subsidence basin and the flow coefficient of the flow area of the basin and the study area were obtained as follows:(10)R=Sbasin−Swater⋅Pday⋅α⋅t1000,whereR is the surface production flow/m3 in the unflooded area, Sbasin is the sinking basin area/m2, Swater is the water surface area/m2, Pday is the daily precipitation/mm, and α is the surface runoff coefficient. According to the relevant hydrological data in the literature [29], the surface runoff coefficient in Huaibei area is about 0.15–0.25. According to the actual situation of the study area, the underlying surface of the mining area is mainly farmland, crops, and vegetation, and there are small areas of the hardened surface. The surface runoff coefficient is taken as 0.2 based on the actual situation, and t is the number of days/d.(3)
Groundwater replenishment refers to the process of an aquifer or an aquifer system acquiring water from the outside world, which is represented as follows:(11)Wunderground=0.042⋅P+R,whereP is the water precipitation/m3 and R is the surface production flow/m3 in the unflooded area.(4)
Water evaporation refers to the evaporation of the water area in the subsidence basin. The product of the water area and evaporation is obtained as follows:(12)E=Swater⋅Eday⋅t1000,whereE is the water evaporation/m3, Swater is the water surface area/m2, Eday is the daily evaporation/mm, and t is the number of days/d.(5)
Water leakage from the subsidence basin into the groundwater is computed as follows:(13)Qseepage=0.09E=Swater⋅Eday⋅t1000⋅0.09,whereQseepage is the water seepage (m3), E is the water evaporation (m3), Eday is the daily evaporation (mm), Swater is the water surface area (m2), and t is the days (d).Supplementary items and excretion items (also known as expenses and receipts) are the principles of water balance. Water balance refers to the difference between the amount of water in revenue and the amount of expenditure in any area (or water body) during any period of time, and it is equal to the amount of change in water storage, as shown in formula (14).(14)Vponding=V+P+R+Wunderground−E,where V is the original water amount, P is the water surface precipitation, R is the surface runoff, Wunderground is the groundwater recharge, and E is the water surface evaporation amount.The iterative model of water balance was established using the principle of water balance, as shown in formula (15), to calculate the volume of water accumulation at different times T in the coal mining subsidence area, which is as follows:(15)VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt,where VwT+Δt is the expected water volume in the subsidence basin in cubic meter, VwT is the original water volume in the subsidence basin, Δt is the time interval, ΔP is the daily precipitation per meter, ΔE is the daily evaporation quantity in meter, α is the surface runoff coefficient, SbT is the subsidence basin boundary contour section area of 10 mm, its unit is square, SwT is the water area in the subsidence basin, its unit is square, QgwΔt is the groundwater seepage, its unit is cubic meter, RΔt is the surface river and subsidence area water exchange, its unit is cubic meter, and MΔt is the manually extracted volume in cubic meter.


### 4.4. 3D Dynamic Spatial Information Calculation

According to the water volume in the sinking basin at different times corresponding to different contour line sections, the relationship between the subsidence basin volume and water volume was established. The contours of the sinking basin are equal to the volume of the accumulated water in the sinking basin. Three-dimensional dynamic spatial information, such as water depth, water area, water range, water volume, and maximum storage capacity of the subsidence basin at different timesT as the working surface is predicted, as shown in Figure 7.Figure 7
Schematic diagram of spatial three-dimensional coordinate system.According to the water volume in the sinking basin at different times and the sinking basin volume corresponding to the different contour sections, a mathematical relationship between the volume of the subsidence basin and the water volume in the basin was established as follows:(16)VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt=d2+d′3⋅Sn+d⋅Sn−1+…+d2⋅Sk=Vbk,1≤k≤n,where VwT+Δt is the expected water volume in the sink basin, (m3) and Vbk is the sink basin volume corresponding to the contour 10 + (k − 1)d section, (m3). From VWT+Δt=Vbk, the contour of the sinking basin 10 + (k − 1)d is equal to the volume of water in the sinking basin.Therefore, the 3D spatial information of the surface water area at different timesT as the working surface advances was obtained by formula (17).(17)HWT=Wm−10−k−1d,SwT=12∑i=1Nxi10+k−1d⋅yi+110+k−1d−xi+110+k−1d⋅yi10+k−1d,VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt,VmaxT=V1+V2+V3+...+Vn+V′=d2S+dS1+dS2...+dSn+d2+d′3Sn,where HwT is the water depth of the subsidence basin at dynamically predicted time T, (m), Wm is the maximum sinking value of the sinking basin at time T, (m), SwT is the water area and unit square of the subsidence basin at dynamically predicted time T, xi10+k−1d and yi10+k−1d are the contour point coordinates with a sinking value of 10 + (k − 1)d, VwT+Δt is the expected water volume in the subsidence basin, (m3), VwT is the original water volume of the subsidence basin, (m3), Δt is the time interval, ΔP is the daily average precipitation, (mm), ΔE is the daily average evaporation, (m), α is the surface runoff coefficient, SbT is the cross-sectional area of the boundary contour of the subsidence basin with a subsidence of 10 mm, (m2), SwT is the area of water in the subsidence basin, (m2), QgwΔt is the groundwater seepage flow, (m3), RΔt is the exchange volume of the surface river with the accumulated water in the subsidence area, (m3), MΔt is the manual withdrawal volume, (m3), VmaxT is the maximum reservoir capacity of the subsidence basin at the dynamic predicted time T, (m3), V1 is the volume between S and S1, V2 is the volume between S1 and S2, Vn is the volume between Sn−1 and Sn, and V′ is Sn the volume between the maximum sinking points Wm. The boundary of the stagnant water range is the curve enclosed by the coordinate points xi10+k−1d and yi10+k−1d on the isoline with a sinking value of 10 + (k − 1)d.


## 4.1. Dynamic Sinking Is Expected

### 4.1.1. Random Medium Theory

The theoretical basis of the probability integration method is the stochastic medium theory. Thus, this method is also known as the stochastic medium method. Littvinishen, the founder of random media theory, is a doctor in mathematics, who is familiar with probability theory. It is impossible to identify whether the rock mass medium is elastic, plastic, continuous, loose, or whole. The mine rock mass medium is generally called random medium. The basic assumptions and corollary of random medium theory indicate that the rock medium comprises many sufficiently small rock block particles. These particles are completely disconnected and can move relatively. Their movement is a stochastic process, as shown in Figure5. The laws of rock strata and surface movement caused by mining are macroscopically similar to those described in the granular medium model acting as a random medium. The mutual verification of theory and practice was carried out, and the effect function based on the random medium theory is very similar to the one-Gaussian function of the Knott–Bradlake theory based on the measured data.Figure 5
Random medium model. (a) Theoretical model. (b) Stochastic model of granular medium movement.
(a)(b)


### 4.1.2. Probabilistic Integral Method

The randomastic medium theory was initially introduced into the rock layer movement study by the Polish scholar Li Tvinishen in the 1950s and was then developed into the probability integration method by Chinese scholars Liu Baochen and Liao Guohua. After more than 20 years of research by mining and subsidence workers in China, the probability integration method has become one of the more mature and most widely used prediction methods in China. Probability integral is named after the probability integral (or its derivatives) in the movement and deformation prediction formula.This method is based on the movement law of the surface and is generally similar to the random medium model. It is a comprehensive response to geographical and mining conditions to mining settlements, such as mining depth, mining thickness, rock mass properties, and coal seam inclination, regardless of the effect of special structures, such as faults, folds, and collapse columns.The probabilistic integration model is one of the most widely used models in the field of mining subsidence currently as it can predict the degree of surface movement deformation through the working surface geological mining condition parameters and the expected parameters. The probability integral method model is represented by equation (1).(1)Wx,y=WmaxCxCy,Wmax=mqcosα,Cx=122π∫0πx⋅tanβ/H0e−λ2dλ−2π∫0πx−l⋅tanβ/H0e−λ2dλ,Cy=122π∫0πy⋅tanβ/H1e−λ2dλ−2π∫0πy−L⋅tanβ/H0e−λ2dλ,l=L0−S1−S2,L=L1−S3−S4×cosα+sinα⋅cotθ0,where x and y are the coordinates of any point in the working surface coordinate system within the influence range of underground mining, and Wx,y is the final sinking amount of the point (x,y) in the working surface coordinate system. Wmax is the surface maximum subsidence, when the goaf is critical or supercritical in size [26]. m is the mining thickness. q is the sinking coefficient. H1 and H2 are the mining depths of downhill and uphill direction, respectively. H0 is the average mining depth. θ0 indicates the propagation angle influence.tanβ is the positive tangent value of the main influence angle. L0 and L1 are the working surface direction and tendency length, respectively. S1, S2, S3, and S4 are the inflection point deviation distances corresponding to both ends of the uphill and downhill directions. α is the seam inclination, and the practical significance of the parameters can be used. The parameters m, H0, H1, H2, L0, L1, and α are called geological mining condition parameters and are generally provided by the production unit, and they are the known quantities: q, θ, tanβ, S1, S2, S3, and S4 are the expected parameters of the probability integration method, and they can be obtained from the measured data. With the known geological mining condition parameters of the working surface, the sinking amount of any point in the working surface can be expected by equation (1) as long as the rate integration prediction parameter is accounted into the model.


### 4.1.3. Knothe Time Function

To more accurately reverse and predict the dynamic process of surface subsidence, mine staff and scholars attempt to simulate the process of the dynamic subsidence of surface points using mathematical models, amongst which the most widely used model is the Knothe function of time. According to the working surface information of the mining area, the dynamic subsidence in the basin was obtained by the probability integration method based on the Knothe time function.In 1952, the Polish scholar Knothe established the Knothe time function by analyzing the relationship between the sinking speed at a certain surface point and the sinking amount at the final subsidence at that time. This function can partly reflect the relationship between the surface subsidence and time [24]. The Knothe model represents the relationship between the surface subsidence velocity and subsidence volume, as expressed by equation (2).(2)dWtdt=vt=CW0−Wt,where t is the surface subsidence moment, W0 is the final surface subsidence amount, W(t) is the subsidence amount of the surface at time t and V(t) is the subsidence speed of the surface point at time t, C is the time factor influence coefficient or lithology factor and is generally considered to be related to the mechanical nature of the overlying rock layer [27].The surface sink is 0 att = 0, namely, W(0) = 0. Thus, equation (3) was integrated with respect to t as follows:(3)Wt=W01−e−Ct,where t is the independent variable, representing the time calculated from the moment the surface begins to sink, W0 is the final amount of surface subsidence and can be obtained or estimated by the probabilistic integration model, C value is the time factor influence coefficient and is related to the geological mining conditions of the working surface, and equation (3) is the Knothe time function. Substituting the final surface subsidence W0 and the time factor influence coefficient C into equation (3) can lead to the surface subsidence at any time. In addition, the horizontal movement, slope, and parameters like curvature can be calculated.Whent tends to zero, W(t) = 0; when t⟶∞, W(t) = W0. Thus, W0 controls the upper and lower boundary limits of the Knothe time function and is obtained by changing the time factor. C is the Knothe time function impact factor.


### 4.1.4. Dynamic Sinking Is Expected

The probability integral method is based on the Knothe time function, and geological mining condition parameters and equation (3) are used to calculate the dynamic subsidence prediction value of the subsidence basin, which is represented as follows:(4)wxi,yi,T=mqcosαr2⋅1−e−cT⋅∬De−πxi−ς2+yi−τ2/r2dςdτ,where wxi,yi,T is the expected value of the dynamic subsidence at any point T on the surface of the subsidence basin, m is the average thickness of the working surface per m, α is the inclined angle of the coal seam, q is the sinking coefficient, r is the main influence radius and is measured per meter, c is the time coefficient, e is the natural constant and is approximately 2.718, dξ and dξ are the binary integral variables, ξ and τ are the integral variables in the s and y directions, respectively, and D is the range of the sinking basin at time T.


## 4.1.1. Random Medium Theory

The theoretical basis of the probability integration method is the stochastic medium theory. Thus, this method is also known as the stochastic medium method. Littvinishen, the founder of random media theory, is a doctor in mathematics, who is familiar with probability theory. It is impossible to identify whether the rock mass medium is elastic, plastic, continuous, loose, or whole. The mine rock mass medium is generally called random medium. The basic assumptions and corollary of random medium theory indicate that the rock medium comprises many sufficiently small rock block particles. These particles are completely disconnected and can move relatively. Their movement is a stochastic process, as shown in Figure5. The laws of rock strata and surface movement caused by mining are macroscopically similar to those described in the granular medium model acting as a random medium. The mutual verification of theory and practice was carried out, and the effect function based on the random medium theory is very similar to the one-Gaussian function of the Knott–Bradlake theory based on the measured data.Figure 5
Random medium model. (a) Theoretical model. (b) Stochastic model of granular medium movement.
(a)(b)


## 4.1.2. Probabilistic Integral Method

The randomastic medium theory was initially introduced into the rock layer movement study by the Polish scholar Li Tvinishen in the 1950s and was then developed into the probability integration method by Chinese scholars Liu Baochen and Liao Guohua. After more than 20 years of research by mining and subsidence workers in China, the probability integration method has become one of the more mature and most widely used prediction methods in China. Probability integral is named after the probability integral (or its derivatives) in the movement and deformation prediction formula.This method is based on the movement law of the surface and is generally similar to the random medium model. It is a comprehensive response to geographical and mining conditions to mining settlements, such as mining depth, mining thickness, rock mass properties, and coal seam inclination, regardless of the effect of special structures, such as faults, folds, and collapse columns.The probabilistic integration model is one of the most widely used models in the field of mining subsidence currently as it can predict the degree of surface movement deformation through the working surface geological mining condition parameters and the expected parameters. The probability integral method model is represented by equation (1).(1)Wx,y=WmaxCxCy,Wmax=mqcosα,Cx=122π∫0πx⋅tanβ/H0e−λ2dλ−2π∫0πx−l⋅tanβ/H0e−λ2dλ,Cy=122π∫0πy⋅tanβ/H1e−λ2dλ−2π∫0πy−L⋅tanβ/H0e−λ2dλ,l=L0−S1−S2,L=L1−S3−S4×cosα+sinα⋅cotθ0,where x and y are the coordinates of any point in the working surface coordinate system within the influence range of underground mining, and Wx,y is the final sinking amount of the point (x,y) in the working surface coordinate system. Wmax is the surface maximum subsidence, when the goaf is critical or supercritical in size [26]. m is the mining thickness. q is the sinking coefficient. H1 and H2 are the mining depths of downhill and uphill direction, respectively. H0 is the average mining depth. θ0 indicates the propagation angle influence.tanβ is the positive tangent value of the main influence angle. L0 and L1 are the working surface direction and tendency length, respectively. S1, S2, S3, and S4 are the inflection point deviation distances corresponding to both ends of the uphill and downhill directions. α is the seam inclination, and the practical significance of the parameters can be used. The parameters m, H0, H1, H2, L0, L1, and α are called geological mining condition parameters and are generally provided by the production unit, and they are the known quantities: q, θ, tanβ, S1, S2, S3, and S4 are the expected parameters of the probability integration method, and they can be obtained from the measured data. With the known geological mining condition parameters of the working surface, the sinking amount of any point in the working surface can be expected by equation (1) as long as the rate integration prediction parameter is accounted into the model.


## 4.1.3. Knothe Time Function

To more accurately reverse and predict the dynamic process of surface subsidence, mine staff and scholars attempt to simulate the process of the dynamic subsidence of surface points using mathematical models, amongst which the most widely used model is the Knothe function of time. According to the working surface information of the mining area, the dynamic subsidence in the basin was obtained by the probability integration method based on the Knothe time function.In 1952, the Polish scholar Knothe established the Knothe time function by analyzing the relationship between the sinking speed at a certain surface point and the sinking amount at the final subsidence at that time. This function can partly reflect the relationship between the surface subsidence and time [24]. The Knothe model represents the relationship between the surface subsidence velocity and subsidence volume, as expressed by equation (2).(2)dWtdt=vt=CW0−Wt,where t is the surface subsidence moment, W0 is the final surface subsidence amount, W(t) is the subsidence amount of the surface at time t and V(t) is the subsidence speed of the surface point at time t, C is the time factor influence coefficient or lithology factor and is generally considered to be related to the mechanical nature of the overlying rock layer [27].The surface sink is 0 att = 0, namely, W(0) = 0. Thus, equation (3) was integrated with respect to t as follows:(3)Wt=W01−e−Ct,where t is the independent variable, representing the time calculated from the moment the surface begins to sink, W0 is the final amount of surface subsidence and can be obtained or estimated by the probabilistic integration model, C value is the time factor influence coefficient and is related to the geological mining conditions of the working surface, and equation (3) is the Knothe time function. Substituting the final surface subsidence W0 and the time factor influence coefficient C into equation (3) can lead to the surface subsidence at any time. In addition, the horizontal movement, slope, and parameters like curvature can be calculated.Whent tends to zero, W(t) = 0; when t⟶∞, W(t) = W0. Thus, W0 controls the upper and lower boundary limits of the Knothe time function and is obtained by changing the time factor. C is the Knothe time function impact factor.


## 4.1.4. Dynamic Sinking Is Expected

The probability integral method is based on the Knothe time function, and geological mining condition parameters and equation (3) are used to calculate the dynamic subsidence prediction value of the subsidence basin, which is represented as follows:(4)wxi,yi,T=mqcosαr2⋅1−e−cT⋅∬De−πxi−ς2+yi−τ2/r2dςdτ,where wxi,yi,T is the expected value of the dynamic subsidence at any point T on the surface of the subsidence basin, m is the average thickness of the working surface per m, α is the inclined angle of the coal seam, q is the sinking coefficient, r is the main influence radius and is measured per meter, c is the time coefficient, e is the natural constant and is approximately 2.718, dξ and dξ are the binary integral variables, ξ and τ are the integral variables in the s and y directions, respectively, and D is the range of the sinking basin at time T.


## 4.2. Volume Calculation of the Subsidence Basin

The Kriging interpolation method is an unbiased optimal estimation of regionalization variables in finite regions based on the variant function theory and structural analysis [27] and was established back in the 1950s with the prototype algorithm and is expressed as follows:(5)z0∧=∑i=0nλizi,where z0 is the value of the estimated point and zi is the actual measurement of the interpolation element at point i. λi= the unknown weight of the actual measurement at the position. The weights are not only related to the distance between the measurement points and the predicted location but also based on the overall spatial arrangement of the measurement points. In the Kriging method, weight λi is a fit model of the optimal set of coefficients that can satisfy the smallest difference between the estimates and the actual values at the unknown point, satisfying the unbiased estimation conditions, depending on the spatial relationship amongst the measured points, the distance of the predicted position, and the measurements around the predicted location [27]. Therefore, the Kriging interpolation method considers not only the known relationships of the points to be estimated but also the spatial correlations of the variables. At present, this interpolation method can theoretically estimate the model error point-by-point without using the error test model, and it is not theoretically possible by other interpolation methods [28].According to the expected value of the subsidence basin, the interpolation method initially uses formula (6) to calculate the difference between the bottom of the sinking basin and the lowest horizontal section in the basin d′, as follows:(6)d′=Wm−10+nd,where n is the number of horizontal sections, except the edge of the sinking basin, Wm is the maximum subsidence value of the surface at any time, and the contour value of 10, 10 + d, 10 + 2d, …, 10 + nd in mm corresponds to the equal section area of different sinking basins S, S1, S2, S3, …, Sn. The volume of the sinking basin corresponding to the different contour sections was calculated using the area of the isoline section of each sinking basin.Secondly, as shown in Figure6(a), V is the volume of the entire sinking basin, V1 is the volume between S and S1, and V2 is the volume between S1 and S2. Vn is the sinking volume between Sn−1 and Sn, and V′ is the sinking volume between Sn and Wm. Thus, the volume of the sinking basin with different contour cross-sections can be obtained, VmaxT is the maximum volume of the basin that can be obtained by formula (7).(7)VmaxT=V1+V2+V3+⋯+Vn+V'=d2S+dS1+dS2+⋯+dSn+d2+d′3Sn,Vbk=V′+Vn+Vn−1+Vn−2+⋯+Vk=d2+d′3Sn+dSn−1+⋯+d2Sk,1≤k≤n,where VmaxT is the maximum volume of the basin, Vbk is the sinking basin volume corresponding to the cross-section of the contour 10 + (k − 1)d, V1 is the volume between the contour 10 and the contour 10 + d is the cross-section, V2 is the volume between the contour 10 + d and the contour 10 + 2d cross-section, Vn is the volume between the contour 10 + (n − 1)d and the contour 10 + nd cross-section, V′ is the volume between the contour of 10 + nd and the maximum sinking point Wm.Figure 6
(a) Volume calculation of subsidence basin. (b) Schematic diagram of the working surface.
(a)(b)As the face advances, the newly mined rock begins to move. When the scope of the goaf is large enough after the impact of mining reaches the surface, the ground affected by the mining will sink down from the original elevation, thereby the surface above the goaf will form a much larger depression than the area of the goaf. Such depressions are often called surface mobile basins or ground subsidences, as shown in Figure6(b).Finally, the horizontal cross-sectional area of each contour of the sinking basin is calculated using formula (8), as follows:(8)Sn=12∑i=1Nxi10+nd⋅yi+110+nd−xi+110+nd⋅yi10+nd,where xi10+nd and yi10+nd are the contour point coordinates of the sink basin boundary of the sink value 10 + nd (n is a positive integer in unit mm).


## 4.3. Calculation of Water Accumulation in the Subsidence Basin

The supplementary and the excretion items are mainly considered when calculating the water volume changes in the basin, and the factors with less influence on the evolution of the water range are not considered. Supplementary items include water surface precipitation, surface runoff, and groundwater recharge. Excretion items include water surface evaporation and permeable water volume.(1)
Water face precipitation refers to the amount of water that directly falls in the water surface area during the atmospheric precipitation process, and this water is directly supplied to the atmospheric precipitation area, and the water surface precipitation can be determined according to the product of the water area and precipitation, which is as follows:(9)P=Swater⋅Pday⋅t1000,whereP is the water surface precipitation/m3, Swater is the water surface area/m2, Pday is the daily precipitation/mm, and t is the duration/d.(2)
Surface runoff refers to the water flowing into the basin during the atmospheric precipitation. The flow area of the subsidence basin and the flow coefficient of the flow area of the basin and the study area were obtained as follows:(10)R=Sbasin−Swater⋅Pday⋅α⋅t1000,whereR is the surface production flow/m3 in the unflooded area, Sbasin is the sinking basin area/m2, Swater is the water surface area/m2, Pday is the daily precipitation/mm, and α is the surface runoff coefficient. According to the relevant hydrological data in the literature [29], the surface runoff coefficient in Huaibei area is about 0.15–0.25. According to the actual situation of the study area, the underlying surface of the mining area is mainly farmland, crops, and vegetation, and there are small areas of the hardened surface. The surface runoff coefficient is taken as 0.2 based on the actual situation, and t is the number of days/d.(3)
Groundwater replenishment refers to the process of an aquifer or an aquifer system acquiring water from the outside world, which is represented as follows:(11)Wunderground=0.042⋅P+R,whereP is the water precipitation/m3 and R is the surface production flow/m3 in the unflooded area.(4)
Water evaporation refers to the evaporation of the water area in the subsidence basin. The product of the water area and evaporation is obtained as follows:(12)E=Swater⋅Eday⋅t1000,whereE is the water evaporation/m3, Swater is the water surface area/m2, Eday is the daily evaporation/mm, and t is the number of days/d.(5)
Water leakage from the subsidence basin into the groundwater is computed as follows:(13)Qseepage=0.09E=Swater⋅Eday⋅t1000⋅0.09,whereQseepage is the water seepage (m3), E is the water evaporation (m3), Eday is the daily evaporation (mm), Swater is the water surface area (m2), and t is the days (d).Supplementary items and excretion items (also known as expenses and receipts) are the principles of water balance. Water balance refers to the difference between the amount of water in revenue and the amount of expenditure in any area (or water body) during any period of time, and it is equal to the amount of change in water storage, as shown in formula (14).(14)Vponding=V+P+R+Wunderground−E,where V is the original water amount, P is the water surface precipitation, R is the surface runoff, Wunderground is the groundwater recharge, and E is the water surface evaporation amount.The iterative model of water balance was established using the principle of water balance, as shown in formula (15), to calculate the volume of water accumulation at different times T in the coal mining subsidence area, which is as follows:(15)VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt,where VwT+Δt is the expected water volume in the subsidence basin in cubic meter, VwT is the original water volume in the subsidence basin, Δt is the time interval, ΔP is the daily precipitation per meter, ΔE is the daily evaporation quantity in meter, α is the surface runoff coefficient, SbT is the subsidence basin boundary contour section area of 10 mm, its unit is square, SwT is the water area in the subsidence basin, its unit is square, QgwΔt is the groundwater seepage, its unit is cubic meter, RΔt is the surface river and subsidence area water exchange, its unit is cubic meter, and MΔt is the manually extracted volume in cubic meter.


## 4.4. 3D Dynamic Spatial Information Calculation

According to the water volume in the sinking basin at different times corresponding to different contour line sections, the relationship between the subsidence basin volume and water volume was established. The contours of the sinking basin are equal to the volume of the accumulated water in the sinking basin. Three-dimensional dynamic spatial information, such as water depth, water area, water range, water volume, and maximum storage capacity of the subsidence basin at different timesT as the working surface is predicted, as shown in Figure 7.Figure 7
Schematic diagram of spatial three-dimensional coordinate system.According to the water volume in the sinking basin at different times and the sinking basin volume corresponding to the different contour sections, a mathematical relationship between the volume of the subsidence basin and the water volume in the basin was established as follows:(16)VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt=d2+d′3⋅Sn+d⋅Sn−1+…+d2⋅Sk=Vbk,1≤k≤n,where VwT+Δt is the expected water volume in the sink basin, (m3) and Vbk is the sink basin volume corresponding to the contour 10 + (k − 1)d section, (m3). From VWT+Δt=Vbk, the contour of the sinking basin 10 + (k − 1)d is equal to the volume of water in the sinking basin.Therefore, the 3D spatial information of the surface water area at different timesT as the working surface advances was obtained by formula (17).(17)HWT=Wm−10−k−1d,SwT=12∑i=1Nxi10+k−1d⋅yi+110+k−1d−xi+110+k−1d⋅yi10+k−1d,VwT+Δt=VwT+SwTΔPΔt+αSbT−SwTΔPΔt+QgwΔt+RΔt−SwTΔEΔt−MΔt,VmaxT=V1+V2+V3+...+Vn+V′=d2S+dS1+dS2...+dSn+d2+d′3Sn,where HwT is the water depth of the subsidence basin at dynamically predicted time T, (m), Wm is the maximum sinking value of the sinking basin at time T, (m), SwT is the water area and unit square of the subsidence basin at dynamically predicted time T, xi10+k−1d and yi10+k−1d are the contour point coordinates with a sinking value of 10 + (k − 1)d, VwT+Δt is the expected water volume in the subsidence basin, (m3), VwT is the original water volume of the subsidence basin, (m3), Δt is the time interval, ΔP is the daily average precipitation, (mm), ΔE is the daily average evaporation, (m), α is the surface runoff coefficient, SbT is the cross-sectional area of the boundary contour of the subsidence basin with a subsidence of 10 mm, (m2), SwT is the area of water in the subsidence basin, (m2), QgwΔt is the groundwater seepage flow, (m3), RΔt is the exchange volume of the surface river with the accumulated water in the subsidence area, (m3), MΔt is the manual withdrawal volume, (m3), VmaxT is the maximum reservoir capacity of the subsidence basin at the dynamic predicted time T, (m3), V1 is the volume between S and S1, V2 is the volume between S1 and S2, Vn is the volume between Sn−1 and Sn, and V′ is Sn the volume between the maximum sinking points Wm. The boundary of the stagnant water range is the curve enclosed by the coordinate points xi10+k−1d and yi10+k−1d on the isoline with a sinking value of 10 + (k − 1)d.


## 5. Engineering Application

Multiple remote sensing monitoring and field measurement results were used to verify and comprehensively analyze the practicability, accuracy, and application scope of the dynamic prediction model, evaluate the estimated accuracy, and provide scientific basis and technical reference for land reclamation and ecological planning in high diving site coal mining subsidence areas.


### 5.1. Remote Sensing Monitoring and Field Measurement

The surface water in the mining area is spread on long range and for a long period of time, and the location is inconvenient for transportation. Conventional observation methods cannot easily obtain large-scale, long-time series, complete, and unified evolution data of the water range, whereas remote sensing observation methods have short observation periods and large amounts of data and high accuracy. Therefore, remote sensing data are selected as the data source to extract the range of surface water in the coal mining subsidence area [30].The mining time of the study area No. 1031 was from October 2013 to March 2015, and the surface water basin was unaffected by other working surface mining as of October 2016. Thus, the remote sensing images obtained from Landsat-8 from 2013 to 2016 were selected. The remote sensing images generated by a sensor exhibits radiation distortion and geometric distortion because of the change in the satellite operation state, aerosol refraction, random noise in the image generation process, and other factors. This condition affects the quality and application of the image. The preprocessing of remote sensing images, with steps including radiation calibration, atmospheric correction, geometric correction, and cropping, is required to eliminate these errors.(1)
Radiation calibration: Radiation calibration is the process of converting the digital quantitative value (DN) of the image into the radiation luminance value or physical quantities, such as reflectivity or surface temperature. The universal radiation calibration tool (Radiometric Calibration) in ENVI software is used to read the radiation calibration parameters in the landsat-8 remote sensing image metadata and complete the radiation calibration by the radiation luminance method.(2)
Atmospheric correction: Atmospheric correction eliminates the radiation error caused by the atmospheric influence and reverses the real surface reflectivity of the Earth. Initially, the radiation target data is input into the FLAASH module in ENVI to set the sensor-related parameters. Then, the corresponding atmospheric model is selected according to the latitude and imaging time, where the tropical (T) model was selected from July 2014 to September 2014 and the midlatitude summer (MLS) model for the rest of the months. According to the environment of the study area, the rural (R) mode was selected for the aerosol model, and the 2-band (K-T) was selected for the aerosol inversion method. Finally, the remaining relevant parameters were set, and the FLAASH module was run to determine the atmospheric correction results.(3)
Geometry correction: Geometry correction uses ground control points and geometric correction mathematical models to correct the errors caused by nonsystematic factors. The coordinate system is endowed to the image data during the correction process, which includes geographic coding. This step uses the geometry correction module in ENVI. The benchmark image is input, and the images corrected more than three image points of the same name are manually selected. Then, the automatic extraction function extracts more than 50 image points with the same image. The same name image point with the root mean square greater than 1 is eliminated, and the geometric correction module output file is finally run.(4)
Cropping: The appropriate size of the image output is cut according to the position of the study area to obtain the final image.Secondly, the reflection characteristics of water to different bands are determined by the water body index method. The principle of the method is to place the band with the strongest reflectivity and the weakest reflectance in the molecule and denominator, respectively. Then, the characteristics of the ratio operation are enlarged, while the characteristics of other objects are inhibited. Finally, water extraction is achieved [31]. At present, the more commonly used water index is normalized difference water index (NDWI) [32], and the improved normalized differential water index [31] is modified normalized difference water index (MNDWI).NDWI takes advantage of the high reflectivity in visible bands and low reflectivity in near-infrared and the range of midinfrared bands to highlight water information using ratio operation. In addition, attributing to the strong reflectivity of the vegetation in the near-infrared band, the water information in the green light band and the near-infrared band is used to calculate the simultaneously suppressed vegetation information, as shown by formula (18).(18)NDWI=Green−NIRGreen+NIR,where Green refers to the green light band reflectivity, and NIR refers to the near-infrared band reflectivity.The remote sensing images are processed using the MNDWI index to avoid miscalculation and improve the accuracy of water extraction. According to the characteristics of water in the image, the water boundary was visually checked. After determining the accuracy of the extraction results, the threshold was set, the water area was calculated, and the water boundary was converted to a vector file output [32]. Soil and buildings have similar spectral characteristics in the green and near-infrared bands. Thus, the mid-near-infrared band (NIR) in formula (18) is replaced by the midinfrared band (MIR) to obtain MNDWI. This method avoids the misjudgment of soil and buildings as water bodies when calculating the normalized water index band. MNDWI is obtained as follows:(19)MNDWI=Green−MIRGreen+MIR,where Green refers to the green optical band reflectivity, and MIR refers to the midinfrared band reflectivity.Considering the 1031 working face of Huaibei Wugou Coal Mine to be the research object, through remote sensing monitoring and on-site measurement, the evolution of the water accumulation area in the coal mining subsidence area with high phreatic level is divided into four periods to facilitate the study of the factors affecting the evolution process. The four stages include the unformed period, the synchronous growth period, the residual growth period, and the relatively stable period, as shown in Figure8.(1)
Unformed period: This stage is the beginning of the working surface mining to the surface water formation period, and in this period, the surface water does not form.(2)
Synchronous growth period: In this stage, the water on the surface stopped flowing to the working surface for the first time and expanded rapidly along the propulsion direction of the working surface, and the water area on the water surface grew rapidly.(3)
Residual growth period: In this stage, the cessation of mining took approximately five months from the working surface after the shutdown. The water slowly expanded along the propulsion direction of the working surface, slightly increasing the water area.(4)
Relative stability period: In this stage, approximately five months after the working surface stops mining, the boundary of water no longer expands along the mining direction, and the overall water area becomes stable [33].Figure 8
Evolution of the water range of 1031 working area.


### 5.2. Combining the Actual Measurement and the Model

Figure9 shows the comparison of the predicted water accumulation area and the measured water accumulation area of the inland water accumulation area from 2014 to 2017.Figure 9
Comparison of the evolution of the measured and simulated water accumulation area.


#### 5.2.1. Unformed Period

According to the image, the measured data did not include water observation in the study area before March 2014, the model forecast water range was approximately 3000 m2, and the water area in the two observation periods was small. Combined with the measured data, the surface of the mining area is mainly farmland, and the surface runoff in the subsidence basin is affected by soil leakage, evaporation, and other factors. In addition, the remote sensing images could not show the water range. The expected result is consistent with the measured results because of the small range of predicted water accumulation, considering the actual situation.


#### 5.2.2. Synchronous Growth Period

Water was observed in the measured data for the first time in March 2014. The water area was 11,000 m2, the model expected the water range in March 2014 to be 12,615 m2, and the difference was approximately 1615 m2, with an error of approximately 10%. According to the prediction results, the range of water in February 2014 is similar to that in March 2014 and the area of water accumulation is 12,812m2, indicating that a certain range of water may have been already present in the subsidence basin before the water was initially observed.From April 2014 to March 2015, the predicted range of water coincided well with the measured value, and the difference between those did not exceed 10% of the water area. Among them, the predicted result from September to November 2014 was because of the precipitation in September, which was 280 mm, leading to the rapid increase in the water level. Although numerous studies have paid significant attention to rainfall-induced instability of multilayer slopes, the interface between the layers is generally considered to be “zero thickness,” and the layer transition zone between the layers is neglected [10]. From October 2014 to November 2014, the water area decreased, and this observation was consistent with the measured value. The results obtained in December 2014 and March 2015 were greater than the predicted value, with a difference of approximately 1. At 5000 m2, the meteorological data of Huaibei City indicate that the winter temperature is less than 0°C, and the low temperatures result in a reduction in actual evaporation. Thus, the observation value is less than that predicted by the dynamic prediction model.


#### 5.2.3. Residual Growth Period

The measured data in April 2015 and October 2015 were 15,000 and 17,400 m2 and 157,839 and 181,655 m2. The prediction error did not exceed 10% of the measured value. The water range in July 2015 was from 171,000 m2 to 191,8941 m2. The difference accounted for approximately 11% of the meteorological data of the research area. During the July 2015 period, no precipitation occurred, and the highest daily temperature exceeded 30°C. In addition, the evaporation amount was large. As a result, the measured value of water is in the range of the model predictive value with a very low error.


#### 5.2.4. Relative Stability Period

After October 2015, the surface of the coal mining subsidence area became stable. Figure9 shows that the range of water accumulation in this stage was generally stable, and the predicted value has the same trend as the measured value. Considering the increase in human activity after the stability of the surface subsidence basin, a large difference was found between the predicted value and the measured value. For example, from March 2016 to June 2016, the measured area of the basin is greater than the predicted value, with an approximate difference of 20,000 m2. The water from the surrounding farmland irrigation converges into the water area affected by topography and permeability, and the new water supply increases. Thus, the measured value is greater than the model predicted value.In conclusion, the 3D spatial information dynamic prediction method predicts the evolution of the water area, and the water range predicted by the dynamic prediction model of the high diving position coal mining subsidence area is consistent with the measured results. In the period, when stagnant water is not formed, the predicted value of the model is basically consistent with the measured value, the error between the predicted value of the model, and the measured value in the synchronous growth period and the residual growth period is about 10%. Since the empirical data of groundwater exchange and ponding water leakage are all used in this study, the dynamic changes of groundwater depth and ponding leakage are not considered, causing some systematic errors that are difficult to eliminate in the model. In addition, there are still some inevitable errors in the method because of the lack of research data. During the relatively stable period, the prediction results of the model are consistent with the change trend of the measured ponding range. Thus, the model accuracy is high and can be used to predict the dynamic evolution of water in high diving mining subsidence areas.


## 5.1. Remote Sensing Monitoring and Field Measurement

The surface water in the mining area is spread on long range and for a long period of time, and the location is inconvenient for transportation. Conventional observation methods cannot easily obtain large-scale, long-time series, complete, and unified evolution data of the water range, whereas remote sensing observation methods have short observation periods and large amounts of data and high accuracy. Therefore, remote sensing data are selected as the data source to extract the range of surface water in the coal mining subsidence area [30].The mining time of the study area No. 1031 was from October 2013 to March 2015, and the surface water basin was unaffected by other working surface mining as of October 2016. Thus, the remote sensing images obtained from Landsat-8 from 2013 to 2016 were selected. The remote sensing images generated by a sensor exhibits radiation distortion and geometric distortion because of the change in the satellite operation state, aerosol refraction, random noise in the image generation process, and other factors. This condition affects the quality and application of the image. The preprocessing of remote sensing images, with steps including radiation calibration, atmospheric correction, geometric correction, and cropping, is required to eliminate these errors.(1)
Radiation calibration: Radiation calibration is the process of converting the digital quantitative value (DN) of the image into the radiation luminance value or physical quantities, such as reflectivity or surface temperature. The universal radiation calibration tool (Radiometric Calibration) in ENVI software is used to read the radiation calibration parameters in the landsat-8 remote sensing image metadata and complete the radiation calibration by the radiation luminance method.(2)
Atmospheric correction: Atmospheric correction eliminates the radiation error caused by the atmospheric influence and reverses the real surface reflectivity of the Earth. Initially, the radiation target data is input into the FLAASH module in ENVI to set the sensor-related parameters. Then, the corresponding atmospheric model is selected according to the latitude and imaging time, where the tropical (T) model was selected from July 2014 to September 2014 and the midlatitude summer (MLS) model for the rest of the months. According to the environment of the study area, the rural (R) mode was selected for the aerosol model, and the 2-band (K-T) was selected for the aerosol inversion method. Finally, the remaining relevant parameters were set, and the FLAASH module was run to determine the atmospheric correction results.(3)
Geometry correction: Geometry correction uses ground control points and geometric correction mathematical models to correct the errors caused by nonsystematic factors. The coordinate system is endowed to the image data during the correction process, which includes geographic coding. This step uses the geometry correction module in ENVI. The benchmark image is input, and the images corrected more than three image points of the same name are manually selected. Then, the automatic extraction function extracts more than 50 image points with the same image. The same name image point with the root mean square greater than 1 is eliminated, and the geometric correction module output file is finally run.(4)
Cropping: The appropriate size of the image output is cut according to the position of the study area to obtain the final image.Secondly, the reflection characteristics of water to different bands are determined by the water body index method. The principle of the method is to place the band with the strongest reflectivity and the weakest reflectance in the molecule and denominator, respectively. Then, the characteristics of the ratio operation are enlarged, while the characteristics of other objects are inhibited. Finally, water extraction is achieved [31]. At present, the more commonly used water index is normalized difference water index (NDWI) [32], and the improved normalized differential water index [31] is modified normalized difference water index (MNDWI).NDWI takes advantage of the high reflectivity in visible bands and low reflectivity in near-infrared and the range of midinfrared bands to highlight water information using ratio operation. In addition, attributing to the strong reflectivity of the vegetation in the near-infrared band, the water information in the green light band and the near-infrared band is used to calculate the simultaneously suppressed vegetation information, as shown by formula (18).(18)NDWI=Green−NIRGreen+NIR,where Green refers to the green light band reflectivity, and NIR refers to the near-infrared band reflectivity.The remote sensing images are processed using the MNDWI index to avoid miscalculation and improve the accuracy of water extraction. According to the characteristics of water in the image, the water boundary was visually checked. After determining the accuracy of the extraction results, the threshold was set, the water area was calculated, and the water boundary was converted to a vector file output [32]. Soil and buildings have similar spectral characteristics in the green and near-infrared bands. Thus, the mid-near-infrared band (NIR) in formula (18) is replaced by the midinfrared band (MIR) to obtain MNDWI. This method avoids the misjudgment of soil and buildings as water bodies when calculating the normalized water index band. MNDWI is obtained as follows:(19)MNDWI=Green−MIRGreen+MIR,where Green refers to the green optical band reflectivity, and MIR refers to the midinfrared band reflectivity.Considering the 1031 working face of Huaibei Wugou Coal Mine to be the research object, through remote sensing monitoring and on-site measurement, the evolution of the water accumulation area in the coal mining subsidence area with high phreatic level is divided into four periods to facilitate the study of the factors affecting the evolution process. The four stages include the unformed period, the synchronous growth period, the residual growth period, and the relatively stable period, as shown in Figure8.(1)
Unformed period: This stage is the beginning of the working surface mining to the surface water formation period, and in this period, the surface water does not form.(2)
Synchronous growth period: In this stage, the water on the surface stopped flowing to the working surface for the first time and expanded rapidly along the propulsion direction of the working surface, and the water area on the water surface grew rapidly.(3)
Residual growth period: In this stage, the cessation of mining took approximately five months from the working surface after the shutdown. The water slowly expanded along the propulsion direction of the working surface, slightly increasing the water area.(4)
Relative stability period: In this stage, approximately five months after the working surface stops mining, the boundary of water no longer expands along the mining direction, and the overall water area becomes stable [33].Figure 8
Evolution of the water range of 1031 working area.


## 5.2. Combining the Actual Measurement and the Model

Figure9 shows the comparison of the predicted water accumulation area and the measured water accumulation area of the inland water accumulation area from 2014 to 2017.Figure 9
Comparison of the evolution of the measured and simulated water accumulation area.


### 5.2.1. Unformed Period

According to the image, the measured data did not include water observation in the study area before March 2014, the model forecast water range was approximately 3000 m2, and the water area in the two observation periods was small. Combined with the measured data, the surface of the mining area is mainly farmland, and the surface runoff in the subsidence basin is affected by soil leakage, evaporation, and other factors. In addition, the remote sensing images could not show the water range. The expected result is consistent with the measured results because of the small range of predicted water accumulation, considering the actual situation.


### 5.2.2. Synchronous Growth Period

Water was observed in the measured data for the first time in March 2014. The water area was 11,000 m2, the model expected the water range in March 2014 to be 12,615 m2, and the difference was approximately 1615 m2, with an error of approximately 10%. According to the prediction results, the range of water in February 2014 is similar to that in March 2014 and the area of water accumulation is 12,812m2, indicating that a certain range of water may have been already present in the subsidence basin before the water was initially observed.From April 2014 to March 2015, the predicted range of water coincided well with the measured value, and the difference between those did not exceed 10% of the water area. Among them, the predicted result from September to November 2014 was because of the precipitation in September, which was 280 mm, leading to the rapid increase in the water level. Although numerous studies have paid significant attention to rainfall-induced instability of multilayer slopes, the interface between the layers is generally considered to be “zero thickness,” and the layer transition zone between the layers is neglected [10]. From October 2014 to November 2014, the water area decreased, and this observation was consistent with the measured value. The results obtained in December 2014 and March 2015 were greater than the predicted value, with a difference of approximately 1. At 5000 m2, the meteorological data of Huaibei City indicate that the winter temperature is less than 0°C, and the low temperatures result in a reduction in actual evaporation. Thus, the observation value is less than that predicted by the dynamic prediction model.


### 5.2.3. Residual Growth Period

The measured data in April 2015 and October 2015 were 15,000 and 17,400 m2 and 157,839 and 181,655 m2. The prediction error did not exceed 10% of the measured value. The water range in July 2015 was from 171,000 m2 to 191,8941 m2. The difference accounted for approximately 11% of the meteorological data of the research area. During the July 2015 period, no precipitation occurred, and the highest daily temperature exceeded 30°C. In addition, the evaporation amount was large. As a result, the measured value of water is in the range of the model predictive value with a very low error.


### 5.2.4. Relative Stability Period

After October 2015, the surface of the coal mining subsidence area became stable. Figure9 shows that the range of water accumulation in this stage was generally stable, and the predicted value has the same trend as the measured value. Considering the increase in human activity after the stability of the surface subsidence basin, a large difference was found between the predicted value and the measured value. For example, from March 2016 to June 2016, the measured area of the basin is greater than the predicted value, with an approximate difference of 20,000 m2. The water from the surrounding farmland irrigation converges into the water area affected by topography and permeability, and the new water supply increases. Thus, the measured value is greater than the model predicted value.In conclusion, the 3D spatial information dynamic prediction method predicts the evolution of the water area, and the water range predicted by the dynamic prediction model of the high diving position coal mining subsidence area is consistent with the measured results. In the period, when stagnant water is not formed, the predicted value of the model is basically consistent with the measured value, the error between the predicted value of the model, and the measured value in the synchronous growth period and the residual growth period is about 10%. Since the empirical data of groundwater exchange and ponding water leakage are all used in this study, the dynamic changes of groundwater depth and ponding leakage are not considered, causing some systematic errors that are difficult to eliminate in the model. In addition, there are still some inevitable errors in the method because of the lack of research data. During the relatively stable period, the prediction results of the model are consistent with the change trend of the measured ponding range. Thus, the model accuracy is high and can be used to predict the dynamic evolution of water in high diving mining subsidence areas.


## 5.2.1. Unformed Period

According to the image, the measured data did not include water observation in the study area before March 2014, the model forecast water range was approximately 3000 m2, and the water area in the two observation periods was small. Combined with the measured data, the surface of the mining area is mainly farmland, and the surface runoff in the subsidence basin is affected by soil leakage, evaporation, and other factors. In addition, the remote sensing images could not show the water range. The expected result is consistent with the measured results because of the small range of predicted water accumulation, considering the actual situation.


## 5.2.2. Synchronous Growth Period

Water was observed in the measured data for the first time in March 2014. The water area was 11,000 m2, the model expected the water range in March 2014 to be 12,615 m2, and the difference was approximately 1615 m2, with an error of approximately 10%. According to the prediction results, the range of water in February 2014 is similar to that in March 2014 and the area of water accumulation is 12,812m2, indicating that a certain range of water may have been already present in the subsidence basin before the water was initially observed.From April 2014 to March 2015, the predicted range of water coincided well with the measured value, and the difference between those did not exceed 10% of the water area. Among them, the predicted result from September to November 2014 was because of the precipitation in September, which was 280 mm, leading to the rapid increase in the water level. Although numerous studies have paid significant attention to rainfall-induced instability of multilayer slopes, the interface between the layers is generally considered to be “zero thickness,” and the layer transition zone between the layers is neglected [10]. From October 2014 to November 2014, the water area decreased, and this observation was consistent with the measured value. The results obtained in December 2014 and March 2015 were greater than the predicted value, with a difference of approximately 1. At 5000 m2, the meteorological data of Huaibei City indicate that the winter temperature is less than 0°C, and the low temperatures result in a reduction in actual evaporation. Thus, the observation value is less than that predicted by the dynamic prediction model.


## 5.2.3. Residual Growth Period

The measured data in April 2015 and October 2015 were 15,000 and 17,400 m2 and 157,839 and 181,655 m2. The prediction error did not exceed 10% of the measured value. The water range in July 2015 was from 171,000 m2 to 191,8941 m2. The difference accounted for approximately 11% of the meteorological data of the research area. During the July 2015 period, no precipitation occurred, and the highest daily temperature exceeded 30°C. In addition, the evaporation amount was large. As a result, the measured value of water is in the range of the model predictive value with a very low error.


## 5.2.4. Relative Stability Period

After October 2015, the surface of the coal mining subsidence area became stable. Figure9 shows that the range of water accumulation in this stage was generally stable, and the predicted value has the same trend as the measured value. Considering the increase in human activity after the stability of the surface subsidence basin, a large difference was found between the predicted value and the measured value. For example, from March 2016 to June 2016, the measured area of the basin is greater than the predicted value, with an approximate difference of 20,000 m2. The water from the surrounding farmland irrigation converges into the water area affected by topography and permeability, and the new water supply increases. Thus, the measured value is greater than the model predicted value.In conclusion, the 3D spatial information dynamic prediction method predicts the evolution of the water area, and the water range predicted by the dynamic prediction model of the high diving position coal mining subsidence area is consistent with the measured results. In the period, when stagnant water is not formed, the predicted value of the model is basically consistent with the measured value, the error between the predicted value of the model, and the measured value in the synchronous growth period and the residual growth period is about 10%. Since the empirical data of groundwater exchange and ponding water leakage are all used in this study, the dynamic changes of groundwater depth and ponding leakage are not considered, causing some systematic errors that are difficult to eliminate in the model. In addition, there are still some inevitable errors in the method because of the lack of research data. During the relatively stable period, the prediction results of the model are consistent with the change trend of the measured ponding range. Thus, the model accuracy is high and can be used to predict the dynamic evolution of water in high diving mining subsidence areas.


## 6. Conclusion

(1)
In view of addressing the problem of water accumulation in the surface subsidence basin, a dynamic prediction model of the water accumulation area in the high-water-level coal mining subsidence area was constructed by combining the dynamic prediction of the coal mining subsidence area based on the probability integral method of the Knothe time function and the principle of water balance in coal mining subsidence water area. Dynamically predicting the depth, area, volume, maximum storage capacity, and scope of the basin became possible with the advancement of the working face and the passage of time. The prediction results have certain timeliness and can effectively avoid the failure of real-time follow-up or repeated management of the comprehensive treatment of subsidence basins because of the lack of timeliness of the data.(2)
Taking the 1031 working face of Wugou Mine as an example, the applicability of the dynamic prediction model of the water accumulation area in the high-water-level coal mining subsidence area was verified. The error between the predicted value of the model and the measured value is approximately 10%. This newly developed model showed high accuracy and thus can be used for predicting the dynamic evolution of water range in high diving mining subsidence areas.(3)
In the dynamic prediction of the 3D spatial information of the coal mining subsidence water area, the acquired hydraulic information in the coal mining subsidence water area basin can provide a reference for the comprehensive management of the coal mining subsidence area and can also assist in protecting the ecology and surrounding geography. Environment and mining subsidence dynamic prediction results can serve as a guide for mining area reclamation. The predicted water range and water depth can provide data for the establishment of lakes and wetland ecosystem in the coal mining subsidence area. A water storage project can be established using the subsidence basin volume. A widely applicable dynamic prediction model of surface subsidence was constructed in the mining areas with high phreatic levels, and reliable technical support is provided for land reclamation and ecological planning in similar mining areas.(4)
The dynamic prediction model of the water accumulation area in the high-water-level coal mining subsidence area needs to be further improved. The water accumulation area in the subsidence area of high-water-water mining established in this study only considers the flat surface and is only affected by the mining of a single working face. The water volume of the recharge term was calculated, and the mutual conversion between these recharge terms and excretion terms was not analyzed in depth. In addition, the influencing factors, such as human activities and surrounding water systems, are not taken into account in the model. Therefore, this model still needs to be further improved, accounting for many influencing factors to accurately predict the surface water range in the high-water-level coal mining subsidence area with a larger area, undulating terrain, and complex hydrogeological conditions.



---

*Source: 1568050-2022-03-01.xml*