Spatio-Temporal Evolution Law of Surface Subsidence Basin with Insufﬁcient Exploitation of Deep Coal Resources in Aeolian Sand Area of Western China

: Coal is one of the fundamental fossil energy supporting the world’s economy. The synergistic development between efﬁcient coal mining and ecological environment protection is the inevitable requirement for the preservation of global harmony. As the world’s largest coal producer, China has conducted a strategic shift from east to west in terms of the exploitation of its energy resources, posing a serious threat to the fragile ecological environment of the western region. In particular, the surface subsidence caused by coal mining is the root of the ecological deteriotation and the destruction of ground structures. However, it is difﬁcult to reveal the law of large-scale surface subsidence in western mining areas merely by conventional measurement methods such as leveling, on account of the high intensity of coal seam mining, the weakness of the lithology of overlying rock and the large thickness of wind-blown sand strata. In view of this, small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) technology was used in this study to obtain the time series of surface vertical displacement during the whole mining process of the 2401 working face in the Yingpanhao coal mine, Inner Mongolia. Based on the deformation data, the dynamic evolution characteristics of surface subsidence under high intensity mining in the western mining area were analyzed exhaustively. It was found that the surface subsidence is characterized by an extensive coverage range (48.52 km 2 ) with minimal ground settlement (250 mm) in the study area. Meanwhile, the boundary shape of the subsidence basin followed a “circular-parallelogram-trapezoid” changeable process and the coverage area of the basin experienced three stages: a linear increasing period, a temporary stagnation period, and a re-expansion period. Furthermore, there existed an abnormal uplift phenomenon on the east side of the open-off cut in the 2401 working face. Combined with the structure of overlying strata, this paper carried out a preliminary analysis on the reasons of the abovementioned phenomenon. The research results are of vital realistic signiﬁcance for ground buildings and ecological environmental protection in the aeolian sand mining area in Western China.


Introduction
Coal is one of the fundamental fossil energy sources supporting the world's economy. The synergistic development of efficient coal mining and ecological environment protection is the inevitable requirement for the preservation of harmony. As the world's largest coal producer, China has strategically shifted its intensive mining activities from east to west as a result of the gradual depletion of coal resources in the eastern region. However, largescale and high-intensity coal exploitation poses a serious threat to the fragile ecological environment in the western aeolian sand area. In particular, the surface subsidence caused by coal mining would destroy the original equilibrium structure of underground aquifer, threaten the safety of infrastructure including railways, highways, as well as natural gas transmission pipelines, and weaken the carrying capacity of the ecological environment. Therefore, it has become a hot research topic to master the mining subsidence law of the aeolian sand area in Western China, in order to solve scientifically the contradiction between coal resources exploitation and the surface human living environment.
Early studies have shown that drastic strata movement and discontinuous surface damage would occur rapidly after coal mining in western aeolian sand areas, such as sudden surface collapse [1], permanent ground fissures [2][3][4][5][6], the decline of groundwater level accompanied by the above severe destruction [7], as well as secondary damages such as the loss of soil nutrients and the reduction of microbial diversity [8]. According to the analysis of Yan et al. [9][10][11], the above severe failure phenomenon has an apparent relationship with the characteristics of the strata structure featured as a large coal seam thickness, shallow burial depth, and thin overlying bedrock. Under the condition of highintensity mining, the overlaying strata in these areas are prone to sliding instability, which causes synchronous failure of the aeolian sand layer. Consequently cracks connecting mining fields and the surface are formed. However, the conclusions of the above research are only applicable to describe the surface damage characteristics of shallow coal seams.
The latest research [12] indicates that the severity of surface subsidence is greatly reduced along with the increased mining depth. Moreover, this reduction trend is particularly obvious for the first working face in the state of insufficient mining. To understand this phenomenon, preliminary studies [12,13] on the overburden deformation characteristics of deep mining in Western China were undertaken with similar material simulation experiments and numerical simulation experiments. It was considered that the giant thick Cretaceous strata weakly cemented to the stratum located above the deeply buried coal seam has superior resistance to deformation and damage because of its good integrity and large thickness [14][15][16][17], thus controlling the development of the overburden movement to the surface. However, these investigations did not specifically discuss the ground settlement characteristics after deep coal mining due to the lack of comprehensive surface deformation observation data.
With the rapid development of satellite-based earth observation technology, a new generation of surface deformation monitoring technology represented by differential interferometry synthetic aperture radar (D-InSAR) [18][19][20], persistent scatterers InSAR (PS-InSAR) [21] and small baseline subset InSAR (SBAS-InSAR) [22][23][24] has been proposed. With great advantages such as a wide monitoring range, a high observation density, historical backtracking, a high accuracy of deformation monitoring and a low influence by the weather, these techniques have been widely used in the fields of volcanoes [25][26][27][28], earthquakes [29][30][31][32][33][34], landslides [35][36][37][38], and other geological hazards [39][40][41]. At the same time, it has been widely used in the field of surface subsidence monitoring in mining areas, including the analysis of the relationship between surface deformation and underground mining [42], the evaluation of potential mine geological hazards by calculating the prediction parameters of surface subsidence [43], and the analysis of the influence of additional factors, such as groundwater migration [44,45] and underground coal fire [46], on surface displacement in mining areas. Nevertheless, there are few cases that employ time series InSAR (TS-InSAR) technology to analyze mining-induced surface subsidence in the aeolian sand area of Western China.
The SBAS-InSAR technique adopts the method of acquiring interferometric image pairs with short spatial and temporal baselines from multiple master images, which can use as many interferograms as possible in the same observation time to participate in the deformation sequence solution. Compared with the PS-InSAR technology, it is more suitable for the aeolian sand area in Western China without permanent scatterers such as houses, roads, railways, bridges, and exposed rocks. Therefore, the time series of surface subsidence during the advancement of the 2401 working face in the Yingpanhao coal mine, located in the central part of the Ordos Basin, China, were obtained by the SBAS-InSAR technique, and the relationship between the two geometric features (shape and coverage area) of the surface subsidence basin and the advancing length of the working face were analyzed in this study. Moreover, the influence mechanism of the aeolian sand layer and the weakly cemented covered rock combination structure on surface subsidence was further analyzed. The purpose of this study was to grasp the characteristics of surface subsidence under the insufficient condition of deep coal mining in the aeolian sand area of Western China so as to provide insights for analyzing the movement mechanism of giant thick weakly cemented overburden.

