Recent Ground Subsidence in the North China Plain, China, Revealed by Sentinel-1A Datasets

Groundwater resources have been exploited and utilized on a large scale in the North China Plain (NCP) since the 1970s. As a result of extensive groundwater depletion, the NCP has experienced significant land subsidence, which threatens geological stability and infrastructure health and exacerbates the risks of other geohazards. In this study, we employed multi-track Synthetic Aperture Radar (SAR) datasets acquired by the Sentinel-1A (S1A) satellite to detect spatial and temporal distributions of surface deformation in the NCP from 2016 to 2018 based on multi-temporal interferometric synthetic aperture radar (MT-InSAR). The results show that the overall ground displacement ranged from −165.4 mm/yr (subsidence) to 9.9 mm/yr (uplift) with a standard variance of 28.8 mm/yr. During the InSAR monitoring period, the temporal pattern of land subsidence was dominated by a decreasing tendency and the spatial pattern of land subsidence in the coastal plain exhibited an expansion trend. Validation results show that the S1A datasets agree well with levelling data, indicating the reliability of the InSAR results. With groundwater level data, we found that the distribution of subsidence in the NCP is spatially consistent with that of deep groundwater depression cones. A comparison with land use data shows that the agricultural usage of groundwater is the dominant mechanism responsible for land subsidence in the whole study area. Through an integrated analysis of land subsidence distribution characteristics, geological data, and previous research results, we found that other triggering factors, such as active faults, precipitation recharge, urbanization, and oil/gas extraction, have also impacted land subsidence in the NCP to different degrees.


Introduction
Groundwater constitutes the main water supply and meets more than 70% of the water needs in the North China Plain (NCP), a major agricultural base and industrial base in China [1]. Similar to several areas around the world [2][3][4][5][6][7][8][9][10][11], the NCP has experienced rapid land subsidence in association with extensive groundwater withdrawal. Moreover, rapid urbanization, increased dynamic and static loads, and the development of underground spaces have also impacted the occurrence of land subsidence to various degrees [12][13][14][15], and these factors also contribute to land subsidence in the NCP [16][17][18][19][20]. Land subsidence threatens the safety of urban infrastructure, increases the risk of urban waterlogging, exacerbates the impacts of sea-level rise on coastal areas, and even affects regional water security. Hence, subsidence monitoring and analysis are vital for detecting potential hazards and for maintaining sustainable regional development.
Subsidence in the NCP has previously been monitored using extensometers, levelling measurements, and global positioning system (GPS) measurements with high precision and spacetime restrictions [19,21,22]. Compared with these point-based geodetic methods, interferometric synthetic aperture radar (InSAR) has the ability to monitor large-scale ground subsidence throughout the day under all weather conditions with millimeter-scale precision. The development of multitemporal InSAR (MT-InSAR) addressed the limitations of traditional InSAR technology, such as temporal decorrelation, spatial decorrelation, and atmospheric delay. Accordingly, variants of the MT-InSAR technique, such as persistent scatterer InSAR (PS-InSAR) [23,24] and small baseline subset (SBAS) [25], have been successfully leveraged to analyze the time-series characteristics of land subsidence. Hence, considering the advantages of InSAR, this technique has been widely used to monitor subsidence areas throughout the NCP [26][27][28][29][30][31][32][33]. However, detailed InSAR subsidence studies always focus on specific sites, often Beijing, Tianjin, and Cangzhou. As a result, the spatial extent, magnitude, and temporal evolution of land subsidence elsewhere in the NCP have not been investigated.
The purpose of this paper is to fill this knowledge gap regarding land subsidence in the NCP by using the InSAR technique. The Sentinel-1A (S1A) satellite provides free C-band data with a swath width of 250 km; with these data, it is feasible to overcome the sparse data coverage and to perform large-scale subsidence monitoring. In this study, we apply the PS-InSAR technique to 3 ascending tracks acquired from 2016 to 2018 by the S1A satellite to obtain time-series displacements across the NCP and to detect the distribution of surface deformation therein. Then, to provide new insights into land subsidence in the NCP, we analyze the temporal and spatial evolution characteristics of and the main factors influencing land subsidence in this region.
The remainder of this paper is organized as follows. Section 2 introduces the history of land subsidence in the NCP and describes the study site. The Synthetic Aperture Radar (SAR) datasets and InSAR time-series method used in this study are presented in Section 3. Section 4 shows the InSAR results and the validation results. Finally, the discussion and conclusions of this study are provided in Sections 5 and 6, respectively.

