Combined GRACE and MT-InSAR to Assess the Relationship between Groundwater Storage Change and Land Subsidence in the Beijing-Tianjin-Hebei Region

Abstract

Beijing-Tianjin-Hebei (BTH) has been suffering from severe groundwater storage (GWS) consumption and land subsidence (LS) for a long period. The overexploitation of groundwater brings about severe land subsidence, which affects the safety and development of BTH. In this paper, we utilized multi-frame synthetic aperture radar datasets obtained by the Rardarsat-2 satellite to monitor land subsidence’s temporal and spatial distribution in the BTH from 2012 to 2016 based on multi-temporal interferometric synthetic aperture radar (MT-InSAR). In addition, we also employed the Gravity Recovery and Climate Experiment (GRACE) mascon datasets acquired by the Center for Space Research (CSR) and Jet Propulsion Laboratory (JPL) to obtain the GWS anomalies (GWSA) of BTH from 2003 to 2016. Then we evaluate the accuracy of the results obtained. Furthermore, we explored the relationship between the regional GWSA and the average cumulative subsidence in the BTH. The total volume change of subsidence is 59.46% of the total volume change of groundwater storage. Moreover, the long-term decreasing trend of the GWSA (14.221 mm/year) and average cumulative subsidence (17.382 mm/year) show a relatively high consistency. Finally, we analyze the heterogeneity of GWS change (GWSC) and LS change (LSC) in the four typical areas by the Lorenz curve model. The implementation of the South-to-North Water Diversion Project (MSWDP) affects the heterogeneity of GWSC and LSC. It can be seen that the largest heterogeneity of LSC lags behind the GWSC in the Tianjin-Langfang-Hengshui-Baoding area. The largest uneven subsidence in Beijing and Tianjin occurred in 2015, and the largest uneven subsidence in Hengshui-Baoding occurred in 2014. After that, the heterogeneity of subsidence gradually tends to stable.

1. Introduction

As China’s “capital circle”, BTH is one of the regions with the most rapid economic development and the highest population density in the world. The area of this region is 2.27% of the whole country, but the water resources account for only 1% [1]. With the city’s rapid development, about 70% of the water resources comes from groundwater [2]. The extraction volume has reached 20 billion m³/year [3,4], which has formed the world’s largest groundwater funnel area. Excessive exploitation of groundwater has led to a severe regional land subsidence [5,6,7]. Thus, it is essential and necessary to explore the relationship between groundwater and land subsidence in the BTH.

To monitor land subsidence, traditional methods (e.g., leveling, layer-wise mark measurements, and global positioning system) can obtain the high temporal resolution and measurement accuracy. However, it is difficult to obtain large-scale information for land subsidence monitoring because of the expensive equipment and the lower spatial resolution. The interferometric synthetic aperture radar (InSAR) technology that appeared in the 1980s made up for the shortcomings of traditional measurements, which has the characteristics of all-weather, wide monitoring range, and high monitoring accuracy [8]. In 1989, Gabriel et al. demonstrated that synthetic aperture radar differential interferometric (D-InSAR) technology could be used to detect sub-centimeter-level minute surface deformations for the first time [9]. However, D-InSAR technology is limited by the temporal and spatial decoherence [10], the evaluation of parameters in phase unwrapping [11], the influence of atmospheric phase [12]. Multitemporal synthetic aperture differential interferometry (MT-InSAR) technology can overcome the limitation of D-InSAR. With the transition of InSAR technology from D-InSAR to MT-InSAR, from the monitoring of large deformations to the monitoring of small deformations, the impact of atmospheric errors has been effectively weakened. The monitoring accuracy has been increased from cm to mm, and the technology has become more mature, making it an effective means to obtain urban land subsidence [13].

The traditional monitoring methods generally include selecting representative wells, springs, and operation wells for measurement [14]. Although it can accurately measure groundwater level changes, it is tough to obtain the changes of groundwater in the entire area with limited measurement point data due to the complicated hydrogeological conditions. The GRACE satellite issued jointly by the National Aeronautics and Space Administration (NASA) and the German Space Agency (DLR) was successfully launched in March 2002, bringing new opportunities for groundwater monitoring [15]. Many scholars have tried to use GRACE to retrieve water storage change information. They have achieved great success, which has dramatically promoted progress in water storage research and played a significant role in hydrological research. Existing research results indicate that GRACE can invert large-scale changes in surface quality, which is sufficient to reveal changes in land water storage with an average water height of less than 1 cm [16,17,18]. The changes in groundwater storage can be inverted by combining the land water storage obtained by GRACE and external data. It proves that GRACE data is adequate for monitoring the changes in groundwater storage. The method has been widely used in many countries, such as India [19,20], the United States [21,22], the Middle East [23,24], Russia [25], Canada [26], and China [4,27].