Study Area
As shown in Figure 1, the study area is located in the inland of Northwest China with a flat terrain and sandy landform, which has a typical midtemperate semiarid continental climate with little rainfall and high surface evaporation, resulting in an extreme shortage of surface water resources. However, the groundwater resources in this area are relatively abundant, especially in shallow strata. Monsoons are prevalent in the study area with an average wind speed of 13.4 m/s over the years, especially in the northwest monsoon in spring. The extreme maximum wind speed can reach 19.0 m/s, allowing dust storms to form frequently. The study area is located in the central part of the Ordos platform, which is the most complete and stable structural unit existing in China. It provides advantageous conditions for the formation of a relatively complete and simple stratigraphic structure and the avoidance of large fault development. According to drilling data, the strata in the area from bottom to top are: Upper Triassic Yanchang Formation (T 3 y), Middle-Lower Jurassic Yan'an Formation (J 1-2 y), Middle Jurassic Zhiluo Formation (J 2 z), Anding Formation (J 2 a), Lower Cretaceous Zhidan Group (K 1 zh), and Quaternary (Q 4 ). Figure 2 describes the stratigraphic distribution of the study area in detail. Overall, the lithology of the strata is mostly sandstone. Such lithologies have poor mechanical properties due to weak cementation, but have a large thickness. Specifically, the high-lying Cretaceous Zhidan Group formations have compressive strengths between 10-20 Mpa, tensile strengths between 0.29-1.61 Mpa, a modulus of elasticity between 1.76-2.91 Gpa, and a cohesion between 4.8-8.6 Mpa with an average thickness of more than 300 m. The 2-2 coal seam currently being mined in the Yingpanhao coal mine was formed during the Jurassic period and is located in the third section of J 1-2 y, with an average burial depth of 722.88 m, a thickness of 6.41 m, and a dip angle of 1 • . The study area is scarce of surface water, but is abundant in groundwater resources, which can meet the water demand of local residents for living and production. According to the borehole pumping data of the Yingpanhao coal mine, there are mainly four aquifers in the overburden of the 2-2 coal seam, as shown in Figure 2. From the surface downward, they are: Q 4 pore phreatic aquifer PW, K 1 zh confined aquifer CW1, Middle Jurassic (J 2 ) confined aquifer CW2, and J 2 z bottom to 2-2 coal seam confined aquifer CW3. The main information of each aquifer is shown in Table 1. The unit water inflow q and permeability coefficient K of each aquifer were determined by steady flow single-hole tests (drop depth of 10 m, steady pumping time of 8 h, and caliber of 91 mm). The phreatic aquifer PW is mainly supplied by atmospheric precipitation and lateral runoff in adjacent areas. At the same time, the phreatic aquifer has a certain hydraulic connection with the lower confined aquifer CW1. The confined aquifer CW1 has a high porosity, strong water permeability and conductivity, which forms a good storage space for groundwater and is the main aquifer in this area. The hydraulic connection between the CW1 and the lower confined water aquifer is small, and further down, the confined aquifers CW2 and CW3 are weakly water-rich with poor permeability and hydraulic conductivity. In the formation J 2 a, J 2 z, and J 1-2 y, there are multiple layers of mudstone, sandy mudstone, and siltstone with thin thickness, dense lithology, and poor permeability. Consequently, they can be regarded as the relatively impermeable layers. The 2401 working face is located in the northeast of the industrial square. The average depth of the working face is about 742 m. The average mining thickness is about 5.6 m. The mining width is 300 m. The working face is mined by comprehensive mechanization with a longwall arrangement and full height of one mining. The mining sequence is backward and the roof management method is a caving method. The 2401 working face was advanced from east to west in early December 2019 from the open cut. The advancing length was 1022 m until 31 January 2021. There is no other working face adjacent to the 2401 working face, except for the 2202 working face in the second mining section at a distance of 2027 m to its south. Since they are far away from each other, it can be initially considered that the two working faces do not produce obvious mining disturbance impact on each other's subsidence area. Figure 3 is a layout of the 2401 working face.

