Study of Ground Movement in a Mining Area with Geological Faults Using FDM Analysis and a Stacking InSAR Method

Underground coal mining activities and ground movement are directly correlated, and coal mining-induced ground movement can cause damage to property and resources, thus its monitoring is essential for the safety and economics of a city. Fangezhuang coal mine is one of the largest coalfields in operation in Tangshan, China. The enormous amount of coal extraction has resulted in significant ground movement over the years. These phenomena have produced severe damages to the local infrastructure. This paper uses the finite difference method (FDM) 3D model and the stacking interferometric synthetic aperture radar (InSAR) method to monitor the ground movement in Fangezhuang coalfield during 2016. The FDM 3D model used calibrated Fangezhuang geological parameters and the satellite InSAR analysis involved the use of ascending C-band Sentinel-1A interferometric wide (IW) data for 2016. The results show that the most prominent subsidence signal occurs in mining panel 2553N and the area between panel 2553N and fault F0 with subsidence up to 57 cm. The subsidence observed for the FDM 3D model and stacking InSAR to monitor land deformation under the influence of fault are in close agreement and were verified using a two-sample t-test. It was observed that the maximum subsidence point shifted towards the fault location from the centre of the mining panel. The tectonic fault F0 was found to be reactivated by the coal mining and controls the spatial extent of the observed ground movement. The impact of dominant geological faults on local subsidence boundaries is investigated in details. It is concluded that ground movement in the study area was mainly induced by mining activities, with its spatial pattern being controlled by geological faults. These results highlight that the two methods are capable of measuring mining induced ground movement in fault dominated areas. The study will improve the understanding of subsidence control, and aid in developing preventive measures in Fangezhuang coalfield with fault reactivation.


