Land Subsidence Response to Different Land Use Types and Water Resource Utilization in Beijing-Tianjin-Hebei, China

The long-term overexploitation of groundwater leads to serious land subsidence and threatens the safety of Beijing-Tianjin-Hebei (BTH). In this paper, an interferometric point target analysis (IPTA) with small baseline subset InSAR (SBAS-InSAR) technique was used to derive the land subsidence in a typical BTH area from 2012 to 2018 with 126 Radarsat-2 and 184 Sentinel-1 images. The analysis reveals that the average subsidence rate reached 118 mm/year from 2012 to 2018. Eleven subsidence features were identified: Shangzhuang, Beijing Airport, Jinzhan and Heizhuanghu in Beijing, Guangyang and Shengfang in Langfang, Wangqingtuo in Tianjin, Dongguang in Cangzhou, Jingxian and Zaoqiang in Hengshui and Julu in Xingtai. Comparing the different types of land use in subsidence feature areas, the results show that when the land-use type is relatively more complex and superimposed with residential, industrial and agricultural land, the land subsidence is relatively more significant. Moreover, the land subsidence development patterns are different in the BTH areas because of the different methods adopted for their water resource development and utilization, with an imbalance in their economic development levels. Finally, we found that the subsidence changes are consistent with groundwater level changes and there is a lag period between land subsidence and groundwater level changes of approximately two months in Beijing Airport, Jinzhan, Jingxian and Zaoqiang, of three months in Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo and Dongguang and of four months in Shengfang.


Study Area
The BTH includes Beijing, Tianjin, and the typical region of Hebei. The land area is approximately 2.18*10 5 km 2 , with a total population of approximately 110 million people. The automobile industry, electronics industry, machinery industry, and metallurgical industry are the main industries in this region, forming the main high-tech and heavy industrial base of the country. BTH is located in the northern North China Plain, which is the most extensive land subsidence area in China.
Land subsidence in the Beijing Plain first appeared in 1935; before 1952, the maximum land subsidence was only 58 mm. However, since the 1950s, land subsidence developed rapidly because of the development of industry and agriculture in the Beijing Plain area caused by the increase of groundwater exploitation and the significant decrease in groundwater level [38,39]. By 2015, the cumulative subsidence had reached 1.4 m and the area of land subsidence exceeding 100 mm had Remote Sens. 2020, 12, 457 3 of 22 reached 4000 km 2 . The areas affected by land subsidence comprise more than 60% of the Beijing Plain and form six land subsidence funnels [40][41][42][43][44].
In 1932, land subsidence first occurred in the Tianjin area. The development of land subsidence in Tianjin is divided into four main stages [45]. The first and second stages occurred from 1923 to 1957 and 1958 to 1966, respectively. Due to the low intensity of groundwater exploitation, the subsidence reached only 7.1 to 12 mm/year and 30 to 40 mm/year, respectively. However, from 1966 to 1985, the subsidence reached 80 to 100 mm/year because the amount of groundwater exploitation increased. After 1985, because the reduction in groundwater exploitation, land subsidence in the urban area was controlled at approximately 15 mm/year. Through 2015, the average subsidence was 26 mm, and the maximum subsidence was 117 mm [46]. No obvious land subsidence occurred in most parts of the north, but the land subsidence situation in most parts of the south is still grim.
Land subsidence in Hebei appeared later than that in Beijing and Tianjin. From 1950 to 1975, the land subsidence rate was less than 10 mm/year. From 1975 to 1985, land subsidence continuously increased due to the exploitation of deep groundwater in the central and eastern plains, reaching 18 to 104 mm/year [47]. Since 1986, land subsidence has developed rapidly, except for relief in Cangzhou [48]. Land subsidence occurred mainly in the middle and eastern Hebei Plain area. By 2009, there were 14 land subsidence centers, including Cangzhou, Renqiu, Suning, Hejian, Xianxian, Dongguang, Hengshui, Nangong, Pingxiang, Fengnan, Tanghai, Langfang, Baoding and Handan. By the end of 2009, the cumulative land subsidence of the Cangzhou subsidence center had reached 2580 mm, and the area with cumulative land subsidence greater than 1000 to 2000 mm had exceeded 11.3 km 2 [49].