Materials and Methods
In this section, the SAR images employed in this paper, the basic principles of the SBAS-InSAR technology, and the main data processing workflow are presented.

SAR Images
The SAR image dataset used in this study was obtained by C-band high resolution synthetic aperture radar on the Sentinel-1 A satellite of ESA in wide swath mode, spanning the period from 29 July 2019 to 19 May 2021 with a total of 53 scenes. The coverage of the Sentinel-1 A SAR images used is shown in the blue rectangle in Figure 1b, and the main parameter is shown in Table 2. The study area is located in the eastern part of the original SLC image coverage. In order to minimize the computation time and storage space required for the SBAS-InSAR processing, all images were first uniformly cropped in this study. The cropped SLC image coverage is shown in the red rectangular in Figure 1d. In addition, ASTER GDEM V2 data with a spatial resolution of 30 m were also used in this study to assist in the interferogram alignment and remove the topographic phase in the interferogram. At the same time, the corresponding precise orbit data of each SAR image were used to weaken the orbit error in the interference phase.

SBAS-InSAR Technology
Several subsets were generated by the free combination of SAR images covering the study area through the setting of temporal baseline and spatial baseline thresholds in the SBAS-InSAR technique [22]. The baselines were smaller within the same subset and larger between subsets. The terrain phase simulation and removal were performed using external DEM, and then the selection of high-quality coherence points was based on the coherence coefficients of the interferogram. Finally, the least-squares solution (LS) within a subset and the singular value decomposition (SVD) between subsets were obtained on the interference phases of these highly coherent points to obtain high-precision deformation time series and annual average deformation rates. The basic principle of SBAS-InSAR technology is described below.
One super master image is selected from N + 1 SAR images covering the same area, and then other SAR images are registered with it. Then, M multiview differential interferograms are generated.
Assume that the images at moments t A and t B generate the jth interferogram. The interferometric phase of the pixel with azimuthal coordinate x and distance coordinate r in the map is: where j (1, 2, · · · , M); λ is the central wavelength of the radar signal; d(t A , x, r) and d(t B , x, r) are the cumulative deformations at moments t A and t B (t B > t A ), respectively, with respect to the radar line-of-sight direction at moment t 0 ; ∆φ j topo (x, r) represents the terrain phase of the residual in the differential interferogram. If the DEM introduced in the differential interferogram has high accuracy, most of the terrain phase can be removed, then the residual terrain information in the differential interferogram is small and can be ignored in the calculation process; ∆φ j APS (x, r, t B , t A ) represents the atmospheric delay phase; ∆φ j noise (x, r) represents the coherent noise (such as the system thermal noise). If the influence of atmospheric delay phase, residual terrain phase and noise phase are not considered, Equation (3) can be simplified as: In order to obtain the deformation sequence with physical significance, the average phase velocity of two adjacent images is set as: Thus, the following matrix equation can be obtained: In the formula, D is an M × N matrix, where each row corresponds to an interferogram and each column corresponds to a SAR image obtained at a certain time. The value of the main image column is +1; the value of the auxiliary image column is −1; and values of the remaining images are 0.
Since the SBAS-InSAR technology adopts the multimain image strategy, matrix D is generally a rank-deficient matrix. The generalized inverse matrix of D is obtained by the SVD method, and the minimum norm solution of the velocity vector is obtained. Finally, the deformation of each time period can be obtained by integrating the velocity of each time period. SBAS-InSAR not only inherits the advantages of traditional D-InSAR, but also improves the time resolution of deformation monitoring, effectively overcoming the decoherence problem caused by a long time baseline and spatial baseline in surface deformation monitoring.