INTRODUCTION
As China's primary source of energy, coal it is critical to the country's social and economic growth. Large-scale coal mining, however, has the potential to substantially harm the environment in mining areas, as well as causing a variety of geological issues and societal problems (Diao et al., 2019). As one of the main effects, a thorough understanding of mining-induced subsidence is essential for preventing or mitigating such issues.
With regards to measuring surface movement, traditional methods like levelling, Global Navigation Satellite System (GNSS), 3D laser scanning, and similar provide accurate information, but can be time-consuming, expensive, and have low spatial extent and are therefore unsuited to surveying large (basin level) areas. A better alternative is to use the satellite Interferometric Synthetic Aperture Radar (InSAR) technique, which provides weather independence, sunlight independence (active sensor), high (basin-level or greater) spatial coverage, and is less tedious and more economical. Therefore, InSAR offers a spatial resolution comparable to other traditional methods of land surveying. Many researchers have used InSAR to study land movement resulting from earthquakes (Fialko, Simons, and Agnew 2001), glacial movements , landslides (Sun et al., 2015;Kang et al., 2017;Zhang et al., 2020), volcanic bulging (Fournier et al., 2010;Albino et al., 2020), groundwater extraction (Bell et al., 2008;Motagh et al., 2017;Castellazzi et al., 2018) and coal mining Yang et al., 2017;Diao et al., 2018;Dong et al., 2021).
Extensive research has been carried out to study ground movement induced by underground mining activities (Park et al., 2012;Yang et al., 2019). In recent years, several research tools, such as numerical analysis (Deck and Anirudh 2010; Shabanimashcool and Li 2012;Shi M et al., 2021) and InSAR time-series analysis (Sowter et al., 2016;Gee et al., 2017;Grebby et al., 2019;Ghayournajarkar and Fukushima 2020) have been widely employed to analyse the mechanism of ground movement. Furthermore, several studies have also focused on ground movement associated with fault activation. Bell et al. (2005) reported mining subsidence from Great Britain, Germany and Colombia, and stated that mining area experiencing reactivation of faults should be surveyed properly before construction and a safe gap of at least 10 m should be maintained between the fault zone's edge and any structures. Moreover, Mohammady et al. (2019) employed Random Forest theory to analyse subsidence susceptibility and found that gap from the fault, elevation, slope angle and water table had the largest influence on ground deformation. Gumilar et al. (2015) and Pacheco-Martínez et al. (2013) discovered a direct link between fault position and high land deformation rates.
In terms of numerical analysis, the finite difference method (FDM) is commonly utilised within the FLAC 3D software because of its efficacy as a tool for solving rock mechanics and geomechanical issues. It can handle material heterogeneity, nonlinearity, complicated boundary conditions, ground condition pressures and gravity. The FDM model idealises the rock mass as a continuous medium that deforms according to a given constitutive law, satisfying compatibility and equilibrium criteria, and yields approximate partial differential equation solutions. Due to the discontinuities (e.g., geological fault) in the study area, a discontinuous numerical method was chose to study the ground movement based on the FDM analysis with FLAC 3D software. In FLAC 3D , the ground movement value of any point in the model can be monitored; it is typically more applicable for continuous and uniform subsidence. FDM analysis is widely used to study the ground movement in coal fields (Cheng et al., 2019;Parmar et al., 2019;Sikora and Wesołowski 2021;Yan et al., 2021).
Only a few attempts have been made to study 3D model simulation of ground movement in conjunction with InSARderived ground movement to understand the effect of a geological fault (Jeanne et al., 2019;Perry et al., 2020;Francioni et al., 2021;Shi Y et al., 2021). Although these previous studies suggest that faults do affect the ground movement, the nature and magnitude of effects have not been fully explored. The current approach uses the hybrid method, which is the combination of the numerical method and stacking InSAR analysis to study how the fault affects the ground movement induced by mining. To study the ground movement associated with subsurface mineral extraction at Fangezhuang coal mine, a fully 3-D elastoplastic FDM model was constructed. The results of the developed FDM model were verified by comparing the outputs of the model with the ground movement revealed by stacking InSAR analysis. Finally, the ground movement under the influence of fault in Fangezhuang coal mine is discussed.
This paper is divided in five sections. The section 1 provides an Introduction, and the section 2 describes the study area. Section 3 describes the methodology for studying the ground movement with the 3D FDM model and stacking InSAR technique. In section 4, we reveal the ground movement in Fangezhuang coal mine using FDM model and InSAR technique and perform a comparison between FDM model and InSAR result and the ground movement trend in Fangezhuang coal mine is revealed. Finally, conclusions and future scope of the study are discussed in section 5.
the drilling data, the fault extends about 3,000 m, and the drop ranges from 14 m to 37 m. The fault is a high-angle normal fault, inclined to SWW, with a dip angle of 70°-84°. F0 fault runs through the whole monocline structural area and has a significant impact on mining production. Meanwhile, there are associated faults of a specific scale with a significant drop on both sides of the F0 fault, and most of the associated faults here develop along the strike of the F0 fault. 1) The O stratum, 500-900 m thick, comprises fine-grained dolomite, stratified limestone.

Coal Seam and Rock Strata
2) The C stratum, 140-290 m thick, comprises siltstone, mudstone and fine sandstone. This formation also contains 1-3 layers of unstable thin coal seams.
3) The P stratum has a mean thickness of 550 m. From bottom to top, the sequence is split into the four forms showed below.
• The Damiaozhuang formation mainly contains siltstone, mudstone and medium sandstone. Four layers of the coal seam can be mined, namely No. 5 coal seam, No. 7 coal seam, No. 8 coal seam and No. 9 coal seam. • The Tangjiazhuang formation is mainly composed of coarse-to-medium sandstone, followed by fine sandstone. The lower strata are interbedded with 1-4 layers of thin coal lines. • The Guye formation mainly comprises of medium-coarse sandstone with a small amount of mudstone and siltstone. • The Wali formation is mainly composed of medium-coarse sandstone, fine sandstone and siltstone. At the bottom, a layer of aluminum mudstone with a thickness of about 4-5 m is developed. 4) The Q stratum is mainly composed of clay, sand and gravel layer. The Quaternary alluvium covers the whole Fangezhuang mine field. The thickness of alluvium varies from 54 m to 424 m and gradually thickens from north to south.
In this research, No. 5 coal seam, located at around 560 m deep in Permian strata, is the main workable seam. The mining panel 2553N which started in January 2016 and terminated in December 2016 along the northwest-southeast direction, is chosen to analyse ground movement. For the geological structure, F0 fault cut No. 5 coal seam along the southeast direction, resulting in a decrease of the depth of No. 5 coal seam at the east side of F0.