Overexploitation of groundwater will cause the water level to drop, and then the soil layer will be compressed, resulting in ground subsidence. Some scholars have experienced the transition from point scale to surface scale in the related research on groundwater and land subsidence. In the point-scale study, the connection between groundwater level at the well and land subsidence in the monitoring location have been investigated [28]. In the larger-scale study, the spatial consistency between the regional groundwater contour and the land subsidence center have been explored [29]. Further research, relevant scholars have begun to integrate the InSAR and GRACE technologies to better find out the connection between regional groundwater and land subsidence [30,31].

The MSWDP will improve the state of groundwater and ground subsidence in the BTH. Although some researchers have explored the groundwater storage and land subsidence in North China Plain [32,33]; however, previous research usually focused on describing the secular trend between groundwater storage consumption and land subsidence. Few researchers have explored the spatial relations and the seasonal monthly scale relationship between land subsidence and groundwater storage. Moreover, the characteristics of land subsidence and groundwater storage change in the BTH have not been quantitatively described.

In this paper, the MT-InSAR and GRACE was used to acquire the land subsidence (LS) and the groundwater storage anomalies (GWSA) in the BTH region. Then, we obtained and analyzed the evolution characteristics of LS and GWSA (Section 4). Furthermore, we explored the relationship between the regional GWSA and the average cumulative subsidence in the BTH (Section 5). Finally, we discussed the heterogeneity of groundwater storage changes (GWSC) and land subsidence changes (LSC) in four typical areas by the Lorenz curve model (Section 5.3).

2. Study Area and Dataset Materials

2.1. Study Area

The BTH, located in (36°05′–42°40′N, 113°27′–119°50′E), which is shown in Figure 1, belongs to a temperate monsoon climate. The climate is characterized by high temperature and is rainy in summer due to the influence of ocean water vapor; in winter, it is cold and dry due to the effect of cold inland air. From the foothills to the coast, the BTH Plain can be divided into quaternary piedmont alluvium, alluvial slope plain (piedmont plain), central alluvial, and lacustrine multilayer plain (central plain), and eastern coastal plains (coastal plains) with alluvial and lacustrine interspersed with the sea [34]. The groundwater storage in the piedmont plain has strong permeability and is the primary infiltration and replenishment area for regional groundwater. The central plain is distributed with a continuous impermeable layer with a large thickness, resulting in its pressure-bearing characteristics and poor ability to receive atmospheric precipitation. The coastal plain has poor water permeability, water storage, and water supply capabilities [35]. The piedmont plains and the central plains have slow water circulation. The quaternary loose sediments are distributed widely in the plain area, with a thickness ranging from 150 to 600 m [36]. As China’s “capital circle,” the BTH is one of the regions with the most rapid economic development and the highest population density in the world. However, the area covers less than 2.27% of the country’s area, and the whole water resources account for only 1%, the population accounts for 8% of the country [1]. Because of the deficit of water resources, the BTH region has long been exploiting excessive groundwater to maintain economic development and human water use. Over-utilized groundwater has caused severe land subsidence in BTH [2,37,38].

In 1923, Tianjin discovered a decline in the level of elevation in the urban area. It became the earliest area where the groundwater was mined to produce land subsidence in the BTH region [39]. After that, subsidence occurred in Beijing and Hebei area one after another. Before 1950s, land subsidence occurred in some areas of the BTH region, and the amount of land subsidence was relatively weak. However, the groundwater of BTH has been used as a major source of water since 1970; thus, the groundwater level has dropped rapidly due to the excessive exploitation [40]. Due to the complex influencing factors of land subsidence, the exploitation and utilization of groundwater vary from place to place, and the development rate and influence range of land subsidence have significant differences in different spaces and times.

From the 1960s to the 1980s, the land subsidence of the BTH entered a period of rapid development. Two subsidence centers, Dajiaoting and Laiguangying, were formed in the eastern suburbs of the Beijing Plain, and the subsidence range gradually expanded. Tianjin’s rapid growth of land subsidence has formed subsidence centers, such as the Central City, Tanggu District, Hangu District, and the Lower Haihe River Industrial Zone. Due to the large-scale exploitation of deep groundwater in Hebei’s east and central plains, a large area of deep groundwater fall funnel has appeared, and the scope of ground subsidence has continued to expand [41]. From the 1980s to the 1990s, due to the different groundwater exploitation management measures in different regions, there were regional differences in the development of land subsidence. The scope of the groundwater funnel in the Beijing Plain has further expanded. The subsidence area has migrated to the north and south, initially forming five main subsidence centers. In the context of Tianjin’s government vigorously controlling subsidence, the rate of land subsidence has shown a slowing trend in the Central City and Binhai New Area. With the economic development of Hebei Province, the population has overgrown, and the amount of groundwater extraction has also increased significantly. The land subsidence area has continued to expand, and it has shown a trend of regional development [38]. Since the 21st century, the regional differences in land subsidence development are still quite significant. The Beijing Plain has formed two major subsidence areas, one north and one south, and seven subsidence centers [42]. The scope of the two subsidence areas continues to expand outward. The land subsidence in Tianjin has entered the stage of comprehensive treatment. The subsidence of the land subsidence centers in Wuqing, the Central urban area, and Tanggu has been effectively controlled. Still, the land subsidence around the metropolitan area and suburban counties have developed rapidly [43]. Except for the relief of land subsidence in Cangzhou City, Hebei Plain, land subsidence in other areas is still developing rapidly [44].