Data Source
In this paper, we applied 126 Radarsat-2 and 184 Sentinel-1 images to derive land subsidence information in the BTH area. Radarsat-2 is an Earth observation satellite which was successfully launched on 14 December 2007, with a revisiting time of 24 days. Radarsat-2 operates in the C-band with the HH polarization mode and a wavelength of 5.6 cm. Sentinel-1 as the first satellite developed by the European Commission (EC) and the European Space Agency (ESA) for the Copernicus global earth observation project which launched in April 2014. This satellite carries a C-band SAR, with a revisit period of 12 days. For this paper, we chose an Interferometric Wide swath (IW) imaging mode with a medium resolution (5 m × 20 m). The location and information of the Radarsat-2 and Sentinel-1 images are shown in Figure 1a and Table 1. The track number of Radarsat-2 and Sentinel-1 is 26 for descending track and 142 for the ascending track, respectively. The external DEM data we selected are the Shuttle Radar Topography Mission DEM (SRTM DEM) data with 90 m resolution. Table 1. Information of the Radarsat-2 and Sentinel-1 images. R2 and S1 are the abbreviations of Radarsat-2 and Sentinel-1, respectively. The numbers in S1 are the frame numbers of the master image. "Numbers" indicates the numbers of SAR images. Furthermore, we chose the groundwater monitoring wells near the subsidence features, mentioned in Section 5.2 to compare subsidence changes with groundwater level changes. Because we only had permission to obtain the groundwater level before 2015, only the subsidence obtained by Radarsat-2 from 2012 to 2015 is included in this analysis. The information of groundwater monitoring wells is shown in Table 2.  Table 1. Information of the Radarsat-2 and Sentinel-1 images. R2 and S1 are the abbreviations of Radarsat-2 and Sentinel-1, respectively. The numbers in S1 are the frame numbers of the master image. "Numbers" indicates the numbers of SAR images.

IPTA
IPTA is a method developed by the GAMMA Company to acquire time series deformation from point targets that have stable spectral characteristics or high backscattering characteristics for single-look complex (SLC) data. Although the IPTA method can accurately obtain the results of nonlinear deformation information in a small area, its ability is often limited when processing large-scale data. However, the SBAS-InSAR technique can reduce the effects of spatio-temporal decorrelation by generating interferograms pairs based on a spatio-temporal baseline and the Doppler centroid frequency shift threshold [22]. Considering the respective advantages of the SBAS-InSAR and IPTA methods, this paper adds SBAS-InSAR into the IPTA method to process SAR images. The specifics steps of the algorithm are as follows: (1) Master image selection The master image is selected on the basis of the temporal and spatial baseline, and the Doppler centroid frequency, and its Doppler center is as close as possible to the average Doppler center of the other SAR images. In this paper, the master images for the three Radarsat-2 data and Sentinel-1 data were taken on 23 April 2014 and 25 June 2017, respectively.
(2) SAR images registration All SAR images are registered using the same radar coordinates. N interferograms are formed based on N + 1 scenes of SAR data in the SLC format, where one is chosen as a master image, and the other images are registered to the master image. This registration includes two steps, which include rough registration based on an image parameter file and precise registration based on the least square method.
(3) Small baseline differential interferograms generation We selected interferometric image pairs with spatial and temporal baseline shorter than 300 m and 300 days, respectively. The final interferometric image pairs are showed in Figure 2. The final interference pairs for R2-a, R2-b and R2-c are 205, 284 and 233, respectively, and the final interference pairs for S1-116, S1-121 and S1-126 are 177, 174, 183, respectively. The external digital elevation model (DEM) data of SRTM DEM with a 90 m resolution is used to remove the topographic phase signature. The differential interferometric phase is as follows: where ϕ dint is the differential interferometric phase, ϕ de f is the deformation phase including the line deformation and non-line deformation phases, ϕ topo is the topographic phase, ϕ atm is the atmospheric phase, ϕ orbit is the orbital phase, and ϕ noise is the noise phase.
(4) Interferometric point target extraction Coherent point target extraction is the foundation and premise of interferometric point analysis technology. The point target extraction method we used is the coherence coefficient method which refers to setting standards to extract points with coherence information [50]. The coherence coefficient was set to 0.3 in this paper.
(5) Interferometric point analysis Based on the extracted interferometric point targets, we performed a regression analysis to acquire the linear deformation and residual phases. The residual phase includes the atmospheric, non-linear deformation, and noise phases [51]. The non-linear phase shows a high correlation in the spatial domain, while the atmospheric delay shows a high correlation in the spatial domain and low correlation in the temporal domain. Noise shows a low correlation in both the spatial and temporal domains. First, the low-pass filtering in the spatial domain is used to suppress the noise phase from the residual. Secondly, in the temporal domain, the low-pass filtering is again used for spatially filtered residual phases to identify the non-linear phase. Thirdly, we superpose the non-linear phase to the linear (4) Interferometric point target extraction Coherent point target extraction is the foundation and premise of interferometric point analysis technology. The point target extraction method we used is the coherence coefficient method which refers to setting standards to extract points with coherence information [50]. The coherence coefficient was set to 0.3 in this paper.
(5) Interferometric point analysis Based on the extracted interferometric point targets, we performed a regression analysis to acquire the linear deformation and residual phases. The residual phase includes the atmospheric, non-linear deformation, and noise phases [51]. The non-linear phase shows a high correlation in the spatial domain, while the atmospheric delay shows a high correlation in the spatial domain and low correlation in the temporal domain. Noise shows a low correlation in both the spatial and temporal domains. First, the low-pass filtering in the spatial domain is used to suppress the noise phase from the residual. Secondly, in the temporal domain, the low-pass filtering is again used for spatially filtered residual phases to identify the non-linear phase. Thirdly, we superpose the non-linear phase to the linear phase to acquire the complete deformation phase. Finally, phase unwrapping is used to acquire the deformation with point targets.

