Spatial and Temporal Evolution of Ground Subsidence in the Beijing Plain Area Using Long Time Series Interferometry

Due to the overexploitation of water resources, ground subsidence is becoming increasingly problematic in Beijing, China's political, economic, and cultural capital. This article aims to investigate the relationship between ground subsidence and changes in groundwater depth, and water supply from a long-term point of view. Multisource synthetic aperture radar (SAR) data using the interferometric SAR (InSAR) technique were adopted in this research, combined with a set of leveling and ground subsidence data in the Beijing Plain area from 2003 to 2020. The InSAR results demonstrate that ground subsidence in the plain area increased steadily from 2003 to 2015, expanding from sporadic to continuous laminar dispersion and producing five major subsidence centers. The South-to-North Water Diversion Project (SNWDP) that was completed in 2008 and 2015 considerably reduced the demand for groundwater supply in the Beijing Plain area. Since then, the groundwater level depth has continued to increase. However, since 2016, the ground subsidence rate has dramatically slowed down. The obtained results showed that, thanks to the SNWDP, which resulted in a decline in groundwater exploitation and an increase in renewable water recycling, the ground subsidence in Beijing's plain area has been effectively managed.


I. INTRODUCTION
G ROUND subsidence is a major geological phenomenon that threatens the survival and long-term development of millions of people around the world, particularly in coastal and highly developed urban areas. It includes both gradual downward curvature and sudden segmental subsidence [1]. Ground subsidence can harm subsurface and aboveground infrastructure (such as buildings, railways, and pipelines), modify the surface configuration [2] and induce soil fractures [3]. The detrimental effects of ground subsidence have already been recognized by government agencies, industry, and academics. Exact knowledge about the spatial and temporal evolution of ground subsidence, as well as its developmental processes, is critical for mitigating the impact of disasters.
Ground subsidence, mainly as a result of groundwater exploitation, which decreases the aquifer's water level, has become a worldwide geoenvironmental issue. Many regions have faced serious ground subsidence problems due to groundwater overdraft, such as Fars Province, Iran [4], Bangkok, Thailand [5], [6], [7], the southern United States [8], Mexico City [9], and Indonesia [10]. Ground subsidence has also occurred in various areas of China as a result of increasing urbanization, including the North China Plain and southern China [11], [12], [13], [14]. In these areas, groundwater is the primary supply of water for industrial, agricultural, and domestic usage, and it has been severely mined in recent decades.
Beijing, one of China's largest cities, has experienced significant population growth, industrialization, and urbanization during the past two decades. Groundwater has provided twothirds of Beijing's water supply in recent decades. According to research released by the China Geological Survey and the Beijing Water Authority, due to the long-term overextraction of groundwater, the aquifer system's groundwater level declined constantly from 2000 to 2015. The largest decline in water level in the 20 years was 10.39 m, and the area of major groundwater decline (where the groundwater depth is deeper than 10 m) grew by 2972.   LAST FIVE YEARS to investigate the relationship between the ground subsidence and excessive groundwater extraction.
Ground subsidence research was carried out in Beijing two decades ago [15], [16], [17], [18], [19]. Previous research, as shown in Table I, concentrated on trend monitoring, subsidence correlation analysis, groundwater spatial and temporal dynamics, and geological formation correlation analysis. Only a few ground subsidence locations were found in the eastern and northern plains of Beijing between 1992 and 2000, according to early InSAR and ground monitoring equipment monitoring. Ground subsidence in plains continued to expand from 2003 to 2010. When compared to the preceding era, the region where ground subsidence exceeded 50 mm per year grew by 63% [13]. Early investigations found a positive association between surface subsidence and groundwater exploitation in 1971-1987 and 2003-2015 [20], with the pattern of spatial development controlled by geological fault distribution [21]. Some researchers have determined that rainfall in Beijing will rise over the next decade, making rapid groundwater recovery more plausible [22]. Previous studies, on the other hand, have not thoroughly examined the relationship between water resources and subsidence, that is, how subsidence expands, how changes in urban water supply and water use structure affect the temporal and spatial evolution of subsidence, and how changes in groundwater level affect the temporal and spatial evolution of subsidence. Furthermore, even if there is a link between subsidence and groundwater level variation, there is little comprehensive study on the impact of water supply structure evolution on land subsidence.
Land subsidence is a long-term and irregular surface movement phenomenon. However, due to the lack of long-term or continuous monitoring analysis of current subsidence monitoring data in Beijing, it is difficult to identify the evolution law of the ground subsidence. The purpose of this research is to monitor land subsidence in the Beijing Plain area macroscopically through long-term SAR data to achieve continuous land subsidence update. The small baseline time series InSAR measurement technique [23] is used with multiple SAR images and leveling datasets spanning the plain area of Beijing from 2003 to 2020. First, by combining the time series results from three different SAR datasets, ERS, COSMO-SkyMed, and Sentinel-1, the spatial and temporal changes in ground subsidence in Beijing were revealed. Second, the InSAR results were cross-checked with leveling data. Furthermore, the annual fluctuation of ground subsidence in the plain region was studied. The interrelationships between ground subsidence and groundwater level depth, water supply structure were finally discussed. The remainder of this article is organized as follows: Section II describes the study area's geographic location as well as the data used for processing and validation; Section III presents the results and precision validation of the three types of ground subsidence data; and Section IV discusses the spatial and temporal variation in ground subsidence, its interrelationships with water supply structure, groundwater depth, as well as the factors that affect ground subsidence. Finally, Section V concludes this article.