From the perspective of the overall development of land subsidence in the North China Plain, it is still in a rapid development stage. The BTH region presents the characteristics of inter-regional contiguous distribution, and the subsidence rate is still relatively large.

2.2. Dataset Materials

In this research, 126 Rardarsat-2 images were acquired from the descending orbit between January 2012 to October 2016. The Rardarsat-2 satellite has a revisiting cycle of 24 days and the operation mode is at C-band with VV polarization and a 5.6 cm wavelength. The spatial resolution is 30 m, and the single scene image coverage area is 150 km × 150 km. In order to cover the BTH region, we used three sets of frame images. Besides, we chose the Shuttle Radar Topography Mission data as the external DEM.

In addition, Monthly GRACE Release-06 mascon products (2003–2016) from the Center for Space Research (CSR) and Jet Propulsion Laboratory (JPL) were employed. GRACE mascon can obtain the total terrestrial water storage (TWS) anomalies. In order to estimate GWS anomalies, the averaged soil moisture storage (SMS) and the snow water equivalent storage (SWES) simulated by Noah models from the Global Land Data Assimilation System (GLDAS) were subtracted from TWS anomalies (TWSA).

Moreover, we collected the BTH water resources from the Beijing Water Resources Bulletin, Tianjin Water Resources Bulletin, and Hebei Water Resources Bulletin, including total water supply, groundwater supply, and MSWDP supply (as shown in Table 1). Plus, we also obtained monthly 0.5° × 0.5° precipitation data during 2012–2016 from China Meteorological Administration.

3. Methods

3.1. MT-InSAR

In order to solve the problem about the spatial-temporal decorrelation and atmospheric delay of the D-InSAR, Ferretti et al. proposed the PS-InSAR method in 2001 [45]. As, a kind of MT-InSAR, PS-InSAR method, has the advantages to obtain the land subsidence of BTH.

In the beginning of the PS-InSAR processing, the SAR images covering the study area were processed by D-InSAR to form a time-series interference image pair.

Each pixel contains the following components after D-InSAR processing:

(1)${\mathsf{\Phi }}_{InSAR}={\mathsf{\Phi }}_{def}+{\mathsf{\Phi }}_{top}+{\mathsf{\Phi }}_{atm}+{\mathsf{\Phi }}_{oribit}+{\mathsf{\Phi }}_{noise}$

Among them: ${\mathsf{\Phi }}_{InSAR}$ is the interference phase of the point target; ${\mathsf{\Phi }}_{def}$ is the deformation phase of the radar line-of-sight direction; ${\mathsf{\Phi }}_{top}$ is the terrain phase; ${\mathsf{\Phi }}_{atm}$ is the atmospheric delay phase; ${\mathsf{\Phi }}_{oribit}$ is the orbit error phase; ${\mathsf{\Phi }}_{noise}$ is the noise phase.

Then, we selected the PS point, which has strong and stable scattering characteristics, such as buildings with dihedral angle, road edges, bridges, and exposed rocks. Then the extracted PS point phase information was filtered to remove the influence of the atmospheric phase and acquire the deformation phase information to obtain the land subsidence rate.

Subsequently, merging the InSAR results in space. There are three sets of frame data covering the BTH. In order to better analyze the results, it is necessary to combine the three groups to obtain the overall subsidence results by using the point averaging method of the overlapping area. The specific calculation equation is as follows:

(2)$\overline{x}=\frac{{\displaystyle \sum _{i=1}^n(x_{1i}-x_{2i})}}{n}$

where n is the number of PS point in the overlapping area, and the $\overline{x}$ represents the average different between different frame data. We merge for the three sets of frame data according to Equation (2) from top to bottom.

Finally, it is essential and meaningful to convert the line-of-sight deformation to vertical for more accurate research and analysis [46].

(3)$d_u=\frac{d_{los}}{\cos \theta }$

where θ is the central incident angle of the SAR satellite, and $d_{los}$ is the deformation in the line-of-sight direction.

3.2. GRACE Data Analysis

Based on satellite observations of time-varying gravity, GRACE provides important ideas and possibilities for changes in global water storage [47]. Relevant studies have shown that the changes in groundwater storage can be accurately extracted from total terrestrial water storage combined with external auxiliary data [48,49].