Data Processing
In this study, 53 Sentinel-1 A SAR images were processed by the SBAS-InSAR module in SARScape. The workflow of data processing is shown in Figure 4.
First, the temporal baseline and spatial baseline thresholds were set to 60 days and 2% of the critical limit to generate as many connected image pairs with a good coherence as possible. The generated spatio-temporal baseline connection maps are shown in Figure 5.
Then, the image pairs were registered before interference processing. After interference calculation, the terrain phase information was eliminated in combination with ASTER GDEM V2. In order to reduce the computational effort and suppress the noise, the interferograms were multilook-processed with an azimuthal direction of 1 and a range direction of 4. Then, the filtering process was performed by the Goldsteid method to further eliminate the effect of noise in the multilook differential interferogram. The minimum-cost flow method was used for unwrapping, and the threshold of the unwrapping coherence coefficient was set to 0.3. After checking the coherence map, the image pairs with poor coherence were removed. Some of the coherence and interferograms maps of the remaining pairs are shown in Figures 6 and 7. The number and the date of each pair are shown in Table 3.    Table 3. The layout of the 2401 working face is marked by the red rectangle.
In Figure 6a-d, the coherence values in the area which are directly above the 2401 working face are lower than those in the surrounding area. It indicates that the surface subsidence basin is expanding during this period. The interferometric stripes in Figure 7a-d also prove that is indeed occurring at this stage. In Figure 6e,f, the coherence values of the area around the 2401 working face is significantly reduced. As can be seen from Table 3, these coherence images correspond to dates ranging from 6 April 2020 to 27 October 2020. Combined with the historical meteorological records of the study area, the distribution of annual precipitation in the area where the Yingpanhao mine is located is extremely uneven, concentrated in July, August, and September. Moreover, the frost-free period is short throughout the year with the annual freeze starting in October and thawing in April of the following year. Studies [47,48] have been conducted to show that rainfall can cause loss of coherence. Therefore, changes in surface water content due to precipitation and low-temperature permafrost are likely to be the main factors contributing to the reduction in coherence over this time scale. From the statistical distribution of the coherence values shown in Figure 8, it can be more intuitively seen that the coherence values of the pixels in   Table 3. The layout of the 2401 working face is marked by the red rectangle.  Table 3.
The high coherence shown in Figure 6i and the absence of interference fringes shown in Figure 7i indicate that no significant subsidence has occurred at the surface after mining stopped. The interference fringes and coherence variations shown in the interferograms and coherence maps described above are consistent with the mining activity on the 2401 working face.
Then, 51 geodetic control points (GCP) were selected for orbit refinement and reflattening in the region without residual topographic streaks, deformation streaks, and phase jumps. After that, two inversions were performed. The first inversion estimated the displacement rate and residual topography, and the second inversion performed phase unwrapping to optimize the results.
The research in [21] shows that the atmospheric delay phase is generally a lowfrequency signal in the spatial domain and a high-frequency signal in the time domain. The deformation phase shows low-frequency signals in both spatial and temporal domains; the noise phase is a high-frequency signal in the space and time domain. Therefore, the phase caused by the atmospheric delay can be obtained by high-pass filtering in time and low-pass filtering in space. The atmosphere high pass size and atmosphere low pass size were set to 365 days and 1600 m, respectively.
Finally, the line of sight (LOS) deformation was converted to a vertical deformation based on the formula: where W Los (x, y) is the deformation value at any point in the LOS direction and θ represents the angle of incidence.

Results
After processing, the surface subsidence time series during the mining of the 2401 working face in the Yingpanhao coal mine were finally obtained. Figure 9 shows the distribution of the surface subsidence rate in the study area. The maximum settlement rate is 180 mm/year. In this section, the subsidence data measured by SBAS are first verified using level survey data. Then, the dynamic evolution process of the subsidence basin, its shape, and the relationship between the area and the advancing length of the working face are described.