METHODOLOGY AND DATA USED
In this study, FLAC 3D was adopted for numerically predicting ground movement caused due by coal mining. FLAC 3D is a 3-D finite-difference computer program to solve geological problems (Kumar et al., 2016;Shi M et al., 2021). In FLAC 3D , the initial model is the geomechanical model with the engineering scale, which is able to simulate the real geological condition in the mining coalfield. The null model plugged in FLAC 3D is used to delete the elements to simulate the mining extraction activities. When the elements representing the coal are deleted to leave a void, the overlying elements will have a free displacement boundary and the stress field will redistribute, resulting in the overlying elements caving. The caving is spread upwards to the top surface, which represent the ground movement. The InSAR analysis used 25 Sentinel-1 ascending, Interferometric Wide (IW)

Establishment of the 3D FDM Model
The mining panel 2553N shown in Figure 1 is selected as the computational model domain on the basis of spatial distribution of the considered mining area. Based on the geological properties of the Fangezhuang coalfield, a 3D FDM simulation model of panel 2553N and F0 fault is set up, which is shown in Figure 3. The parameters of the FDM model are set as follows: the size of panel 2553N is 900 m × 180 m, No.5 coal thickness is 3 m, mining depth is 560 m, angle of the fault F0 is 70°, and step excavation distance is 75 m. The coal mineral extraction was simulated as long-wall mining with step by step excavation. The ground surface deformation map was obtained for the study area from January 2016 to December 2016 when the ground surface was stable (where the unbalanced force is less than 10 -5 of the maximum unbalanced force) (Du et al., 2019). The dimensions of the model are 1800 m long, 900 m wide, and 590 m high. The top surface does not show any geomorphic feature and is assumed horizontal in the model. The shortest horizontal distance between the mining boundary and the model boundary is 360 m to remove the boundary effect.

Generalisation of Strata and Faults
Mining panel 2553N lies in the Permian strata, and No. 5 coal seam is the main workable coal seam. As a result, two geological strata Quaternary and Permian, as well as No. 5 coal seam are simulated in a FDM model. The structure of each strata is made with reference to the geological report of Fangezhuang coal mine and borehole data. The fault could be designed with interface elements in FLAC 3D if the thickness of the fault is small (Cai et al., 2021). However, when the thickness of a fault is more than 10 m, it is appropriate to use a layer with a certain thickness to simulate a fault (Xu et al., 2013). In this paper, we choose a layer with two boundary surfaces to model the F0 fault around 10-20 m thick.

Computational Mesh
FLAC 3D supports numerous element shapes, including hexahedrons, tetrahedrons, pyramids, triangular prisms, and others (Abbasi et al., 2013). With a consideration of achieving a balance between processing time and simulation accuracy, a set of octahedral elements was adopted to mesh the simulation model. The meshing approach is as follows. It can be seen in Figure 3, from top to bottom, the model is divided into 13 layers which are loose layer, fine sandstone, siltstone, medium coarse sandstone, mudstone, medium coarse sandstone, siltstone, fine sandstone, siltstone, No. 5 coal, siltstone, No. 7 coal and siltstone. Each layer has a thickness of 3-160 m. The closer a layer is to the coal seam roof, the thinner the layer is set; the farther it is from the roof, the thicker it is set.

The Constitutive Model and Boundary Conditions
The Mohr-Coulomb model is used as the constitutive model and Mohr-Coulomb yield criterion are used in this study, along with displacement boundary conditions. No vertical displacement is allowed at the bottom of the simulation model, and no displacement is allowed in the direction perpendicular to the lateral boundaries. Trapezoidal distributed load is applied in the horizontal direction. Gravity force is applied to the model. A numerical model is obtained from these results, which reflects the mining and geological condition of the region after long-wall panel exploitation.

Monitor Points
The 3D  to obtain the ground movement during the coal mining process.