In this study, we used the monthly GRACE release-06 mascon products from the CSR and JPL. GRACE mascon solutions can acquire the total terrestrial water storage (TWS) anomalies. GRACE data is the total change in TWS, including all the surface water, snow, soil water, canopy water content, and groundwater in a region. The groundwater storage (GWS) can be obtained when combined with auxiliary hydrological datasets. Considering the arid climate of the North China Plain, there is little surface water. About 40% of the rivers are dry throughout the year, and the surface water storage of rivers, lakes, and reservoirs is negligible [50,51]. In addition, the amount of canopy water obtained is very small, we ignored its influence in our research. Hence, we ignore the surface water and canopy water when calculating changes in groundwater storage.

(4)$\mathsf{\Delta }GWS=\mathsf{\Delta }TWS-\left(\mathsf{\Delta }SMS+\mathsf{\Delta }SWES\right)$

where ΔGWS is the change of groundwater storage, ΔTWS is the change of whole terrestrial water storage, ΔSMS is the change of soil moisture on the surface, and ΔSWES is the change of snow water equivalent storage. The 0–10 cm, 10–40 cm, 40–100 m, and 100–200 cm SMS and SWES were provided by the 0.25° × 0.25° GLDAS-Noah hydrological model.

3.3. Spatio-Temporal Data Interpolation

We resampled the GLDAS data and the GRACE data using spatial kriging interpolation. The purpose is to ensure the consistency of spatial grid position (0.5° × 0.5°) between the TWS change acquired by GRACE data and the change of SMS and SWES acquired by GLDAS data. In the spatial trend analysis, we resample the grid points obtained by GRACE and the PS points obtained by InSAR into 2 km × 2 km raster pixel images, respectively. In addition, to better study the seasonal changes in groundwater storage and land subsidence from 2013 to 2016 and ensure the same time resolution (1 month), the cubic-spline interpolation method was used to interpolate some missing months of land subsidence.

3.4. Loren Curve Model

The Lorenz curve is a method proposed by the American economic statistician M. Lorenz to describe the unequal distribution of social wealth. The Lorenz curve can be expressed mathematically as [52]:

(5)$\begin{array}{l}L(y)=\frac{{\displaystyle {\int }_0^{y}xdF(x)}}{\mu }\\ \end{array}$

where F(x) is the cumulative distribution function of ordered individuals, and $\mu$ is the mean. The Lorenz curve uses a graph to visually express wealth concentration, indicating the percentage of accumulated wealth in a given percentage of the population. Under ideal conditions, the relationship between wealth and population is 1:1. The degree of curvature of the Lorenz curve is of great significance, the greater the degree of curvature, the more unequal income distribution. The Lorenz curve is shown in Figure 2.

The Gini coefficient was proposed by the Italian economist Corrado Gini in 1922 based on the Lorenz curve. It uses the curve drawn by the cumulative number to describe the degree of inequality (centralization or dispersion). The Gini coefficient is the ratio of the area between y = x and the Lorenz curve in the Lorenz curve (A) and the area between the complete uniform line and the whole uneven line (A + B). The Gini coefficient is expressed by Gini and is obtained by the integral method:

(6)$\mathrm{Gini}=1-2{\displaystyle {\int }_0^{1}f(x)dx}$

The variation range of the Gini coefficient is 0 ≤ Gini ≤ 1. When Gini = 0, all individuals are equal; as the Gini coefficient increases, the research individuals become more uneven. The Gini coefficient is used to evaluate and quantify the uniformity of distribution. In addition to being used in economics and social income distribution, it can also be applied to other disciplines to evaluate the evenness of distribution. For example, it has been successfully applied to analyze uneven precipitation distribution [53,54] and temporal inhomogeneity of runoff [55], When it is applied to heterogeneity of the land subsidence and groundwater storage change, when Gini = 0, there is no difference in changes in various periods during the year. As the Gini coefficient increases, the degree of difference in changes is significant, the greater the heterogeneity.

3.5. Technical Flow Chart

The entire flowchart of this study is shown in Figure 3:

Firstly, the PS-InSAR method was adopted to process RADARSAT-2. We validated the InSAR results by leveling data.

Secondly, GRACE mascon products were adopted to obtain the total terrestrial water storage. Then, the soil moisture storage (SMS) and the snow water equivalent storage (SWES) were simulated by GLDAS Noah models. Furthermore, we obtained the groundwater storage anomalies (GWSA) by combining TWS anomalies (TWSA) and the SMS and SWES. As a result, we obtained the GWS anomalies.

Thirdly, temporal and spatial evolution characteristics and correlation between land subsidence and groundwater storage in BTH were explored. On the one hand, the long-term trends of cumulative subsidence and groundwater storage were analyzed first, and then the volume of land subsidence and groundwater storage changes in different regions were calculated. On the other hand, the short-term seasonal changes between the land subsidence monthly changes and groundwater storage were analyzed. Moreover, the heterogeneity of land subsidence changes and groundwater storage changes was obtained by the Lorenz curve model.

4. Results