Verification of Settlement Monitoring Results
In order to verify the SBAS-InSAR settlement monitoring results, leveling measurement data were collected at the surface around the 2401 working face in this study. The leveling observation line was located at the northwest side of the 2401 working face, as shown in the Figure 9, which consisted of 22 measurement points (B1-B22) in total and was laid along the road. A total of four observations were conducted. The first observation date was 2 December 2019, and the remaining observation dates were 6 March 2020, 24 July 2020 and 8 March 2021, respectively.
For comparison, the cumulative subsidence values recorded at 26 November 2019, 1 March 2020, 23 July 2020, and 8 March 2021 were extracted from the SBAS subsidence monitoring sequences containing pixels with level measurement points B1-B22 in this study. The cumulative subsidence for the other dates was then calculated by taking 26 November 2019 as the first observation. The comparison results of SBAS-InSAR and level settlement monitoring are shown in Figure 10. It can be clearly seen that the subsidence curves composed of each point also show consistent subsidence characteristics. However, it should be noted that the monitoring results of B1-B13 in July 2020 and March 2021 are slightly different. Specifically, the monitoring results of the level survey are larger than those of SBAS-InSAR by about 17 mm. Surface deformation in mining areas is usually a threedimensional phenomenon [49]. Therefore, the discrepancy shown in Figure 10 is related to the fact that the equation for converting subsidence from the LOS direction to the vertical direction ignores subsidence perpendicular to the line-of-sight direction. However, since the differences are not significant, it can be assumed that the subsidence monitoring results of SBAS-InSAR have limited influence on the analysis of surface subsidence characteristics.

Dynamic Development Process of Surface Subsidence Basin
The vertical deformation time series of the surface in the study area from 29 July 2019 to 19 May 2021 were obtained by the above SBAS-InSAR process. The spatio-temporal evolution of the surface subsidence basin in the study area is shown in Figure 11 with a 2-month interval.
Before 26 November 2019, the 2401 working face had not been mined; meanwhile the ground surface did not show subsidence, as shown in Figure 11a,b. On 6 February 2020, the working face was advanced by a length of 100 m, while slight subsidence had occurred on the surface, as shown in Figure 11c. As the advancing length increased, the area where surface subsidence occurred expanded rapidly to the west, as shown in Figure 11d-k. During this period, the settlement amount also gradually increased. In particular, the maximum cumulative surface subsidence exceeded 200 mm on 31 January 2021, when the advancing length was approximately 1022 m. Afterwards, the mining was stopped and the extent and sink values of the surface subsidence basin increased slightly. Figure 11k shows the final formation of the surface subsidence basin. The area where surface subsidence occurs is basically located at the west side of the working face open cut. On the east side of the open cut, there is a small-scale uplift of the ground surface instead of subsidence. However, the uplift value is relatively small, not exceeding 50 mm at maximum. Such an anomalous phenomenon indicates that the movement of the overburden above the open cut is probably not in the form of a simple downward bending, but a rotation towards the center of the mining area. A further analysis of this anomaly is carried out in Section 5.1. According to the degree of sinkage, the surface settlement area can be divided into an edge zone and a central zone. The edge zone is located on the west side of the working face, and the sinkage value of each place is small, but the coverage is extremely extensive. The central zone is located directly above the working face, slightly on the side of the coal wall, and is axially distributed with the axis of orientation of the working face, which is consistent with the surface subsidence characteristics after longwall mining of a horizontal coal seam. Meanwhile, the local area where a slight uplift occurs on the east surface of the open cut is called the uplift zone.

Variations of Surface Subsidence Basin Shape
Shape is the most intuitive geometric feature of a surface subsidence basin, which can be delineated by boundary contour lines. In a previous study [50], considering the influence of measurement error on monitoring results, 10 mm was used as the determination value of subsidence basin boundary. It has been demonstrated that the accuracy of surface subsidence monitoring by the SBAS-InSAR technology can reach subcentimeter level [51]. Moreover, the smaller the deformation, the higher the accuracy of the monitoring results. Therefore, this study used the 10 mm contour of subsidence to delineate the influence range of surface mining above the 2401 working face in each period, as shown in the dark green contour C 10 in Figure 12. In addition, to distinguish the edge zone and the center zone of the subsidence basin, this study also determined the areas where the subsidence exceeded 30 mm and 50 mm, as shown by the blue contour C 30 and the black contour C 50 , respectively.
First of all, it can be obviously recognized that the contour C 10 is more cluttered, which is increasingly serious as time passes backwards. The contours C 30 and C 50 are relatively distinct. Secondly, all three contours experienced different degrees of change during the advancement of the 2401 working face from east to west.
Among them, the contour shape of C 10 changed most significantly, roughly from circular to parallelogram and then finally to trapezoid. Accordingly, the outward expansion of the surface subsidence basin was divided into three stages as follows.  The shapes of C 30 and C 50 were less variable, being circular first and then oval. Moreover, the variation process was obviously directly associated with the mining activity of the 2401 working face. Assuming that 30 mm of subsidence is exactly the threshold value to distinguish the edge zone and the center zone of the subsidence basin, then it can be considered that the area within C 30 is the center zone of the subsidence basin, and the area between C 10 and C 30 is the edge zone.
From the above analysis, it can be seen that the settlement in the central zone was directly related to the mining of the 2401 working face, while the settlement in the edge zone was not directly affected by the underground mining, but more likely induced by the decline of groundwater caused by mining. A further analysis of this phenomenon is carried out in Section 5.2.