Merging the InSAR Monitoring Results
(1) Merging the InSAR results in space In this study, we only consider the difference of deformation in point targets, caused by a different selections of reference points [41]. We first acquired the deformation for each SAR image with the IPTA technique. Secondly, we extracted the same interferometric point targets in

Merging the InSAR Monitoring Results
(1) Merging the InSAR results in space In this study, we only consider the difference of deformation in point targets, caused by a different selections of reference points [41]. We first acquired the deformation for each SAR image with the IPTA technique. Secondly, we extracted the same interferometric point targets in overlapping regions. Thirdly, the adjacent strip data results are adjusted to splice different strip data results. The calculation method for the difference is shown in Equation (2).
where n represents the number of point targets in the overlapping area, and a represents the average difference between two strip data results. For the Radarsat-2 data, first, we let a 1i be the subsidence at i point targets of Radarsat-2a, and let a 2i be the subsidence at the i point targets of Radarsat-2b. Then, we use Equation (2)  results of Radarsat-2a and Radarsat-2b with those of Radarsat-2c. Afterwards, we merge the InSAR results of S1-116, S1-121, and S1-126 using the same steps.
(2) Merging the InSAR results in a time-series First, the InSAR results derived by the Radarsat-2 and Sentinel-1 images are transformed from line of sight (LOS) to the vertical direction by Equation (3) to ensure the InSAR results have a common direction: where θ is the central incident angle, and d LOS is the deformation in the LOS direction. Secondly, the InSAR results of Rasarsat-2 were resampled based on the InSAR results of Sentinnel-1. Thirdly, the mean difference between the InSAR results of Rasarsat-2 and Sentinel-1 is carried out in overlapping observation periods from January 2016 to October 2016. Then, the corresponding time-series InSAR result for Sentinnel-1 is resampled based on the InSAR result of Rasarsat-2 with the previously-calculated difference. Finally, the long time-series subsidence from 2012 to 2018 in BTH is derived.

Impulse Responses Function
The impulse response function (IRF) method is the corresponding innovation accounting method in the vector autoregressive (VAR) model which was proposed by Sims in 1980 [52]. The IRF is used to measure the impact of a standard deviation shock from the random disturbance term (innovation) on the current and future values of variables. This model can intuitively describe the dynamic interaction between variables and their effects [53,54]. In this paper, the IRF method was employed to estimate the dynamic effect of groundwater level changes on land subsidence changes. The IRF curve is drawn to analyze the impact of groundwater level changes to land subsidence changes. The period corresponding to the peak value of the curve is the lag time of land subsidence. It is necessary to test the stability of the time series data before the IRF because this method is used for stationary data. In our case, an augmented Dickey-Fuller (ADF) test method was used for test, and the confidence coefficient was 5%. The test results are shown in Table 3. For the non-stationary time series values in the original time series data, a first-order difference is carried out, and the ADF test is carried out again for the sequence after this difference to ensure that the sequences used for the IRF analysis are all stationary sequences. The p-value is used to test a statistical hypothesis. The D(p-value) refers to the p-value of an ADF test based on the difference of the time series data. In Table 3, when the p-value or D(p-value) is less than 0.05, the time series data have passed the ADF test.  Figure 3 shows the distribution of land subsidence in a part of the BTH plain, which includes Radarsat-2 and Sentinel-1 data overlapped with the average land subsidence rate from 2012 to 2018. The land subsidence in BTH is very serious, and eleven subsidence features have been identified in Beijing, Langfang, Tianjin, Baoding, Hengshui, and Xingtai. We will investigate these subsidence features in the following sections. Figure 3 shows the distribution of land subsidence in a part of the BTH plain, which includes Radarsat-2 and Sentinel-1 data overlapped with the average land subsidence rate from 2012 to 2018. The land subsidence in BTH is very serious, and eleven subsidence features have been identified in Beijing, Langfang, Tianjin, Baoding, Hengshui, and Xingtai. We will investigate these subsidence features in the following sections.