4.1. Land Subsidence in BTH

4.1.1. Spatial Distribution and Evolution Characteristics of Land Subsidence in BTH

Land subsidence is distributed widely, and the spatial difference is obvious in BTH. The average subsidence rate in part of the BTH Plain from 2012 to 2016 is shown in Figure 4. The maximum average subsidence rate is 143.35 mm/year, and the maximum accumulated subsidence is up to 708 mm (as shown in Figure 5), which is located in Jinzhan, Chaoyang District, Beijing. In 2012–2013, the maximum subsidence of the study area increased from −122.63 mm to −156.91 mm, and it started to decrease from −148.86 mm to −123.88 mm after 2014. The maximum alarm value of land subsidence is 50 mm according to the National Norms. The part of the study area experiencing an average subsidence rate greater than 50 mm/year amounts to 1610 km². The area subsided more than 50 mm/year was expanded from 1804 km² to 2556 km² from 2012 to 2013. Since then, the subsidence area has been continuously reduced after MSWDP (as shown in Figure 6). The largest subsided area in Tianjin is located near Wangqingtuo Town in Wuqing District, and its location is marked with a black triangle in Figure 4, with an average subsidence rate of 140 mm/year. In addition, there are large-scale subsidence zones connected to the southeast of Baoding, northern Hengshui, and western Cangzhou in the BTH Plain. The rate is mostly larger than 50 mm/year.

4.1.2. Accuracy Verification

To assess the accuracy of the subsidence information from the Radarsat-2, we used the 18 known benchmarking points from 2015 to 2016, the locations of which can be found in Figure 1. We transformed the deformation of the line-of-sight of the InSAR radar to vertical and neglected the horizontal component of movement. Among them, we used the benchmarking point as the center, using ArcGIS for buffer analysis, obtaining the 200-m buffer range of all benchmarking points, extracting the PS points in the buffer, respectively. Finally, the PS point values of the respective buffer ranges were averaged and then verified with the corresponding benchmark points. As a result, a higher correlation value was obtained, with a Pearson correlation coefficient of 0.97 and a Spearman correlation coefficient of 0.95 (Figure 7). Thus, it can be suggested that the results of land subsidence are accurate and usable.

4.2. Groundwater Storage Change in BTH

In this paper, the CSR and JPL data from 2003 to 2016 are processed, and the two data results are averaged to obtain the final abnormal value of groundwater storage in BTH. The difference of the CSR and JPL is mainly due to the different data processing procedures of the two official institutions [56,57]. It can be seen that the TWSA results from CSR and JPL in BTH are roughly consistent and have same fluctuation from Figure 8. Next, the SMSA and SWESA were obtained from GLDAS. The specific and detailed timing information is shown in Figure 9. Then we obtained the GWSA (As shown in Figure 10) by combining the results of GRACE and GLDAS. The two different data have the effect of mutual verification, and we can better use the average value to ensure the accuracy of the overall data results. The GWSA derived by GRACE and in situ groundwater level results on the monthly scale of 2003–2016 from the existing literature [33] have been compared and verified in Figure 11, with a Pearson correlation coefficient of 0.87 and a Spearman correlation coefficient of 0.88, which has a high correlation. Thus, it can illustrate that the GWSA results of GRACE are feasible.

5. Discussion

5.1. The Spatial Trend of GWSA and Land Subsidence in BTH

When calculating the spatial trend, we used the 13-point sliding filter to deal with the seasonality. Figure 12 shows the spatial trend distribution of groundwater storage and land subsidence, and it can be seen that there is a certain degree of overlap and difference in space. In order to visually show the spatial difference, we will resample the grid data obtained from GRACE and the PS point data obtained from InSAR into a raster image of 2 km × 2 km pixels. In the study area, the largest groundwater consumption trend is 28.056 mm/year, located in the southwest of Beijing-Tianjin-Hebei. The largest land subsidence trend is 142.78 mm/year. To better understand the spatial relationship between subsidence and groundwater storage, we calculated their spatial volume. The total volume of the subsidence in the past five years is 330 million km³, and the corresponding total volume of groundwater storage consumption is 555 million km³. The total volume of land subsidence accounts for 59.46% of the total volume of groundwater consumption. Overall, from 2012 to 2016, the volume of groundwater consumption was more remarkable than the volume of land subsidence (Figure 13).

5.2. Time Series Trend of Regional GWSA and Subsidence in BTH

It shows that the whole regional average GWSA and accumulated subsidence in the whole research area reveal a declined tendency continuously from May 2012 to December 2016 in Figure 14. In addition, according to the linear fitting function, the slope of the regional average GWSA is −1.1851 and the slope of the regional average accumulated subsidence is −1.4485. Thus, the regional average GWSA and subsidence show a similar secular trend. We can obtain the annual trend of the regional average GWSA (14.221 mm) and the regional average accumulated subsidence (17.382 mm) through the slope of the linear fitting function.