Relationship between Coverage Area of Subsidence Basin and Advancing Length of Working Face
Since the continuous exploitation of the 2401 working face is the direct cause of the expansion of the surface subsidence basin, it is necessary to have an in-depth investigation of the quantitative relationship between the geometric characteristics of the subsidence basin and the length of the working face. Accordingly, the area of the circled areas of C 10 , C 30 , and C 50 curves, denoted as A 10 , A 30 , and A 50 respectively, was first counted in this study, as shown in the Table 4 and Figure 13.  As can be seen from Figure 13a, A 10 is much larger than A 30 and A 50 . When the mining length L of the 2401 working face reaches 1022.3 m, A 10 tends to be 50 km 2 , while both A 30 and A 50 are only within 10 km 2 . Specifically, A 10 is about 48.52 km 2 , A 30 is about 7.47 km 2 , and A 50 is about 2.45 km 2 . It indicates that the central zone was only one fifth of the entire subsidence basin area.
Combined with Figure 13a,b, it can be found that although the magnitudes of A 10 , A 30 , and A 50 are quite different from each other, all of them increase with L through three stages: a linear increasing period, a temporary stagnation period, and a re-expansion period. All of them have basically the same working face advance length corresponding to the stage shift.
• Linear increasing period: 0 < L ≤ 732.6 m. A 10 , A 30 , and A 50 all increase linearly with L, with rates about 57.46 km 2 /m, 10.68 km 2 /m, and 4.08 km 2 /m, respectively. It demonstrates that the subsidence basin was rapidly expanded during this period. • Temporary stagnation period: 732.6 m < L ≤ 971.5 m. A 10 , A 30 , and A 50 all change extremely slowly as L increases, and it can be assumed that the subsidence basin is not expanding outward temporarily during this period. • Re-expansion period: 971.5 m < L ≤ 1022.3 m. A 10 , A 30 , and A 50 increased with L again. Due to the limited amount of data, the increasing relationship between area and length cannot be clearly given during this period.
The above analysis shows that the outward expansion process of the surface subsidence basin in the study area did not increase completely with the mining of the working face, but there was a pause in the middle of the development. However, there is no relevant theory to explain this phenomenon reasonably yet. This anomaly is analyzed in Section 5.3 of this paper.

Discussion
In Section 4 of this paper, there were three abnormal phenomena found during the elucidation of the surface subsidence characteristics of the 2401 working face, which were: (1) a slight uplift appears on the east surface of the open cut eye; (2) obvious changes in the shape of the edge zone of the subsidence basin; (3) the expansion process of the surface subsidence basin undergoes a suspended expansion stage. This section discusses the possible explanations for the above anomalies.