Mechanical Parameters of Rock Strata
Intact rock and discontinuities make up rock layers. The mechanical properties of complete rock samples acquired through laboratory testing differ substantially from the mechanical properties of rock strata (Hoek and Brown 1997). However, to some extent, the reliability of numerical simulation results depends on the choice of rock mechanical parameters. The methods for determining mechanical parameters mainly include empirical reduction method, engineering rock mass classification method and displacement back analysis method. This research adopts the back analysis with orthogonal test and numerical simulation to determine the rock mechanical parameters. The detailed procedure is as follows (Xu et al., 2013): 1) Based on the geological report of Fangezhuang coal mine, the initial rock parameters are shown in  Table 1. Therefore, each strata mechanical parameters could be given as The maximum ground movement W max is chosen as the test indicator.
2) The orthogonal test with four factors and five levels L 25 (5 4 ) was set up to test the rock parameters ( Table 2). All the testing schemes are listed in an orthogonal table (Table 3).
3) 25 numerical simulations were conducted using FLAC 3D and the results for each scheme were listed in Table 3. According to the result in Table 3, it was found that the maximum ground movement of the 4 th scheme is closest to the measured ground movement of 16.2 cm (geological report of Fangezhuang coal mine). Therefore, the parameters of the four factors are E 3.6GPa, μ 0.31, C 6.3MPa and ∅ 37°. The final parameters for each layer shown in Table 4.

Stacking InSAR Analysis
An InSAR stacking method was employed to map the ground deformation due to the decimetre scale rates of deformation expected to occur following the long-wall extraction of panel 2553N. The processing chain utilised is summarized in Figure 4. The Sentinel-1 SLC data were initially deburst and merged before coregistration to the slant-range coordinate system of the master image (January 2, 2016). Phase ramps attributed to orbital errors were subtracted using the precise orbit determination and topographic phase was simulated and removed using a DLR digital elevation model (DEM) from the TanDEM-X mission (Rizzoli et al., 2017). The data were multi-looked by a factor of 4 in range and 1 in azimuth and interferograms, with approximately a 10 m resolution, were generated between consecutive SAR acquisitions irrespective of the perpendicular baseline. The interferometric fringes (phase cycles that correspond to a displacement of half of the sensor wavelength) are related to the surface deformation. Generating interferograms only over the shortest epochs the minimises like likelihood of quantitatively underestimating the deformation due to the ambiguity effect when deformation gradients are high, such as over active mining sites. Further, an accurate estimation of the deformation depends upon high coherence (or phase correlation) between the two forming SAR images, hence, utilizing consecutive image pairs helps maintain coherence. In addition, the interferograms were filtered using a modified Goldstein filter to further improve coherence and the quality of phase. The interferograms were unwrapped from modulo-2π phase to relative deformation using a statistical-cost network-flow algorithm (Chen and Zebker 2001) with respect to a reference point located at Tangshan (39.6309°N, 118.1802°E). The unwrapped interferograms were subsequently stacked and an average rate of motion was derived from a least squares covariance analysis of the unwrapped phase. Once linear velocities had been generated, the relative height change for each image acquisition was calculated in accordance with that of Berardino et al. (2002). Finally, the line-of-sight time-series were projected into the vertical, by means of dividing by the cosine of the incidence angle (∼ 0.639 radians), to facilitate an appropriate comparison with the model.   (Zheng et al., 2018). From the perspective of ground movement, both the FDM model and InSAR measurements confirm that the area is unstable due to the coal mineral extractions. The ground movement basin is located within the mining panel 2553N and the area between the panel and fault F0. It can be seen from Figure 5A that the shape of the subsidence basin is not symmetrical over the mining area (mining panel 2553N), and the maximum subsidence point is located at point B, which is about 75 m from point A (centre of panel 2553N). The subsidence map shows a clear boundary nearby Fault F0. The boundary corresponds to the fault F0 and clear ground movement differences can be seen either side of fault F0. Figure 5B shows the similar subsidence basin revealed by FDM modelling. The maximum subsidence point is also located at point B but with a marginally smaller magnitude  Figure 5C shows the residual map of the FDM model and the InSAR result. The red colour represents where the ground movement of FDM model is less than the InSAR result. In other words, negative residual measurements represent underestimation of the FDM model with respect to the InSAR result. It can be seen that the average ground movement obtained by InSAR is smaller, especially in the area of severe ground movement between panel 2553N and fault F0. Nonetheless, the surface subsidence caused by long-wall mining is significant. InSAR only detects the medium-scale ground movement. As a result, the 3D FDM model proposed in this study is a valuable tool for ground movement analysis. Since the numerical model does not take the rock fissure and other geological properties of rock mass into account, the simulated results appear relatively smoother than the InSAR result. In addition, the InSAR measurements are relative to a reference point and are therefore not absolute. Thus, when comparing the InSAR results directly to the results from the FDM-model, it is worth noting that there may be some offset if the InSAR reference point was not entirely stable (i.e., it has a non-zero displacement) during the observed period. Furthermore, InSAR results may be affected by noise such as atmospheric, ionospheric (Liao et al., 2018) and unwrapping correlation error (Yunjun et al., 2019). Thus, any differences between the measurements could be due to several aspects relating to the InSAR processing, e.g., noise and the fact that the InSAR measurements are relative to a reference point that may be moving slightly.