5.3. Analysis of GWS and LS in Typical Areas of BTH

There are four subsidence centers in the study area with subsidence rates greater than 50 mm/year and subsidence areas greater than 80 km², including Jinzhan and Heizhuanghu in Beijing, Shengfang in Langfang and Wangqingtuo in Tianjin, and Jingxian in Hengshui and Dongguang in Cangzhou. The specific locations are shown in Figure 12 and Figure 15.

Corresponding analysis was carried out between the LS and the GWS in the four typical areas in the 0.5° × 0.5° grid calculated by GRACE. In this paper, in order to better analyze the impact of the MSWDP, the Lorenz curve model was utilized to explore the heterogeneity of LSC and GWSC from 2013 to 2016. Moreover, the Gini coefficient was used to characterize the size of the heterogeneity quantitatively.

5.3.1. Analysis of GWS and LS in Beijing

Figure 16 shows the relationship between GWSA and the monthly change of LS during 2013–2016 in the typical area of Beijing, which has a certain trend correlation. The trend of GWSA increased from 2013 to 2016, during which the trend of GWSA increased from −4.089 mm to −1.207 mm. Among these, the trend decreased from −4.089 mm in 2013 to −7.829 mm in 2014. Moreover, the trend increased from −7.829 mm in 2014 to 3.439 mm in 2015. Although the trend is slowing down, the largest seasonal consumption of GWSA is around July each year. The maximum consumption has continued to increase from 2013 to 2015 and has slowed down in 2016. Due to precipitation replenishment, the GWSA has rebounded in the second half of each year. In addition, the land subsidence showed an increasing trend from 2013 to 2014. After MSWDP in 2014, the monthly subsidence showed a significant slowing trend. It is worth mentioning that in June of each year, when the consumption of groundwater storage continues to increase, the land subsidence has a trend of slowing down. We found that it may be related to the increase in precipitation.

The Lorenz uneven curve model of GWSC and LSC in Beijing is shown in Figure 17. The heterogeneity of GWSC decreased from 0.36492 in 2013 to 0.23879 in 2014, and the heterogeneity has continued to increase since 2015. The heterogeneity of LSC in 2015 was 0.3109, which was obviously uneven compared to other years. The uneven changes in other years are more consistent. Since the implementation of the MSWDP Project in 2014, the heterogeneity of GWSC continued to increase, and the heterogeneity of land subsidence changes in the first year of South Water entering Beijing showed obvious heterogeneity and then weakened again in the following year.

5.3.2. Analysis of GWS and LS in Tianjin and Langfang

Figure 18 shows the relationship between GWSA and the monthly change of LS during 2013–2016 in the typical area of Tianjin and Langfang, which has a certain trend correlation. The trend of GWSA increased from 2013 to 2016, during which the trend of GWSA increased from −1.310 mm to 5.809 mm. The trend decreased from −1.310 mm in 2013 to 5.809 mm in 2014. Moreover, the trend increased from −5.809 mm in 2014 to 6.004 mm in 2015. Although the trend is slowing down, the enormous seasonal consumption of GWSA is around July each year. The maximum consumption has continued to increase from 2013 to 2015 and has slowed down in 2016. Due to precipitation replenishment, the GWSA has rebounded in the second half of each year. In addition, the monthly change in land subsidence showed an increasing trend from 2013 to 2014, and the monthly average subsidence increased from −7.449 mm in 2013 to −7.617 mm in 2014. However, after implementing the MSWDP Project, the monthly subsidence showed a significant slowdown, and the average monthly subsidence decreased to −6.81 mm in 2015.

The Lorenz uneven curve model of GWSC and LSC in Tianjin and Langfang is shown in Figure 19. The heterogeneity of GWSC increased from 0.30933 in 2013 to 0.53044 in 2014, and the heterogeneity has continued to decrease since 2015. The heterogeneity of LSC in 2015 was 0.44094, which was uneven compared to other years. The uneven changes in other years are more consistent. In this area, the largest heterogeneity of GWSC appeared in 2014, and the largest heterogeneity of LSC appeared in 2015. It can be seen that the heterogeneity of LSC lags behind the GWSC.

5.3.3. Analysis of GWS and LS in Hengshui and Cangzhou

Figure 20 shows the relationship between GWSA and the monthly change of LS during 2013–2016 in the typical area of Hengshui and Cangzhou. The trend of GWSA increased from 2013 to 2016, during which the trend of GWSA increased from −1.247 mm to 2.571 mm. The trend decreased from −1.247 mm in 2013 to −12.159 mm in 2014. Moreover, the trend increased from −12.159 mm in 2014 to 2.295 mm in 2015. Although the trend is slowing down, the immense seasonal consumption of GWSA is around July each year. The maximum consumption has continued to increase from 2013 to 2016. Due to rainfall replenishment, the GWSA has rebounded in the second half of each year. In addition, the monthly change in land subsidence showed an increasing trend from 2013 to 2014, and the monthly average subsidence increased from −4.868 mm in 2013 to −7.617 mm in 2014. After implementing the MSWDP Project, the monthly subsidence showed a significant slowdown, and the average monthly subsidence rose to −5.222 mm in 2015. It continued to rise to −4.261 mm in 2016. In July 2013, the groundwater consumption was severe, but the subsidence slowed down, which may be related to precipitation.