Land Subsidence in Beijing
In Figure 4, we find that the spatial distribution of the land subsidence rates in Beijing is quite different. The annual subsidence increased from 2012 to 2013, during which the maximum subsidence increased from 120 mm to 148 mm. The maximum subsidence decreased from 148 mm in 2016 to 106 mm in 2018. Furthermore, the area of land subsidence exceeding 50 mm (Figure 5a) is calculated and four subsidence features (Figure 5b) are identified. From 2012 to 2014, the area of land subsidence over 50 mm shows an increasing trend from 308 km 2 to 374 km 2 even though the maximum subsidence decreased from 148 in 2013 to 131 in 2014. After 2014, the area of land subsidence over 50 mm and the annual subsidence especially show a decreasing trend. As seen in Figure 5b, the subsidence in the Jinzhan and Heizhuanghu area shows a similar trend, which increased from 2012 to 2013 and decreased from 2013 to 2018. Significantly, this subsidence was maintained at approximately 40 mm from 2012 to 2018 in the Beijing Airport area, so we should pay attention to its development.
land subsidence over 50 mm shows an increasing trend from 308 km to 374 km even though the maximum subsidence decreased from 148 in 2013 to 131 in 2014. After 2014, the area of land subsidence over 50 mm and the annual subsidence especially show a decreasing trend. As seen in Figure 5(b), the subsidence in the Jinzhan and Heizhuanghu area shows a similar trend, which increased from 2012 to 2013 and decreased from 2013 to 2018. Significantly, this subsidence was maintained at approximately 40 mm from 2012 to 2018 in the Beijing Airport area, so we should pay attention to its development.

Land Subsidence in Langfang, Tianjin
In Tianjin and Langfang, the maximum subsidence is relatively severe, with 145 mm compared to the subsidence in the Beijing area. In this typical region, three subsidence features were identified: Guangyang and Shengfang in Langfang and Wanqingtuo in Tianjin. As seen in Figure 6, we found   subsidence in this area slowed from 2012 to 2018. Meanwhile, we figure out the area of land subsidence greater than 50 mm, as shown in Figure 9(a). The result is similar to that for Beijing; the area of land subsidence greater than 50 mm decreased from 2012 to 2015 and from 2017 to 2018, increased from 2016 to 2017, with the same trend of subsidence. As seen in Figure 9(b), we found that the land subsidence in Dongguang, Jingxian, Zaoqiang and Julu experienced an overall decreasing trend from 2012 to 2018, especially after 2015.

Accuracy Assessment of InSAR Results vs Leveling Data
To evaluate the accuracy between the InSAR-derived subsidence from the Radarsat-2 and Sentinel-1 data and the leveling data, we converted the InSAR results into vertical displacements using a trigonometric equation and neglected the horizontal component of movement. In this paper, and two from 2012 to 2018 to verify the time series subsidence. The locations of these leveling benchmarks can be found in Figure 1. The precise subsidence was measured at leveling benchmarks from 2012 to 2013 and 2017 to 2018, and we selected the point targets within 200 m of the leveling benchmarks. The results of this comparison are shown in Figure 10. The differences between the two measuring methods vary from 1 mm to 23 mm and from 1 to 20 mm; the standard deviation is 6 mm and 5 mm, respectively. The results of the two leveling benchmarks from 2012 to 2018 were used to further evaluate the accuracy of InSAR-derived time series subsidence and leveling results ( Figure 11). We found that the time series subsidence of two methods is in good agreement. The mean square errors (MSE) was 8 mm and 7 mm, respectively. Meanwhile, the linear regression analysis between InSAR and leveling results in 2012 to 2013 and 2017 to 2018 were shown that the correlation coefficient (R 2 , at a 95% confidence level) reached 0.97 and 0.95, which indicated that InSAR results were highly correlated with leveling values.