Comparison Between FDM Model and Stacking InSAR Result
The ground movement pattern may result from underground mining, fault reactivation and other factors. Figure 6 shows the time series surface deformation of the FDM model and InSAR analysis at points A, B, C and D. It can be seen that the maximum subsidence values gradually increase from 14 th January to 20 th September. After 20 th September, both the FDM model and InSAR result show that the rate of increase of the maximum ground movement value decreases. Although the mining panel is advancing from September to December, the ground movement is increasing slowly at a relatively steady pace. According to mining subsidence theory, when the advancing distance reaches 1.2-1.4 times the average mining depth, the mining activity will reach full mining condition (Wang et al., 2020). This suggests that after 9 months of mineral extraction (panel 2553N advancing 675 m along the strike direction), the mineral extraction reached supercritical mining. From January to March, uplift is observed in the InSAR time series, which might be the noise effect in the stacking analysis. Compared with the ground movement value at point A, point B experiences less ground movement, confirming that the maximum ground movement is not located at the centre of panel 2553N. There is a clear difference in ground movement between C and D, which are located on different sides of fault F0.  Table 5 summarises the two-sample t-test results for points A, B, C and D (marked in Figure 5) to determine, quantitatively, whether the differences between the InSAR and FDM model series at these points are within the acceptable range (Agarwal et al., 2020). The mean difference between the InSAR and the FDM model is not the same at each of these points. As a result, it is necessary to investigate if the difference is significant. For the test, the Null Hypothesis (H0) is taken as the average difference is zero (µ 0) and alternate Hypothesis (HA) as the average difference is not zero (µ ≠ 0).
The t-value computed for the comparison of the InSAR and FDM modelling outputs at point A is 1.32. This means that these values are 1.32 standard deviations off of the mean. A t-value less than the critical t-value (2.069 for point A) is necessary to accept the null hypothesis. The t-value for point A is less than the crucial t-value, meaining that the null hypothesis is accepted, indicating that there is no significant difference between the InSAR result and the FDM model at point A. Another approach to confirm this is to look at the p-value (0.25) and alpha (0.05 for the 95 percent confidence limit), where a p-value greater than alpha means the null hypothesis is accepted. Acceptance of the null hypothesis in both cases demonstrates that the InSAR-derived subsidence exists, and that the FDM model is consistent at point A. From Table 5, similar results can be seen for points B, C, and D. As a result, we have sufficient evidence that InSAR-based subsidence is consistent and agrees well with FDM model values for all four points A, B, C, and D.
In order to examine the ground movement trends, the transect lines M-M′ along the coal seam strike direction, N-N′ along the coal seam dip direction and F-F′ vertically across the fault F0 were extracted (as shown in Figure 5). Figure 7 shows the land subsidence profiles of observation lines M-M′, N-N′ obtained from the FDM model and stacking InSAR analysis. It can be seen that the ground movement trends of the two lines are almost identical. The root-mean-square error along M-M′ in the dip direction is 0.113 m, with a maximum difference 0.141 m. The root-mean-square error along N-N′ in the strike direction is 0.105 m, with the largest difference of 0.210 m. The blue and orange lines show InSAR and FDM-model subsidence, respectively. The grey shading shows the standard deviation of the displacements observed at each location for both InSAR and FDM subsidence. The overlapping darker area for both the curves highlights that the deviation between the methods is small and within the standard deviation error limits. It can be seen that the From the subsidence profile of N-N′ in Figure 7B, it can be seen that the InSAR result shows the ground movement in the southern part of N-N' is slightly greater than that in the northern part. This is likely because the ground movement obtained by InSAR includes some residual subsidence caused by adjacent goaf (Fan et al., 2021). According to the geological report of Fangezhuang coal mine, mining along a panel to the south of panel 2553N stopped in 2015 and may cause the residual ground movement in 2016. Figure 8 shows the land subsidence profiles of observation lines F-F' obtained from the FDM model and stacking InSAR analysis. The grey shading shows the standard deviation of the displacements observed at each location for both InSAR and FDM model subsidence. The overlapping darker area for both the curves highlight that the deviation obtained from both the methods is small and within the standard deviation error limits.