A. Study Area
Beijing is located in the northern part of the North China Plain, with continuous mountains to the west, north, and northeast, as well as vast plains to the southeast. The Beijing Plain is a typical piedmont alluvial plain produced by sediments transported by the Yongding, Chaobai, and Wenyu Rivers. The city's total land  area is 16 807.8 km 2 , with a plain area of 6390.3 km 2 [24], [25]. The plains are the primary location for human production and dwelling. The footprint of human activity has been extending and accumulating in plain areas with increasing populations.
The demand for water resources rises dramatically as the population grows. The amount of water consumed in plain areas has expanded significantly in recent years, well beyond the carrying capacity of water resources. To alleviate the problem of water scarcity, the use of industrial and agricultural water has been reduced over the last 20 years, whereas investment in water conservation facilities, recycling water, and the south-to-north water diversion project (SNWDP) has increased, which has helped to alleviate the problem of water scarcity to some extent.
B. Data Sources 1) InSAR Datasets: Since SAR data from one sensor are not capable of providing a long time series from 2000 to 2020, the InSAR dataset comprises three different types of data: C-band ENVISAT ASAR data acquired from June 2003 to 2010, X-band COSMO-SkyMed data acquired from 2013 to 2015, and C-band Sentinel-1 data acquired from 2015 to 2020. Unfortunately, there were still some gaps between the consecutive SAR datasets. The coverage of SAR images and the temporal distribution of SAR images are shown in Figs. 1 and 2. The parameters of the SAR datasets are summarized in Table II. 2) Ancillary Datasets: The Beijing Water Resources Bulletin (BWRB) (http://swj.beijing.gov.cn/zwgk/szygb/) provided the groundwater level, water supply, and water use statistics. It primarily contains information on water consumption by various industries, water supply statistics from various water sources (surface water, groundwater, SNWDP, renewable water, etc.), and yearly groundwater depth in Beijing from 2000 to 2020. A 30-m digital elevation model obtained from the Shuttle Radar Topography Mission [28] was used for topography phase removal and geocoding.
In this study, level measurements were obtained as validation data to check the accuracy of the InSAR results, and the distribution of points is depicted in Fig. 1. The measurements covered the years 2005 to 2020, and each location reported annual subsidence measurements.

A. Small Baseline Subset (SBAS) Methods
The SBAS method by combining multiple interferograms appropriately in order to limit the temporal decorrelation phenomena [3], [23]. First, data are processed by correspondent multilook, which is designed based on the resolution ratio proportionate relationship between azimuth and range direction of a single look complex. ENVISAT ASAR data adopt 1 × 5 multilook, Sentinel-1 data adopt 10 × 2 multilook, and COSMO-SkyMed data adopt 3 × 3 multilook. The interference network is set using small temporal and spatial baselines thresholds, assuring the coherence of the interferogram. For the three different types of SAR data, different temporal and spatial baseline combinations are used. The ENVISAT-ASAR data spatial and temporal baselines are set to 1000 m and 200 days, respectively. For Sentinel-1 data, the two parameters are 150 m and 36 days; for COSMO-SkyMed data, the spatial baseline is 1000 m, without limited temporal baseline. Since the data mainly cover urban areas. Flat-earth phase reduction is conducted for each interferometric combination based on the SRTM DEM data in 3" radians throughout the differencing process using the interference network. The minimum cost flow method is used for phase unwrapping to ensure the accuracy of the obtained unwrapping phase. A coherence higher than 0.4 is used as the masking threshold for unwrapping. According to the data processing process, unwrapped data are used as the input data. First, the elevation correction phase is calculated, and the orbit error is estimated. Here, the precise orbital parameter information of ENVISAT-ASAR and Sentinel-1 data were used to remove the orbital errors. For the COSMO-SkyMed data, we use the quadratic surface fitting method to extract the orbit phase. Second, the trend errors, such as residual orbital and atmospheric phase, are filtered out utilizing high and low pass filtering. Due to the large spatial coverage and magnitude of the deformation zone, the atmospheric phase has certain similarities with the deformation phase. Therefore, the deformation zone is masked according to the interferogram and phase average rate, and in the filtering process, the atmospheric phase is estimated in the stable region to reduce the phase loss of the deformation zone. The Laplace interpolation method is adopted for the interpolation, through which the phase loss of the deformation area can be effectively reduced and the accuracy of the data processing results can be ensured. Finally, the time series deformation measurements for each image element position are obtained using singular value decomposition. Fig. 3 depicts the overall workflow of the strategy employed in this article.
Three types of radar data are registered before data processing to ensure that each scene image has the same geometric structure. The Sentinel-1 data were imaged through terrain observation by progressive scans mode. The enhanced spectral diversity method was used for the registration [29]. Based on the aforementioned criteria, 280 interferograms were constructed, of which 96 were from C-band ENVISAT, 47 were from X-band COSMO-SkyMed, and 137 were from C-band Sentinel-1. These three datasets were then processed separately by the SBAS method.
The InSAR technique captures the deformation along the line-of-sight (LOS) direction due to side-looking imaging characteristics. In the case of large surface ground subsidence, global navigation satellite system (GNSS) monitoring data from 2006 to 2010 revealed that the vertical displacement of the ground subsidence area was 1 mm per year in the east-west direction and 0.7 mm per year in the north-south direction [21], with the vertical component controlling the majority of the ground displacement. Therefore, the displacement time series in the LOS directions from the three sources of data were converted to the vertical direction to generate a continuous vertical displacement time series. Considering the data discontinuity, during the first data gap between ENVISAT and COSMO-SkyMed in June 2010-January 2013, the subsidence rate is assumed to be equal to the subsidence rate in June 2003-June 2010. Similarly, the subsidence rate in the second data gap is assumed to be the same as the annual subsidence rate in 2015 from the COSMO-SkyMed dataset.

B. Validation of InSAR Results
Many other measurement approaches or phenomena, such as GNSS, leveling, ground fracture, and field surveys, can be used to verify InSAR deformation results. For cross-validation, four leveling point L2, L3, L11, and L18 with large deformation values and their surrounding InSAR deformation points are selected for analysis. As mentioned before, there is a two-year interval between the three SAR datasets, so the obtained cumulative subsidence time series are discontinuous. However, the three discrete time series of cumulative subsidence can be linked together by continuous annual leveling of the plain area, generating continuous land subsidence time series. Fig. 4 shows the subsidence time series of the four leveling points L2, L3, L11, and L18 and the adjacent InSAR points. It is found that the subsidence time series produced by InSAR is in good agreement with the subsidence evolution recorded by leveling measurement, which is consistent with the accuracy of the land subsidence mapped by InSAR. , and Sentinel-1 (August 2015-December 2020) datasets were processed using the data processing methods described in Section III-A. Depending on the temporal and spatial baseline properties of the three data sources, different combinations of spatial baselines are employed to construct interferograms. The COSMO-SkyMed was unable to cover the plain area entirely, but it can cover the subsidence area.  After 2016, the land subsidence rate gradually slowed down, and the signal disappeared or weakened in some severely subsidence regions. The deformation rate decreased dramatically from 2016 to 2020, and the deformation signal vanished in some regions, such as Jinzhan and Chaoyang, where the capital airport is located, as shown in Fig. 4, by comparing the findings of the leveling points and InSAR feature points. The deformation area in the city's plain area gradually shrank.

B. Analysis of the Causes of Land Subsidence
Ground subsidence in the Beijing Plain was caused by the long-term overexploitation of groundwater, according to numerous studies [26]. Water scarcity has long been a major stumbling block in Beijing's socioeconomic development [30], [31]. Water use for agricultural irrigation, industry, and domestic usage continues to rise as the regional economy develops, disrupting the initial input-output balance. Many locations are utilizing deep pressured water to meet rising water demand, and groundwater has been overexploited for years, resulting in a dramatic reduction in water volume due to a lack of water level recovery. Ground subsidence was induced by a drop in regional differential pressure, a decrease in pore pressure, and an increase in effective stress due to long-term overexploitation of groundwater. As a result, the key elements examined in this study are groundwater change and water supply. Fig. 6 shows how land subsidence in Beijing's plain area progresses from a sporadic distribution to a contiguous state, indicating that long-term persistent driving causes are to blame for the development of land subsidence in the plain area. Although we cannot obtain reliable groundwater exploitation information at the ground level, the interaction between them reveals a substantial association between the typical deformation area and average groundwater variation in the plain region. Fig. 7 illustrates that the relationship between the depth of the groundwater level and the change in groundwater storage is related to the time series deformation rates for four typical deformation zones. Fig. 7 shows CYJC, CJ, SHT, and regions, from 2004 to 2016, the groundwater level decline was consistent with land subsidence. With the decline of the water level, the subsidence rate deepened year by year. From 2016 to 2020, the groundwater level rose and the subsidence rate slowed down year by year. In the HG region, area with the lowest water level in 2014 was also that with the smallest subsidence value. After 2014, as the groundwater level rose, the subsidence rate slowed down. On the whole, due to the adjustment of external water supply and water use structure, the groundwater level first decreased and then increased in the Between 2016 and 2020, groundwater levels rose year after year, land subsidence slowed, and subsidence phenomenon in some regions, such as the Chaoyang Airport area, disappeared. The relationship between the change in groundwater storage and the subsidence rate, on the other hand, demonstrates this notion. When the groundwater level drops and the groundwater storage drops, land subsidence accelerates; when the groundwater level rises and the groundwater storage increases, land subsidence slows down. As a result, groundwater alteration is one of the reasons for land subsidence in plain areas. However, we discovered in the correlation study that there was a slight difference between groundwater changes and land subsidence in some regions, which could be due to aquifer characteristics or geological sediment spatial changes.

C. Water Supply Structure Evolution and Subsidence Development
Groundwater exploitation caused by shortage of water resources is an important factor that induces land subsidence in Beijing. In order to solve the contradiction between supply and demand of water resources in Beijing, a series of adjustment measures, including water supply and water use, have been carried out from 2000 to 2020. In terms of water supply, the SNWDP is adopted to expand the supply of water resources, increase the utilization efficiency of recyclable water resources, and reduce the exploitation of groundwater. Before 2000-2015, groundwater supply accounted for the largest proportion of Beijing's total water supply, with an annual average of about 63%, the area of severe groundwater decline is 3180 km 2 . The spatiotemporal development of land subsidence in the plain area and groundwater exploitation shows consistency in this stage. When groundwater exploitation enters a period of large-scale overexploitation, land subsidence also starts a period of rapid development, the subsidence rate of the subsidence center increases, and the speed of spatial expansion accelerates. From 2000 to 2015, the spatial expansion trend of the area affected by land subsidence in Beijing is as follows: the eastern subsidence area expands to the   south, north, and east; The subsidence areas centered on Shahe and Baxianzhuang in Changping, Pinggezhuang, and Tianzhu in Shunyi, Xixiaoying in Haidian, and Laiguangying in Chaoyang are gradually formed, and connected with the eastern subsidence area. In 2015, groundwater water supply accounted for 47% of the total water supply in the whole year. This is the first time that the water supply of surface water, inflow from the SNWDP, and renewable water resources exceeded the annual water supply of groundwater, accounting for 53%. The water supply pattern in Beijing has entered a new stage. The adjustment of the water supply structure mainly focuses on the inflow of the SNWDP and the utilization of renewable water resources. Among them, the SNWDP is the main source of water supply for the six urban areas of Beijing (Dongcheng, Xicheng, Chaoyang, Haidian, Fengtai, and Shijingshan), which is a significant strategy to effectively control the rapid development of land subsidence in the main urban area, gradually realize the prohibition of groundwater overexploitation in the urban area, and at the same time select areas and opportunities with conditions to replenish local groundwater with southern water. After the adjustment of the water supply structure in 2015, the development trend of land subsidence has been significantly slowed down. Figs. 7 and 8 show the curve of annual subsidence rate and the acceleration curve of subsidence rate in the four subsidence funnels of CYJC, CJ, HG, and SHT, respectively. In Fig. 7, from 2015 to 2016, it can be seen that the deformation of each subsidence funnel has slowed down a lot. In Fig. 8, if the acceleration of the subsidence rate is greater than zero, it indicates that the upward supporting force (uplift pressure) of the underground soil layer on the ground is greater than the overall downward pressure on the ground. The sum of the forces, at this time in the subsidence area, are opposite to the direction of the subsidence movement. While the acceleration of the subsidence rate is less than zero, it means that the upward support force of the underground soil layer on the ground is less than the sum of the downward force of the whole ground. At this time, the ground force in the subsidence area is in the same direction as the subsidence movement. The subsidence acceleration of the four subsidence funnels gradually increased after 2015, which is related to the increase of groundwater level and uplift pressure year by year. The accelerated analysis also proves the relationship between land subsidence and water level changes. The reasons for the difference in the time of deformation mitigation between the settlement funnels are analyzed: First, the time when the water from the SNWDP provided for each region is different. Although the annual water supply of groundwater decreases, there is still a certain amount of groundwater overexploitation in local areas. Second, in the long-term groundwater exploitation process, the self-provided groundwater wells were not all closed in 2015, and there was still some illegal groundwater exploitation in some areas. Third, the underground soil layer will undergo two processes of compression and rebound during the period of groundwater level rise. The geological structure of different subsidence funnels leads to inconsistent soil compression and rebound states, resulting in the time difference in the subsidence rebound and the water level rise. Fig. 9 shows the changes in the main water use structure in Beijing from 2000 to 2020. It is worth noting that agricultural and industrial water use continued to decline and reached the lowest value of 300 million cubic meters in 2020, effectively easing the pressure on Beijing's water supply to a certain extent. In the past, most of the agricultural irrigation in the areas not covered by the water conservancy pipe network came from self-provided wells. Disorderly mining of the self-provided wells causes the groundwater level to drop and forms a settlement funnel. Industrial restructuring and technological development have reduced industrial water consumption. While industrial and agricultural water consumption has been decreasing year by year, environmental water consumption and domestic water consumption have continued to increase. As can be seen from Fig. 10, since 2000, Beijing's environmental water consumption has increased, mainly due to the increase in the area of Beijing's parks. Comparing Fig. 10 with Fig. 9, a stable "supply and demand balance" between renewable environmental water and reclaimed water was maintained in 2020. Domestic water consumption has increased due to continued population growth and changes in per capita water consumption. From 2007 to 2020, the biggest pressure on Beijing's water supply comes from the continuous growth of domestic water use, which can be explained by the improvement of living standards, the gradual popularization of water utensils and household bathing equipment, and the rapid growth of industries. Therefore, in terms of water use structure, the increase in domestic water use has a strong correlation with the development of the subsidence. Second, the subsidence centers shown in Fig. 6 are all urban expansion and rapid population growth areas, and the proportion of domestic water in these areas is much larger than that in other areas on the plains.
Through the analysis of Beijing's water supply and water use structure, over the next few years, the general water supply pattern of Beijing is likely to follow a steady development of existing conditions. The groundwater supply will gradually decrease, the supply of renewable water sources will continue to increase and exceed the groundwater supply in the short term, and the depth of the groundwater table is expected to rebound further.

V. CONCLUSION
Multisource radar (ENVISAT ASAR, COSMO-SkyMed, Sentinel-1) data and the small baseline method were utilized to investigate the ground subsidence over the past 17 years in the Beijing Plain in this article. The analysis of the relationship between InSAR results, level observation data, water supply data, and groundwater depth data could draw the following conclusions. First, the cross-validation between InSAR and leveling observations proved that the accuracy of InSAR subsidence monitoring in this study could reach 10 mm per year. Meanwhile, the large area of ground subsidence in the plain area was mainly composed of vertical components, whereas the horizontal components were negligible, which proved the effectiveness of using InSAR results to decompose different types of deformation (landslides, earthquakes, and ground subsidence). Second, the deformation results showed that the Beijing Plain area experienced severe ground subsidence between June 2003 and 2015, forming an obvious continuous large deformation area in the eastern and northern parts of the plain area. The spatial and temporal analysis of ground subsidence showed that since 2003, the ground subsidence rate and development trend in the plain area have presented an upward trend, from the originally sporadic subsidence areas to continuous and large-area subsidence gradually, and the value of subsidence has increased year by year. Since 2016, the deformation rate has decreased significantly, and some areas have stopped subsiding and remain stable. Third, the ground subsidence caused by InSAR was positively correlated with the decline in the groundwater level and surface water. The changes in water supply and water table depth in the plain area of Beijing from 2013 to 2015 indicated that along with the development of ground subsidence, the depth of groundwater, which was the main source of water supply, sunk year by year. The accumulated sinking displacement reached 7.42 m. At the same time, Phase I and Phase II of the SNWDP were put into operation in 2008 and 2015, respectively, which greatly reduced the demand for groundwater supply; hence, the depth of the groundwater level has continued to rebound since 2015, and the rate and development of ground subsidence have slowed down significantly. The groundwater supply continues to decrease year by year, and the utilization rate of renewable water resources is increasing year by year. It is expected that in a short period, renewable water resources will serve as the first water supply resource in Beijing.
In summary, ground subsidence mapping using the SBAS technique and multisource SAR data enabled us to manifest the spatial and temporal patterns of ground subsidence in the Beijing Plain. Therefore, it is recommended to continuously increase the observation dataset to collect and analyze the relationships between ground subsidence and groundwater or other factors in Beijing and the North China Plain and, thus, to propose more effective strategies to mitigate ground subsidence in densely populated northern areas. At the same time, further research is needed in the analysis of the relationship between land subsidence and groundwater dynamics, as well as the relationship between the evolution of subsidence and the geological background of the plain, which is also the main concern in the follow-up research.