Land Subsidence in Langfang, Tianjin
In Tianjin and Langfang, the maximum subsidence is relatively severe, with 145 mm compared to the subsidence in the Beijing area. In this typical region, three subsidence features were identified: Guangyang and Shengfang in Langfang and Wanqingtuo in Tianjin. As seen in Figure 6, we found that the land subsidence trend in Guangyang significantly decreases after 2015. However, the land subsidence in Shengfang and Wangqingtuo shows a decreasing trend, although that trend is not obvious. Furthermore, the time-series characteristics of land subsidence and the area exceeding 50 mm in this region are investigated. For 2012 to 2018, the area of land subsidence greater than 50 mm shows a decreasing trend from 344 km 2 to 147 km 2 . The area of land subsidence greater than 50 mm dropped by almost half from 2015 to 2018, despite a rebound in 2017. From the three subsidence features (Figure 7b), we found that in Guangyang, land subsidence decreased obviously from 60 mm in 2012 to 23 mm in 2018, and the distribution of land subsidence also decreased ( Figure 6). However, in Shengfang, the land subsidence increased from 101 mm in 2012 to 105 mm in 2013 and decreased from 105 in 2013 to 77 in 2018. Nevertheless, although the land subsidence decreased, but the subsidence range in space did not decrease ( Figure 6). In Wangqingtuo, the land subsidence was largely maintained at around 70 mm, except for 2014 (86 mm). Further, there was basically no change in the land subsidence range in space.

Land Subsidence in Heibei (Baoding, Hengshui, Xingtai)
The last subsidence bowl and four subsidence features occur in Baoding, Hengshui, Canzhou, and Xingtai of Heibei province. The maximum subsidence is 131 mm in this area, which is relatively lower compared to the other two areas (Sections 4.1 and 4.2). From Figure 8, we found that the subsidence in this area slowed from 2012 to 2018. Meanwhile, we figure out the area of land subsidence greater than 50 mm, as shown in Figure 9a. The result is similar to that for Beijing; the area of land subsidence greater than 50 mm decreased from 2012 to 2015 and from 2017 to 2018, increased from 2016 to 2017, with the same trend of subsidence. As seen in Figure 9b, we found that the land subsidence in Dongguang, Jingxian, Zaoqiang and Julu experienced an overall decreasing trend from 2012 to 2018, especially after 2015.

Accuracy Assessment of InSAR Results vs Leveling Data
To evaluate the accuracy between the InSAR-derived subsidence from the Radarsat-2 and Sentinel-1 data and the leveling data, we converted the InSAR results into vertical displacements using a trigonometric equation and neglected the horizontal component of movement. In this paper, there are total of 72 leveling benchmarks: 35 from 2012 to 2013 to verify the subsidence derived from  Figure 10. The differences between the two measuring methods vary from 1 mm to 23 mm and from 1 to 20 mm; the standard deviation is 6 mm and 5 mm, respectively.
The results of the two leveling benchmarks from 2012 to 2018 were used to further evaluate the accuracy of InSAR-derived time series subsidence and leveling results (Figure 11). We found that the time series subsidence of two methods is in good agreement. The mean square errors (MSE) was 8 mm and 7 mm, respectively. Meanwhile, the linear regression analysis between InSAR and leveling results in 2012 to 2013 and 2017 to 2018 were shown that the correlation coefficient (R 2 , at a 95% confidence level) reached 0.97 and 0.95, which indicated that InSAR results were highly correlated with leveling values.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22 Figure 11. (a) and (b) Comparison between time series subsidence derived by InSAR and leveling results in two levelings, respectively, whose locations are shown in Figure 1 with levelings-I.

Land Subsidence Response to Different Land Use Types
As seen in Figure 12 and Table 4, the land-use in the subsidence features of Shangzhuang, Beijing airport, Jinzhan, Heizhuanghu, Wangqingtuo and Shengfang is relatively complex. The maximum subsidence in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang areas exceeds 100 mm/year, with 117, 105, 120 and 101 mm/year, respectively. In the Jinzhan and Heizhuanghu areas, the dense residential land and high population density in the area lead to a considerable water demand. Meanwhile, these areas are located in the groundwater funnel area of Beijing, which leads to large land subsidence. Since the MSWDP began operation in 2015, the urban water supply has been alleviated, and land subsidence has been slowed to some degree. As for Wangqingtuo and Shenfang areas, the dense population as well as agricultural and industrial lands there all result in an increasing demand of water and further aggravation of land subsidence. A large number of water-consuming plants and crops are planted in this region, and 85% of the well water is used for residents' lives and agricultural production [55]. Unlike the Wangqingtuo area, Shengfang mainly focuses on industry, mainly the steel industry, glass industry, and sheet metal industry. The steel industry, in particular, involves high water consumption industry. There are three large steel mills in the region, with a steelmaking capacity of 10 million tons. Industries with high water consumption may lead to the serious development of land subsidence in this area. Similarly to Shengfang, in the Guangyang area, the main industrial categories are machinery manufacturing, building materials, electronic appliances, food, seasoning, clothing, glass products, printing, and chemical. In addition, the agricultural water demand is also large (winter wheat, maize, peanut, and cotton are mainly planted here). The crop area of Jinzhou town of Guangyang is one of the top 1000 crop sown areas in China. Figure 11. (a) and (b) Comparison between time series subsidence derived by InSAR and leveling results in two levelings, respectively, whose locations are shown in Figure 1 with levelings-I.