Influence of Fault F0 on Ground Movement
For the profile along F-F′, it can be observed that at the distance of around 155 m there is a distinct change in slope in the subsidence profile. Interestingly, the fault F0 also lies at this location, which appears to be controlling this subsidence pattern. At this location, the fault is attributable for subsidence of 100-300 mm/ year on its right, and a subsidence rate of 300-500 mm/ year on its left. It can also be seen from Figure 5 that fault F0 has a noticeable effect on the spatial distribution of ground movement. This is most likely due to differential vertical compaction of various thicknesses of compressible soil deposited on both sides of the faults, which is greater in the hanging wall in normal faults (Burbey 2002). This observed structural control of land subsidence causes differential subsidence rates on either side of the fault line, which could potentially cause damage to villages, trains, and other infrastructure.
Fault F0 crosses the Fangezhuang coal mine in a north-south direction. The fault trace divides the subsidence rates and the arrows in Figure 8 indicate the relative displacement of the fault. To the west of F0, along F-F', the relative severe subsidence has a maximum value of 425 mm obtained for the FDM model result and 433 mm measured by InSAR. The east of F0 exhibits relatively stable subsidence, with a value up to 225 and 152 mm obtained with the FDM model and InSAR, respectively. In addition, Figure 5 shows the maximum subsidence point (point A) is not above the centre of panel 2553N, but it approaches F0. It can be seen that fault F0 has a blocking effect on overburden movement. The surface subsidence morphology under the influence of fault is not symmetric about    the centre of goaf. The maximum subsidence is shifted towards the fault. The ground movement patterns show a good agreement with the results derived from the recent study (Yu 2020). The precise different ground movement on both sides of F0 indicates that the fault is reactivated by panel 2553N coal mining. These results suggest that the subsidence in the study area is a dynamic process. Overall, the deviation obtained from the FDM model and InSAR method is small and within the standard deviation error limits. The InSAR-derived surface subsidence curve coincides with the FDM model simulation result.

CONCLUSION
This paper studies the ground movement in Fangezhuang coalfield in Tangshan city using a complete 3D FDM analysis in conjunction with stacking InSAR analysis. In order to obtain a reliable model, an orthogonal test was applied to back analyse the strata parameters, and the subsequent FDM model result was compared with the ground movement measured using InSAR. The differences in the maximum ground movement obtained from the FDM model analysis and stacking InSAR analysis were within 25%. Both results show that when the fault is activated by underground mining, the maximum subsidence value will not be located precisely above the centre of the mining panels. Instead, it is shifted towards the fault, and the geological fault clearly affects the spatial distribution of ground movement. The most dominant subsidence occurs in mining panel 2553N and the area between panel 2553N and fault F0, and has subsided by up to 57 cm. Overall, the ground movement patterns and magnitudes obtained by two different methods have a relatively good consistency, which attests the efficacy of the FDM numerical model in simulating the impact of the mining activity. The ground movement pattern in the Fangezhuang coal mine is spatially controlled by the geological fault to some degree, causing differential subsidence that could affect infrastructures. The method proposed in this paper can help to improve the understanding of subsidence control, develop preventive measures in Fangezhuang coalfield that consider fault reactivation, and forecast the impact of future underground mining activities in Fangezhuang coal mine. Further numerical analyses are required at a more realistic mine scale to study the ground movement trends in more detail.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

FUNDING
The project was supported by the University of Nottingham Faculty of Engineering Research Excellence PhD scholarship. The work was supported financially by the Fundamental Research Funds for the Central Universities (2019XKQYMS63).