Study Area
The NCP is an alluvial-proluvial plain formed from sediments carried by the Yellow River, Luanhe River, and Haihe River and covers the region between 34°46′N, 112°30′E and 40°25′N, 119°30′E, spanning a total area exceeding 140,000 km 2 . The NCP has a typical warm, semi-humid continental monsoon climate with a mean annual precipitation of 500-700 mm. The plain consists of the piedmont alluvial-proluvial plain (PP), central alluvial-lacustrine plain (CP), and littoral plain (LP) with quaternary deposit thicknesses ranging from 150 to 600 m [34]. The hydrogeological cross section A-A' from Shijiazhuang to Bohai Bay is illustrated in Figure 1b. The quaternary aquifers underlying the NCP are divided into four groups (Table 1): aquifer 1 is an unconfined aquifer at depths of 10-50 m, and aquifer 2 is a shallow confined aquifer at depths of 120 to 210 m, whereas aquifer 3 (depths of 250-310 m) and aquifer 4 (depths of 350-550 m) are confined and defined as deep groundwater [34][35][36]. Hereafter, the shallow aquifer refers collectively to aquifers 1 and 2. From the PP to the CP and LP, the aquifer systems generally vary from a single aquifer comprising sandy gravel to multiple aquifers composed of sand separated by silt or clay layers [37].
The first record of subsidence in the NCP was discovered in Tianjin by a levelling survey conducted in 1923. Before 1954, the NCP subsided in only parts of Beijing and Tianjin, and the annual subsidence velocity reached only a few millimeters. However, groundwater has been used as a major source of water since 1970; consequently, groundwater levels have fallen at alarming rates due to excessive exploitation [38]. Between 1975 and 1995, the average water table of both the shallow aquifer and the deep aquifer in some areas declined over 6 m and 23 m, respectively [39]. In response, land subsidence developed rapidly, extending from cities to rural areas [37]. The measured maximum accumulative displacement reached 800 mm from 1975 to 1995. As reported by the China Geological Survey (CGS), by the end of 2013, 42.3% and 10.8% of the total area of the NCP exhibited a subsidence velocity in excess of 20 mm/yr and 50 mm/yr, respectively (http://www.cigem.cgs.gov.cn). Among these areas, 90% of subsidence was distributed in Beijing, Tianjin, and Hebei and the cumulative subsidence in some parts exceeded 3 m [40]. In this study, we focus on land subsidence in the northwestern NCP, which includes the entire plain area around Beijing, Tianjin, and Hebei.

SAR Data and MT-InSAR Processing
The SAR datasets employed herein are composed of three tracks of images acquired by S1A, a C-band Earth observation satellite with a 12-day repeat cycle that launched on 3 April 2014. This study utilizes 404 S1A SAR ascending orbit images (tracks 40, 142, and 69) acquired with verticalvertical (VV) polarization in interferometric wide swath (IW) mode between 2016 and 2018. In IW Terrain Observation with Progressive Scan (TOPS) mode, S1A has a spatial resolution of 5 × 20 m and a swath width of 250 km (https://sentinel.esa.int/). Detailed parameters of the S1A data are shown in Table 2.  Figure 2) are as follows: master image selection, SAR data registration, digital elevation model (DEM) simulation, interferogram generation, PS selection with an amplitude stability index (ASI) threshold, atmospheric phase screen (APS) estimation and removal, and time-series deformation estimation. A DEM from the Shuttle Radar Topography Mission (SRTM) (http://dds.cr.usgs.gov/srtm/) with a 90-m resolution is used to remove the topographic phase. In addition, to ensure high coherence and stability, the PS candidates are chosen by employing an ASI better than 0.8. A previous study showed that the horizontal deformation within the NCP occurs with a relatively low velocity of 1.57-1.93 mm/yr with GPS observations [42,43]. Hence, we neglect the horizontal deformation and convert the line-of-sight (LOS) deformation into the vertical displacement by using Equation (1) [44]: where is the vertical displacement, is the LOS displacement, and θ is the incidence angle.
Then, the vertical displacements are validated with ground levelling data.
To generate a broad-coverage deformation velocity map, we need to mosaic the results from the three S1A tracks. However, multi-track SAR datasets are measured along different LOSs, and the locations of the detected pixels in each dataset also vary. Therefore, to accurately construct a broadcoverage subsidence velocity map, we guarantee three considerations. First, the LOS deformation from different tracks should be converted into vertical deformation. In this study, the area covered by track 142 is chosen as the reference area because it shares the same region as that with the data from the other two tracks. The incidence angles (Table S1) on track 142 and track 40 in the overlapping area are approximately 33.72° and 43.75°, respectively; similarly, the incidence angles on track 142 and track 69 in the overlapping area are approximately 43.77° and 33.67°, respectively. Each PS point from the three datasets is converted into the vertical deformation by Equation (1) according to the abovementioned incidence angles. Second, the locations of the PS points detected from each dataset differ in the overlapping area; accordingly, a nearest-neighbor search is used to match the PS points in the overlapping area from different tracks. Third, the least square method is used to calculate the deformation offset in the overlapping area between two neighboring tracks. Finally, the results from the three S1A tracks can be mosaicked by compensating for the mean offset.

Land Subsidence in the NCP during the Observation Period
We applied the PS-InSAR technique to the S1A datasets from 2016 to 2018 to investigate the ground deformation behavior in the NCP. To produce a broad-coverage deformation velocity map, we mosaicked the results from the three S1A tracks by compensating for the offset in the overlapping areas. Figure 3 shows the overall displacements derived from the S1A datasets by using the PS-InSAR technique in the NCP. The results reveal regions of significant subsidence mainly to the southeast and southwest of the PP, the southwest and east of the LP, and most of the CP. The overall displacement rate ranges from −165.4 mm/yr (where negative values denote subsidence) to 9.9 mm/yr (with positive values indicating uplift) with a standard variance of 28.8 mm/yr (during the acquisition of S1A data).
A total of 1,426,967 PS points was detected in the study area. Points with velocities higher than 20, 50, and 100 mm/yr accounted for 53.9%, 23.9%, and 3.6% of the total points, respectively. From 2016 to 2018, the percentage of points with land subsidence rates above 20 mm increased from 53.8% to 57.2%, while the percentage of points with subsidence rates over 100 mm presented a slight decreasing trend (Figure 3b). PS points with subsidence rates over 50 mm/yr accounted for 23.7%, 23.8%, and 24% of the total points during 2016, 2017, and 2018, respectively. During the observation period, the total area with a subsidence rate greater than 50 mm/yr reached 1.4 × 10 4 km 2 . We further calculated the area of subsidence exceeding 50 mm among the three sub-plains, and the results are presented in Figure 3c. In the PP, the area with subsidence greater than 50 mm presented a slight increase with an increase of 1.5% from 2016 to 2018. In the CP, the area with subsidence exceeding 50 mm presented a relatively stable tendency, constituting an increase of 0.3% between 2016 and 2018. Compared with 2016, the area with subsidence greater than 50 mm in the LP increased significantly in 2018 (an increase of 33.8%). From 2016 to 2018, the maximum subsidence rate was relatively stable in the CP and LP (Figure 3d). In the PP, the maximum subsidence rate first increased and then decreased. Figure 4 shows an example of the accumulative time-series displacements at selected points. For example, the accumulative subsidence at CY and WQ exceeded 30 cm in 3 years; the timeseries displacements at most selected points are similarly dominated by a permanent decreasing trend with minor nonlinear variations. In contrast, the time-series displacements at FN are characterized by significant seasonal variations: from September to January of the following year, the subsidence disappears or even rebounds (uplifts), after which the surface continues to subside.  We plot profiles of the displacement rate in the west-east (W-E) and north-south (N-S) directions to illustrate the spatial characteristics of subsidence areas in the NCP. The location of each profile is shown in Figure 5. The behavior of the spatial evolution of subsidence differs significantly at different subsidence funnels between 2016 and 2018. For example, the profiles across the Beijing, Tianjin, Shijiazhuang, Langfang-Tianjin, and Hengshui-Xingtai subsidence centers, labelled BJ-1/2, TJ, SJZ-1/2, LT, and HX ( Figure S1 in the Supplementary Material), respectively, are characterized by steady subsidence. The TS-1/2 profiles across the Tangshan subsidence center exhibit increasing subsidence rates ( Figure 6). In particular, the TS-2 profile indicates significant increasing subsidence rates from west to east and from north to south, with the maximum subsidence rate increasing from 81.3 mm/yr to 105.6 mm/yr and from 82.4 mm/yr to 106.9 mm/yr, respectively. On the other profile lines, we observe decreasing subsidence rates, for example, on the BD-2 profile across Baoding and the CZ-1 profile across Cangzhou (Figure 6). At distances between 11 and 15 km (from west to east) on profile BD-2, the subsidence velocity slows down from 2016 to 2018. The CZ-1 profile indicates that, at distances from 0 to 5 km, a small zone of subsidence slows down to the north and east. Combined with the statistical chart in Figure 3, these results indicate that the development of subsidence in the PP and CP was relatively stable from 2016 to 2018 and even slowed down in some areas, such as Baoding and Cangzhou. In contrast, the subsidence in the LP increased between 2016 and 2018, especially in the coastal areas of Tangshan.

Validation
In this section, we assess the precision and consistency of the subsidence rates derived from the individual SAR datasets. To evaluate the precision of the results, the InSAR-derived vertical displacement rates are compared with independent levelling measurements. In this study, 23 ground levelling survey benchmarks were acquired annually from September 2017 to September 2018 (shown in Figure 1 as red dots). Then, the average displacement of the pixels within a certain radius (100 m in our case) centered on each benchmark was obtained as the corresponding InSAR measurement. Figure 7 shows the comparison of the levelling data with the subsidence velocities obtained by two S1A datasets, showing generally good agreement. However, at some locations, differences of a few millimeters are detected. To investigate these differences, the root-mean-square error (RMSE) of the difference between the InSAR and levelling measurements is calculated for each levelling point, and the maximum and minimum errors between the two measurements are calculated. The overall RMSEs are 5.5 mm/yr and 7.2 mm/yr for tracks 142 and 69, respectively. The maximum and minimum errors for track 142 are 11.6 mm/yr and 0.2 mm/yr, respectively, and those for track 69 are 13.5 mm/yr and 1.0 mm/yr. Hence, the difference between the InSAR and levelling measurements is partially attributable to the differences in the location and acquisition time between the levelling and InSAR observations.

Comparison with Groundwater levels
As illustrated in Figure 9a, areas with high subsidence rates always correspond to low groundwater levels. The groundwater level contours drawn in Figure 9a are obtained from the China Geological Environmental Monitoring Institute (http://geocloud.cgs.gov.cn/). The distribution of land subsidence over the study region matches the spatial distribution of the sinking groundwater table, especially that of the deep groundwater level, as shown in the right panel of Figure 9a. As mentioned in Section 2, the shallow groundwater aquifer refers collectively to aquifers 1 and 2, whereas the deep groundwater aquifer refers collectively to aquifers 3 and 4. According to profiles A-A', B-B', and C-C' shown in Figure 9a, the shallow groundwater levels, deep groundwater levels, and subsidence velocities along the profiles are extracted to further analyze the spatial correlation between the hydraulic head level and land subsidence. Figure 9b shows the comparison between the hydraulic head levels from different aquifers and the InSAR-measured subsidence along selected profile lines. As shown in Figure 9b, the spatial distributions of the subsidence centers and the deep groundwater depression cones are fully aligned with each other with regard to profiles A-A', B-B', and C-C'. In other words, the subsidence centers are correlated with relatively low deep groundwater levels. On profile B-B' at distances between 0 and 120 km and between 240 and 285 km, the distribution of subsidence is similar to that of the shallow groundwater level. Similarly, on profile C-C', the subsidence distribution matches the shallow groundwater level distribution at distances between 0 and 45 km. In contrast, these two distributions are not consistent on the other sections. Groundwater exploitation in the NCP since the 1970s has lowered the hydraulic heads in the aquifers and formed groundwater depressions. A previous study pointed out that groundwater depressions are widespread in the PP, while deep groundwater depressions are distributed mainly in the middle and eastern NCP [39,45,46] (shown in the left panel of Figure 1a). According to Terzaghi's theory, groundwater overexploitation has caused a shift in the pore fluid pressure, causing an equivalent change in the effective stress within the aquifer system that compresses the aquifer system skeleton under the new load [47,48]. Thus, subsidence has occurred. Compared with the shallow aquifer, the groundwater within the deep aquifer has a poor ability to recharge within the study region [49]. Therefore, the contradiction attributed to the unbalanced replenishment of the deep aquifer is more prominent than that of the shallow aquifer, and hence, the spatial response of the former to subsidence is more obvious. Figure 10a depicts a contour map of the 40-mm interval superimposed onto the InSAR-measured subsidence and highlights four subsidence areas with maximum velocities over 100 mm/yr (marked A to D in the left panel). To examine the land use types corresponding to different magnitudes of subsidence, the land use classification is obtained from the Finer Resolution Observation and Monitoring-Global Land Cover (FROM-GLC10) map provided by Tsinghua University (http://data.ess.tsinghua.edu.cn/) with a resolution of 10 m. In this study, we integrated the classification results, and the final outcome is illustrated in the right panel of Figure 10a. The entire region is categorized into five main classes: cropland, vegetation, water, impervious land, and bare land. Figure 10b shows the statistics of these land use types within the ranges of different subsidence contours. In the areas with subsidence greater than 20, 60, and 100 mm, cropland is the most distributed land use type, followed by impervious land. As illustrated in Figure 10c, from 2003 to 2018, the percentage of agricultural water consumption ranged from 68% to 54% of the total water withdrawn. Agricultural water evidently consumes a large portion of the water in the study area. In the NCP, irrigation water is mostly pumped from groundwater, and high-intensity exploitation for irrigation causes the NCP to experience severe groundwater depletion [50,51]. Therefore, throughout the study area, agricultural planting has a major impact on land subsidence.  Figure 11 shows magnified views of the land use types in the four selected zones marked with black rectangles in Figure 10a (left panel). The detailed statistical information of these subsidence zones with velocities exceeding 100 mm/yr is provided in Table 3. In zone A, land subsidence occurs mainly in urban areas, and the main land use type in the areas with subsidence rates over 100 mm/yr is impervious land. In zone B, the main land use types in the area with subsidence rates over 100 mm/yr are cropland and impervious land, and their ratio is close to 1:1. As illustrated in Figure 11, there are considerable differences between the land use types in the two areas in zone B (marked B1 and B2) with subsidence rates exceeding 100 mm/yr. In zone B1, the major land use type is impervious land, but that in zone B2 is cropland. In zones C and D, the main land use types are cropland in the areas with subsidence velocities over 100 mm/yr; the ratios between cropland and impervious land are 1:3.7 and 1:9.4 for zones C and D, respectively.

Other Factors
As mentioned in Section 2, the NCP is divided into three main hydrogeological units. From the PP to the CP and the LP, the aquifer systems change from a single aquifer consisting of sandy gravel to a multilayer structure comprising sand separated by silt or clay layers. The NCP began to exploit and apply groundwater on a large scale since the 1970s; consequently, groundwater overexploitation problems exist in many areas at present. The areas in which shallow groundwater is overexploited are distributed mainly in the PP, and those in which deep groundwater is overexploited are distributed mainly in the CP and LP (Figure 12). Deep groundwater is characterized by slower renewal and less recharge than shallow groundwater [52,53]. Furthermore, precipitation is the main source of groundwater replenishment in the NCP: the spatial distribution gradually weakens from the PP to the LP and then increases [54,55]. The alluvial-proluvial fan area easily accepts recharge in the form of lateral inflow in the PP, and the conditions supplying atmospheric precipitation are relatively good. In contrast, groundwater presents weak rechargeability in the CP and LP due to the relatively poor groundwater recharge conditions therein. Therefore, under the same groundwater exploitation intensity, the CP and LP are more prone to land subsidence. Figure 12. Precipitation and groundwater exploitation information in the study area (the groundwater exploitation information is derived from [41,56]).
The development of subsidence is controlled by faults, and the main deformation direction is basically parallel to the fracture direction in the NCP [27,[57][58][59][60]. Figure 13 demonstrates that faults coincide with the boundaries of subsidence areas, while greater subsidence rates are directly related to the groundwater level. The information of faults in Figure 13 are shown in Table 4. As seen from Figure 13a, there are several active fault zones, and sharp gradients in subsidence are observed across the faults in the study area. Figure 13b shows an example of the average subsidence in four fault zones along profile AA': the Tongxian-Nanyuan fault zone (TNFN), eastern Taihang fault zone (TFZ), Wuqing-Wenan fault zone (WWFN) and Cangdong fault zone (CFZ). The profile of the average displacement rate along AA' shows marked differential differences on both sides of these faults, especially WWFN and CFN, with displacement gradients up to 160 mm/yr. These highvelocity gradients correlate with faults, indicating that these faults coincide with the boundaries of the subsidence areas. We observe a significant subsidence center correlated with a relatively low groundwater level, particularly between the WWFN and CFN. In addition, land subsidence is also connected with compressible quaternary layers, urbanization, and geothermal and oil exploitation [61][62][63].

Conclusions
In this paper, we first applied SARProZ MT-InSAR technology to multi-track Sentinel-1A images collected from 2016 to 2018 to detect the spatial and temporal distributions of land subsidence across the NCP. The vertical displacements were estimated to have a maximum value of approximately −165.4 mm/yr. The InSAR results are in good agreement with ground levelling measurements and were cross-validated, indicating the reliability of our results. During the observation period, the rates and spatial distributions of subsidence in the PP and CP have remained unchanged for at least three years. In contrast, in the LP, especially in the coastal areas of Tangshan and Qinhuangdao, the subsidence rates have increased and the spatial distribution has expanded significantly. An analysis of the observed displacements coupled with information obtained from geological mapping and groundwater level data illustrates that the distribution and development pattern of land subsidence in the NCP is correlated with deep groundwater exploitation and is controlled by geological structures. Moreover, comparing InSAR-derived subsidence maps with maps of land use types demonstrates that groundwater pumping for agricultural activities is the major mechanism causing land subsidence in the northwestern NCP.