Land Subsidence Response to Different Land Use Types
As seen in Figure 12 and Table 4, the land-use in the subsidence features of Shangzhuang, Beijing airport, Jinzhan, Heizhuanghu, Wangqingtuo and Shengfang is relatively complex. The maximum subsidence in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang areas exceeds 100 mm/year, with 117, 105, 120 and 101 mm/year, respectively. In the Jinzhan and Heizhuanghu areas, the dense residential land and high population density in the area lead to a considerable water demand. Meanwhile, these areas are located in the groundwater funnel area of Beijing, which leads to large land subsidence. Since the MSWDP began operation in 2015, the urban water supply has been alleviated, and land subsidence has been slowed to some degree. As for Wangqingtuo and Shenfang areas, the dense population as well as agricultural and industrial lands there all result in an increasing demand of water and further aggravation of land subsidence. A large number of water-consuming plants and crops are planted in this region, and 85% of the well water is used for residents' lives and agricultural production [55]. Unlike the Wangqingtuo area, Shengfang mainly focuses on industry, mainly the steel industry, glass industry, and sheet metal industry. The steel industry, in particular, involves high water consumption industry. There are three large steel mills in the region, with a steelmaking capacity of 10 million tons. Industries with high water consumption may lead to the serious development of land subsidence in this area. Similarly to Shengfang, in the Guangyang area, the main industrial categories are machinery manufacturing, building materials, electronic appliances, food, seasoning, clothing, glass products, printing, and chemical. In addition, the agricultural water demand is also large (winter wheat, maize, peanut, and cotton are mainly planted here). The crop area of Jinzhou town of Guangyang is one of the top 1000 crop sown areas in China. Dongguang, Jingxian, Zaoqiang, and Julu are the main grain production areas in Hebei. The main crops in the region are wheat and maize, which feature the largest total water consumption in BTH, with the areas of 680, 1293, 753, and 607 million km 2 , respectively [56]. The annual groundwater overdraft for crop irrigation ranges from 3 billion to 5.5 billion m 3 [57]. A large amount of groundwater overexploitation results in serious land subsidence in the region. The serious land subsidence could affect the urban public facilities, transportation and buildings, such as the damage of urban underground pipe network, the instability of railway subgrade and cracks of buildings. The response characteristics of groundwater change and land subsidence are discussed in Section 5.2.  Expressway, subway, and urban residential land 9-56