The Lorenz uneven curve model of GWSC and LSC in Hengshui and Cangzhou is shown in Figure 21. The largest heterogeneity of GWSC in 2013 was 0.45291, and the largest heterogeneity of LSC in 2014 was 0.51094, which was, obviously, uneven compared to other years. Thus, it can be seen that the heterogeneity of LSC lags behind the GWSC.

5.3.4. Analysis of GWS and LS in Baoding

Figure 22 shows the relationship between GWSA and the monthly change of LS during 2013–2016 in the typical area of Baoding. The trend of GWSA decreased from 2013 to 2016, during which the trend of GWSA decreased from −5.503 mm to −5.717 mm. The trend decreased from −5.503 mm in 2013 to −7.897 mm in 2014, and the trend increased from −7.897 mm in 2014 to 7.145 mm in 2015. Although it shows volatility in trend, the largest seasonal consumption of GWSA is around July each year. The maximum consumption has continued to increase from 2013 to 2016. Due to precipitation replenishment, the GWSA has rebounded in the second half of each year. In addition, the monthly change in land subsidence showed an increasing trend from 2013 to 2015, and the monthly average subsidence increased from −3.662 mm in 2013 to −5.768 mm in 2014. After the implementation of the MSWDP Project, the monthly subsidence showed a significant slowdown. The monthly subsidence did not slow down in 2015 but showed a clear slowing trend in 2016.

The Lorenz uneven curve model of GWSC and LSC in Baoding is shown in Figure 23. The largest heterogeneity of GWSC in 2014 was 0.48411, and the largest heterogeneity of LSC in 2014 was 0.49946, which was obviously uneven compared to other years. The heterogeneity of the GWSC presents a fluctuating change, while the heterogeneity of LSC after 2014 shows a more reduced trend.

5.4. Uncertainty Analysis

More detailed solution of GRACE mascon can be found in the literature of Scanlon et al. [58], which reduces the leakage from land to sea and does not require extra leakage calibration. We evaluated the uncertainty of TWSA by propagating the errors with combining the GRACE products provided by CSR and JPL. The slope of TWSA-CSR is −1.48, while TWSA-JPL is −1.84. The combined use of two products can better improve the effect of research.

The uncertainty of the land subsidence and groundwater storage volume during 2012–2016 calculated in this study is analyzed by the existing literatures. Zhu Juyan et al. obtained the average annual subsidence volume of the North China Plain as 1073 million m³/year until 2010 [59]. This study obtained the annual subsidence volume of the typical area in Beijing-Tianjin-Hebei from 2012 to 2016 as 660 million m³/year. Zhu Juyan et al. and Bai Lin et al., respectively, calculated the average annual consumption volume of groundwater storage in Cangzhou, finding that it was 198 million m³/year during 1970–2008 and 57–103 million m³/year during 2003–2010 [59,60]. This paper has calculated the annual consumption volume of groundwater storage from 2012 to 2016, which was 173.2 million m³/year.

6. Conclusions

Over-utilized groundwater has caused severe land subsidence in BTH. In this study, we explored the relationship between groundwater and land subsidence.

Firstly, the land subsidence results of the BTH obtained by MT-InSAR were verified with the leveling data, and the correlation coefficient is high. Then, the GRACE inversion results were verified with the in situ groundwater level from the existing literature, which maintain high timing fluctuations and correlation.

Moreover, we analyzed the relationship in time and space between groundwater storage and land subsidence in the BTH region. Within the scope of the study from 2012 to 2016, the total volume of subsidence is 330 million km³, while the total volume of groundwater storage consumption is 555 million km³. The total volume change of subsidence is 59.46% of the total volume change of groundwater storage. The regional groundwater storage is decreasing at a trend of 14.221 mm/year, and meanwhile, the regional cumulative subsidence is decreasing at a trend of 17.382 mm/year. Thus, the secular decreasing trend in the regional average GWSA and cumulative subsidence reveal a relatively high consistency.

Furthermore, we have selected four typical areas in this study. There are the Chaoyang District in Beijing; Wuqing City in Tianjin and Bazhou City in Langfang; Jing County in Hengshui and Dongguang County in Cangzhou; Gaoyang County in Baoding. The land subsidence in four typical areas and the corresponding groundwater storage in the 0.5° × 0.5° grid calculated by GRACE are analyzed. The abnormality of groundwater storage has a certain correlation and a lagging correlation with the monthly change of subsidence.