Slight Uplift of the Ground Surface above the Open Cut
As mentioned in the Section 4.1, during the mining process of the 2401 working face, an abnormal slight uplift of the ground surface occurred on the east side of the open cut, forming a trend of rotation toward the center of the mining area together with the subsidence area on the west side. To further analyze this phenomenon, a vertical profile AA' was made along the strike centerline of the 2401 working face in this study, as shown in Figure 14. Referring to the results of the three-stage division derived from Section 4.2, the cumulative vertical displacement curves on 18 February 2020, 11 July 2020, and 31 January 2021, noted as VC1, VC2, and VC3, respectively, were plotted as shown in Figure 15. The above dates are at the end of stage one, stage two, and stage three, corresponding to the advancing lengths L of the 2401 workings of 117.6 m, 512.5 m, and 1022.3 m, respectively.
Overall, with the increase of L, the range affected by the mining gradually increased, the mining influence range gradually enlarged, and the maximum settlement also increased from 30 mm to 118 mm and finally reached 230 mm. Attention should be focused on the phenomenon that there was a slight uplift near the ground surface directly above the open cutting, as shown in the rectangle (a) and (b) in Figure 15. Figure 16a,b are partial enlarged views of two rectangles respectively.
Within the region shown in Figure 16a, the vertical displacement from the ground surface is dominated by a slight uplift of less than 20 mm. There exists an intersection point of vertical displacement curves VC1 and VC2, abbreviated as CP1. During the transition from VC1 to VC2, the relative vertical displacement of the surface on the left side of CP1 is in the form of subsidence, while on the right side, it is in the form of an uplift. It indicaes that the surface in region (a) underwent a rotation with CP1 as the center from 18 February 2020 to 11 July 2020. The above rotation phenomenon also existed between 11 July 2020 and 31 January 2021. Only the center of rotation was shifted from CP1 to CP2 along the advancing direction of the working face, as shown in Figure 16b. Therefore, it can be concluded that the uplift on the east side of the open cut is actually part of the overburden rotation. Studies [52,53] have shown that the thicker the rock stratum is, the higher the bearing capacity it has and it is more prone to rotation. The thickness of the Cretaceous Zhidan Group rock formation in the strata of this study area is generally over 300 m. Therefore, it is likely that the ultrathickness of the weak stratum caused the rotational motion characteristics during the overburden movement.
In contrast, there is no uplift on the ground in front of the working surface, but an overall downward bending occurs, as shown in the rectangle (c) in Figure 15. The characteristics of the surface vertical displacement in those two regions indicate that there is an apparent difference in the state of the stress distribution and overburden structure between the overburden in front of the working face and that directly above the open cut.

Influence of Groundwater and Aeolian Sand
Based on the SBAS-InSAR settlement monitoring results, it can be seen that the subsidence in the study area is very slight. According to the key stratum theory [54,55], it can be presumed that the stratum, which plays an overall controlling role on the surface settlement, is in a state of bending deformation but not yet destroyed. By analyzing the stratigraphic structure of the study area, it is clear that the giant thick weakly cemented sandstone layer with an average thickness of 341.33 m is the key strata that plays a role in controlling surface subsidence. However, due to the loose structure and high porosity, it is strongly permeable and has a certain hydraulic connection with the aeolian sand layer Q. Therefore, when the bending occurs, the groundwater in the two rock layers is transferred from all around to the bottom of the bending center along the original channels or newborn fine fissures. During the process, the groundwater in the nearby subsidence area migrates to the central zone of the basin, especially the area in front of the working face. On the other hand, after the migration of groundwater occurs, under the effect of gravity, the pores in the aeolian sand layer Q that used to be water-filled are immediately occupied by the surrounding sand particles, causing the compaction phenomenon of the aeolian sand layer, which in turn causes the surface subsidence, as shown in Figure 17. This is a typical soil consolidation phenomenon, similar to the surface settlement caused by direct groundwater pumping [56]. The latest research [57,58] shows that in the excavation process of geotechnical engineering such as deep buried pits, the pumping and drainage of groundwater will cause a large area of ground settlement, and if groundwater is not supplied effectively in time, permanent compression deformation will occur in the overlying strata under the consolidation of soil. Consequently, it can explain the widespread coverage of surface subsidence basins in the study area, and the sporadic subsidence along with the haphazard boundary.