Land Subsidence Response to Different Water Use
As seen in Figure 13, we found that the Beijing and Hebei supplied water mainly from the groundwater, while Tianjin mainly supplied surface water. In Beijing and Hebei, the total groundwater supply amount was 2.04 and 15.13 to 1.63 and 10.61 billion m 3 from 2012 to 2018, and the total surface water supply amount was 0.8 and 4.13 to 1.23 and 7.04 billion m 3 from 2012 to 2018, respectively. However, in Tianjin, the total groundwater and surface water supply amount were 0.55 to 0.44 and 1.6 to 1.95 billion m 3 from 2012 to 2018, respectively. From the point of view of water consumption, Beijing needs more water for residential use, while Tianjin and Hebei need more water for agricultural use, especially Hebei, which used more than 10 billion m 3 water for agricultural use from 2012 to 2018. The water for residential use in BTH showed an increasing trend, which ranged from 1.6 to 1.84 billion m 3 in Beijing, 0.5 to 0.74 billion m 3 in Tianjin, and 2.34 to 2.78 billion m 3 in Hebei, from 2012 to 2018, with an increase in population. The population of Beijing, Tianjin, and Hebei increased from 20.69, 14.13, and 72.88 million to 21.54, 15.6, and 75.56 million, respectively. Due to major water transfer projects such as the MSWDP, the surface water supply in BTH shows an increasing trend, while the groundwater supply shows a decreasing trend from 2012 to 2018. Further investigations revealed that Beijing and Tianjin received 34.9 and 33.8 million m 3 water from MSWDP as surface water for the urban water supply. In summary, after 2015, the exploitation of underground water decreased and the supply of surface water increased, which slowed the land subsidence in BTH.
Hebei increased from 20.69, 14.13, and 72.88 million to 21.54, 15.6, and 75.56 million, respectively. Due to major water transfer projects such as the MSWDP, the surface water supply in BTH shows an increasing trend, while the groundwater supply shows a decreasing trend from 2012 to 2018. Further investigations revealed that Beijing and Tianjin received 34.9 and 33.8 million m 3 water from MSWDP as surface water for the urban water supply. In summary, after 2015, the exploitation of underground water decreased and the supply of surface water increased, which slowed the land subsidence in BTH.  (c) is Beijing, Tianjin and Hebei, respectively. Total surface water supply refers to the water intake of surface water engineering, which are water storage, water extraction and water diversion. Total groundwater supply refers to the exploitation of water wells. Total industrial water refers to the water used by industrial and mining enterprises in manufacturing, processing, cooling, air conditioning, purification, washing and other aspects in the production process. Total agricultural water use includes farmland irrigation Figure 13. Water resources utilization in BTH from 2012 to 2018. (a-c) is Beijing, Tianjin and Hebei, respectively. Total surface water supply refers to the water intake of surface water engineering, which are water storage, water extraction and water diversion. Total groundwater supply refers to the exploitation of water wells. Total industrial water refers to the water used by industrial and mining enterprises in manufacturing, processing, cooling, air conditioning, purification, washing and other aspects in the production process. Total agricultural water use includes farmland irrigation water, forest and fruit land irrigation water, grassland irrigation water, fish-pond recharge water and livestock and poultry water. Total residential water use includes urban domestic water and rural domestic water.
In Figure 14, the groundwater level shows a drop with a significant fluctuation trend, especially in wells h to j. The groundwater level is the lowest around May to August in wells h to j, possibly because these wells are located in Hengshui of the Hebei area, which is one of the national grain bases. The main crop in the area is winter wheat [57]. Because the irrigation period of this crop is February and March, and the precipitation is lower at this time, the groundwater exploitation volume is large, and the groundwater level drops significantly. After June, due to an increase in precipitation, the exploitation of groundwater decreased and the water level gradually rises. Furthermore, the comparison results of land subsidence changes and groundwater level changes reveal that the subsidence changes are consistent with the groundwater level changes, showing a seasonal trend. Meanwhile, we found that there is a time lag between land subsidence and groundwater level.
Furthermore, we applied the IRF method to investigate the lag time between subsidence change and groundwater level change. The ADF test and IRF results are shown in Table 3 and Figure 14, respectively. It can be seen from Figure 15 that groundwater level change has an obvious effect on land subsidence at each subsidence feature, and this effect has a lag effect. For the subsidence features of Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo, and the Dongguang area, the groundwater level changes in the third period (the third months) have the greatest impact on land subsidence, which means that there is a lag period of approximately three months between land subsidence and groundwater level change. However, for the subsidence features of Beijing Airport, Jinzhan, Jingxian and Zaoqiang, the groundwater level change in the second period (the second months) have the greatest impact on land subsidence, which means that there is a lag period of approximately two months between land subsidence and groundwater level change. Moreover, for the Shengfang subsidence feature, the groundwater level changes in the fourth period (the fourth months) have the greatest impact on land subsidence, which means that there is a lag period of approximately four months between land subsidence and groundwater level change.
is large, and the groundwater level drops significantly. After June, due to an increase in precipitation, the exploitation of groundwater decreased and the water level gradually rises. Furthermore, the comparison results of land subsidence changes and groundwater level changes reveal that the subsidence changes are consistent with the groundwater level changes, showing a seasonal trend. Meanwhile, we found that there is a time lag between land subsidence and groundwater level. Jinzhan, Jingxian and Zaoqiang, the groundwater level change in the second period (the second months) have the greatest impact on land subsidence, which means that there is a lag period of approximately two months between land subsidence and groundwater level change. Moreover, for the Shengfang subsidence feature, the groundwater level changes in the fourth period (the fourth months) have the greatest impact on land subsidence, which means that there is a lag period of approximately four months between land subsidence and groundwater level change. Figure 15. Orthogonal impulse response from groundwater level change to subsidence in subsidence features of (a) Shangzhuang Beijing Airport, Jinzhan and Heizhuanghu, (b) Guangyang, Wangqingtuo and Shengfang, (c) Dongguang, Jingxian and Zaoqiang. The horizontal axis represents the lag order of the impact (month), and the vertical axis represents the response degree.