Finally, we discussed the heterogeneity of GWSC and LSC in four typical areas by the Lorenz curve model. The heterogeneity of GWSC in different regions fluctuated between 2013 and 2016, while the heterogeneity of LSC in different regions showed different patterns. It can be seen that the largest heterogeneity of LSC lags behind the GWSC in the Tianjin-Langfang-Hengshui-Baoding area. The largest uneven subsidence in Beijing and Tianjin occurred in 2015 after MSWDP Project, while the largest uneven subsidence in Hengshui-Baoding occurred in 2014. After that, the heterogeneity of subsidence gradually tends to stable.

Author Contributions

Funding acquisition and project administration, H.G. and. B.C.; methodology, W.Y. and Q.Z.; writing—original draft preparation, W.Y.; writing—review and editing, C.Z.; All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by The National Natural Science Foundation of China (No. 41930109/D010702, 41771455/D010702); Beijing Outstanding Young Scientist Program (No. BJJWZYJH01201910028032); Natural Science Foundation of Beijing Municipality (No.8182013), and Beijing Youth Top Talent Project, a program of Beijing Scholars and National “Double-Class” Construction of University Projects.

Data Availability Statement

Not applicable.

Acknowledgments

We are grateful to Center for Space Research and Jet Propulsion Laboratory for providing GRACE mascons solutions. We are grateful to the Canadian Space Agency (CSA) for their great efforts in developing and distributing the remotely sensed SAR data. We are grateful to the manufacturers of the excellent software package SarProz.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. (a) The location of the NCP in China. (b) The location of BTH in NCP. (c) The location of BTH area.

Figure 2. Schematic diagram of Lorenz curve.

Figure 3. The flowchart of this study (5.1, 5.2, and 5.3 refer to the chapters where the corresponding content appears in this paper).

Figure 4. The average land subsidence rate from 2012 to 2016.

Figure 5. The accumulated subsidence rate from 2012 to 2016.

Figure 6. The time-series subsidence changes and area of subsidence exceeding 50 mm from 2012 to 2016 in BTH.

Figure 7. Comparison between subsidence derived by InSAR and leveling results from 2015 to 2016, whose locations are shown in Figure 1.

Figure 8. Time series of total water storage anomalies (TWSA) estimated from the CSR and JPL.

Figure 9. Time series of the SMSA and SWESA from GLDAS.

Figure 10. Time series of the groundwater storage anomaly (GWSA) from GRACE and GLDAS.

Figure 11. Comparison between GWSA derived by GRACE and in situ groundwater level results from 2003 to 2016.

Figure 12. Spatial trend map of GWSA and land subsidence in BTH from 2012 to 2016. (A) The panel A is the spatial trend map of GWSA; (B) The panel B is the spatial trend map of land subsidence.

Figure 13. The relationship between groundwater storage consumption and subsidence volume in the cities of BTH.

Figure 14. Long-term time series trend of abnormal of groundwater storage and accumulated land subsidence.

Figure 15. Typical Areas in BTH (a) Chaoyang District in Beijing; (b) Wuqing City in Tianjin and Bazhou City in Langfang; (c) Jing County in Hengshui and Dongguang County in Cangzhou; (d) Gaoyang County in Baoding.

Figure 16. Chaoyang District in Beijing.

Figure 17. Lorenz curve of Chaoyang District in Beijing.

Figure 18. Wuqing City in Tianjin and Bazhou City in Langfang.

Figure 19. Lorenz curve of Wuqing City in Tianjin and Bazhou City in Langfang.

Figure 20. Jing County in Hengshui and Dongguang County in Cangzhou.

Figure 21. Lorenz curve of Jing County in Hengshui and Dongguang County in Cangzhou.

Figure 22. Gaoyang County in Baoding.

Figure 23. Lorenz curve of Gaoyang County in Baoding.

Table 1

Total water (TW) supply, groundwater (GW) supply, and MSWDP supply and precipitation (P) in BTH from 2012 to 2016.

Year	Beijing				Tianjin				Hebei
	TW Supply(km3)	GW Supply (km3)	MSWDP Supply (km3)	P(mm)	TWSupply(km3)	GWSupply (km3)	MSWDP Supply (km3)	P(mm)	TW Supply (km3)	GW Supply (km3)	MSWDP Supply (km3)	P(mm)
2012	3.59	2.04	0.28	807.63	2.01	0.55	0	820.32	19.53	15.12	0	692.81
2013	3.64	2.01	0.35	516.48	2.38	0.57	0	452.02	19.13	14.46	0	557.83
2014	3.75	1.96	0.08	427.3	2.62	0.53	0.0063	435.66	19.28	14.21	0.00012	388.2
2015	3.82	1.82	0.76	596.03	2.57	0.49	0.39	561.2	18.72	13.36	0.11	526.06
2016	3.88	1.75	0.84	661.88	2.72	0.47	0.89	623.44	18.26	12.50	0.36	599.30