Influence of Giant Thick Weakly Cemented Overburden
The analysis in Section 4.3 shows that the surface subsidence basins in the study area were suspended in the process of outward expansion. Combined with the key strata theory, it can be preliminarily considered that this phenomenon reflects the relative equilibrium stage of strata controlling surface subsidence in the process of overburden movement.
The theory of the three-zone distribution [59] in the extraction area divides the overburden into a collapse zone, a fracture zone, and a bending zone from bottom to top. With the extension of the mining area, the development height of the collapse zone and the fissure zone gradually increases, while the thickness of the bending zone gradually reduces. Considering the stratigraphic structural characteristics of the study area, it is clear that the giant thick Cretaceous weakly cemented sandstone layer is the key strata controlling surface subsidence and is also the major component of the bending zone. Bounded by the rock floor, the bending zone is further divided into two regions. Taking the floor of the rock stratum as the boundary, the bending zone is further divided into two regions, namely the upper bending zone (UBZ) and the lower bending zone (LBZ). Before mining activities, the LBZ zone basically has an upward supporting force on the UBZ. As the working face continuously advances, the supporting force gradually decreases until it reaches zero. Subsequently, the force of the LBZ on the UBZ is transformed into a downward tension. The magnitude of the tensile force increases as L enlarges at first. However, when L reaches a certain value, the deformation of the lowermost rock in the LBZ exceeds the limit, and the breakdown occurs. As a result, the tension effect of the LBZ on the UBZ is released to a certain extent. Combined with the fact of the basin paused expansion, it can be seen that although the lateral range of the UBZ is expanded and the gravity is relatively increased, the sum of the increased gravity and the decreased tension does not change significantly, but is exactly balanced with the support force of the surrounding rock body on the UBZ. After the broken section reaches a certain range, the increase of gravity exceeds the decrease of tension, so the equilibrium state of the UBZ is broken again, and the surface continues to deform synchronously with the giant thick weakly cemented rock layer.
Recently, some researchers have also found similar results during the dewatering of deep buried confined aquifers [60,61]. That is, the strata above the pumped confined aquifer experience tension stress and show some tensile deformations. Moreover, the essence is considered to be the relief of the pore water pressure. The above-mentioned force transformation process of the LBZ on the UBZ described in this paper is essentially the decompression of strata stress during coal mining. In fact, both the decompression of strata stress and the relief of the pore water pressure are essentially the unloading process. Therefore, this study summarizes the mechanism of the suspended expansion of the surface subsidence basin as "the unloading equilibrium of weakly cemented overburden", which is illustrated in Figure 18.

Conclusions
In this paper, the time series of surface subsidence during the mining process of the 2401 working face in the Yingpanhao coal mine, a deep mining area in the aeolian sand area of Western China, were obtained based on the SBAS-InSAR technology. According to the calculation results, the relationship between the geometric characteristics (shape and coverage area) of the surface subsidence basin and the advancing length of the working face in the study area was analyzed, and the following insights were obtained: (1) Under the condition of insufficient deep mining in the western aeolian sand area, the settlement area covered a wide area with a low value of settlement. When the advancing length of the 2401 working face reached 1022.3 m, the maximum surface subsidence was only about 250 mm, and the area of subsidence area reached 48.52 km 2 .
(2) The abnormal uplift phenomenon on the east side of the open cut was essentially part of the overburden rotation during the advancing process of the 2401 working face. By contrast, there merely existed an overall downward bending on the surface in front of the working face. The characteristics of the surface vertical displacement in those two regions indicated that there was an apparent difference in the state of stress distribution and overburden structure between the overburden in front of the working face and that directly above the open cut.
(3) According to the shape characteristics, the surface subsidence basins in the study area were divided into a central zone and an edge zone. The central zone was located directly above the 2401 working face, and its shape underwent a change from circular to oval, covering only one fifth of the total subsidence basin coverage. The edge zone was located in front of the working face, and its shape changed significantly, experiencing the circular-parallelogram-trapezoid change. The difference in shape reflected that the settlement mechanisms of the two zone were essentially different. It was concluded that the subsidence in the central zone was directly related to the mining activities, whereas a large-scale settlement in the edge zone was associated with the soil consolidation induced by the decline of groundwater.
(4) The quantitative relationship between the coverage area of the subsidence basin and the advancing length of the working face indicated that the outward expansion process of the surface subsidence basin in the study area was divided into a linear increasing period, a temporary stagnation period, and a re-expansion period. The temporary suspension of the expansion reflected the temporary relative equilibrium state of the giant thick and weakly cemented sandstone layer that controlled the surface subsidence. Based on the three-zone distribution theory of goaf, the bending zone was divided into the UBZ and the LBZ. When the advancing length of the 2401 working face reached 732.6 m, the rock layer at the bottom of the LBZ was broken, and the tension force of the UBZ by the LBZ was reduced. The UBZ reached a relative mechanical equilibrium state under the action of gravity, tension and the support force of the surrounding rock until the advancing length of the 2401 working face exceeded 971.5 m.
The research work in this paper provides an important reference for a reasonable solution to the contradiction between the exploitation of underground coal resources and the groundwater, as well as the ecological protection of the environment in the aeolian sand mining area located in Western China. Since the 2401 workings had been mined before the study was conducted, there was a lack of groundwater level observation data and overburden internal deformation damage observation data. It is recommended that the abovementioned missing data are collected to conduct an actual measurement test or an experimental verification from other mining activities to study the phenomenon and analytical results mentioned in this paper, so as to further clarify the mechanism of strata movement after deep coal resources mining in the aeolian sand area.