Conclusion
Because of the rapid growth of population and the development of industry and agriculture in BTH, the regional water demand is huge. However, the over-exploitation of groundwater in the BTH is obvious due to the shortage of water resources, which leads to the serious land subsidence in BTH. The serious subsidence in BTH resulted in the elevation loss, damage of roads, bridges, underground pipes and other municipal facilities, urban flooding and hydrops. In this study, the IPTA with SBAS-InSAR method was optimized to investigate land subsidence based on 126 Radarsat-2 and 184 Sentinel-1 images from 2012 to 2018 in BTH. We compared the InSAR-derived results with levelingderived results from 2012 to 2013 and 2017 to 2018. The correlation coefficients (at a 95% confidence level) all exceed 0.9, indicating that our result is reliable. The land subsidence in BTH is very serious, with the maximum subsidence rate exceeding 120 mm/year. Eleven subsidence features in three subsidence areas were identified to be distributed in Beijing, Tianjin and Langfang, Baoding, Hengshui, and Xingtai, with maximum subsidence rates of 148, 145 and 131 mm/year, respectively. The area of subsidence greater than 50 mm in these areas was shown a decreasing trend after 2015. Meanwhile, the analysis of time-series subsidence among eleven features also shows a decreasing trend starting from 2015.
Furthermore, we discussed the relationship between land subsidence and different land-use types in subsidence features. We found that when the land-use type in an area is residential land or overlaps with agricultural and industrial land, land subsidence is relatively serious. For example, in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang area, the maximum subsidence exceeds 100 Figure 15. Orthogonal impulse response from groundwater level change to subsidence in subsidence features of (a) Shangzhuang Beijing Airport, Jinzhan and Heizhuanghu, (b) Guangyang, Wangqingtuo and Shengfang, (c) Dongguang, Jingxian and Zaoqiang. The horizontal axis represents the lag order of the impact (month), and the vertical axis represents the response degree.

Conclusions
Because of the rapid growth of population and the development of industry and agriculture in BTH, the regional water demand is huge. However, the over-exploitation of groundwater in the BTH is obvious due to the shortage of water resources, which leads to the serious land subsidence in BTH. The serious subsidence in BTH resulted in the elevation loss, damage of roads, bridges, underground pipes and other municipal facilities, urban flooding and hydrops. In this study, the IPTA with SBAS-InSAR method was optimized to investigate land subsidence based on 126 Radarsat-2 and 184 Sentinel-1 images from 2012 to 2018 in BTH. We compared the InSAR-derived results with leveling-derived results from 2012 to 2013 and 2017 to 2018. The correlation coefficients (at a 95% confidence level) all exceed 0.9, indicating that our result is reliable. The land subsidence in BTH is very serious, with the maximum subsidence rate exceeding 120 mm/year. Eleven subsidence features in three subsidence areas were identified to be distributed in Beijing, Tianjin and Langfang, Baoding, Hengshui, and Xingtai, with maximum subsidence rates of 148, 145 and 131 mm/year, respectively. The area of subsidence greater than 50 mm in these areas was shown a decreasing trend after 2015. Meanwhile, the analysis of time-series subsidence among eleven features also shows a decreasing trend starting from 2015.
Furthermore, we discussed the relationship between land subsidence and different land-use types in subsidence features. We found that when the land-use type in an area is residential land or overlaps with agricultural and industrial land, land subsidence is relatively serious. For example, in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang area, the maximum subsidence exceeds 100 mm/year. The distribution of different land use types leads to different water use structures. There is a large number of residential areas in Jinzhan and Heizhuang, with a large population and a correspondingly large water demand. In Guangyang, Wangqingtuo, and Shengfang, many high-water consumption industries and agriculture increase the regional water demand. However, a large amount of agricultural land is distributed in the Dongguang, Jingxian, Zaoqiang and Julu areas, and most of that land is cultivated wheat. Due to the lack of precipitation during the irrigation period of wheat, a large amount of groundwater exploitation leads to land subsidence in the region.
The relationship between land subsidence changes and groundwater level changes in the subsidence features was examined by using ten groundwater monitoring wells. We determined that the land subsidence changes are consistent with the groundwater level changes and show a seasonal trend. Nevertheless, there is a time lag between land subsidence changes and groundwater levels. The IRF method was used to investigate the time lag between subsidence changes and groundwater level changes. Our analysis shows that the lag time is two months in Beijing Airport, Jinzhan, Jingxian, and Zaoqiang, three months in Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo, and Dongguang, and four months in Shengfang.