Next Article in Journal
Multi-Modality and Multi-Scale Attention Fusion Network for Land Cover Classification from VHR Remote Sensing Images
Next Article in Special Issue
A Meta-Learning Approach of Optimisation for Spatial Prediction of Landslides
Previous Article in Journal
Tensor-Based Reduced-Dimension MUSIC Method for Parameter Estimation in Monostatic FDA-MIMO Radar
Previous Article in Special Issue
Risk Factor Detection and Landslide Susceptibility Mapping Using Geo-Detector and Random Forest Models: The 2018 Hokkaido Eastern Iburi Earthquake
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Key Laboratory of the Ministry of Education Land Subsidence Mechanism and Prevention, Capital Normal University, Beijing 100048, China
2
College of Resources Environment and Tourism, Capital Normal University, Beijing 100048, China
3
College of Geospatial Information Science and Technology, Capital Normal University, Beijing 100048, China
4
Observation and Research Station of Groundwater and Land Subsidence in Beijing-Tianjin-Hebei Plain, MNR, Beijing 100048, China
5
Beijing Laboratory of Water Resources Security, Capital Normal University, Beijing 100048, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(18), 3773; https://0-doi-org.brum.beds.ac.uk/10.3390/rs13183773
Submission received: 11 August 2021 / Revised: 16 September 2021 / Accepted: 17 September 2021 / Published: 20 September 2021

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 m3/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:
Φ I n S A R = Φ d e f + Φ t o p + Φ a t m + Φ o r i b i t + Φ n o i s e
Among them: Φ I n S A R is the interference phase of the point target; Φ d e f is the deformation phase of the radar line-of-sight direction; Φ t o p is the terrain phase; Φ a t m is the atmospheric delay phase; Φ o r i b i t is the orbit error phase; Φ n o i s e 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:
x ¯ = i = 1 n ( x 1 i x 2 i ) n
where n is the number of PS point in the overlapping area, and the 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].
d u = d l o s cos θ
where θ is the central incident angle of the SAR satellite, and d l o s 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.
Δ G W S = Δ T W S Δ S M S + Δ S W E S
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]:
L ( y ) = 0 y x d F ( x ) μ
where F(x) is the cumulative distribution function of ordered individuals, and μ 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:
Gini = 1 2 0 1 f ( x ) d x
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 km2. The area subsided more than 50 mm/year was expanded from 1804 km2 to 2556 km2 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 km3, and the corresponding total volume of groundwater storage consumption is 555 million km3. 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 km2, 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 m3/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 m3/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 m3/year during 1970–2008 and 57–103 million m3/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 m3/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 km3, while the total volume of groundwater storage consumption is 555 million km3. 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.

References

  1. Li, X.; Ye, S.; Song, F.; Zhou, P. Quantitative Identification of Major Factors Affecting Groundwater Change in Beijing-Tianjin-Hebei Plain. J. China Hydrol. 2018, 38, 21–27. [Google Scholar]
  2. Gong, H.; Li, X.; Pan, Y.; Zhu, L.; Zhang, Y.; Mi, C.; Chen, B.; Ke, Y.; Wang, Y.; Gao, M. Groundwater depletion and regional land subsidence of the Beijing-Tianjin-Hebei area. Bull. Natl. Nat. Sci. Found. China 2017, 31, 72–77. [Google Scholar]
  3. Cao, G.; Zheng, C.; Scanlon, B.; Jie, L.; Li, W. Use of flow modeling to assess sustainability of groundwater resources in the North China Plain. Water Resour. Res. 2013, 49, 159–175. [Google Scholar] [CrossRef]
  4. Feng, W.; Zhong, M.; Lemoine, J.; Biancale, J.; Hsu, H.; Xia, J. Evaluation of groundwater depletion in North China using the Gravity Recovery and Climate Experiment (GRACE) data and ground-based measurements. Water Resour. Res. 2013, 49, 2110–2118. [Google Scholar] [CrossRef]
  5. Lei, K.; Luo, Y.; Chen, B.; Guo, G.; Zhou, Y. Distribution characteristics and influence factors of land subsidence in Beijing area. Geol. China 2016, 43, 2216–2228. [Google Scholar] [CrossRef]
  6. Zhang, Y.; Hongan, W.; Kang, Y. Ground Subsidence over Beijing-Tianjin-Hebei Region during Three Periods of 1992 to 2014 Monitored by Interferometric SAR. Acta Geod. Et Cartogr. Sin. 2016, 45, 1050–1058. [Google Scholar] [CrossRef]
  7. Zhang, J.; Liu, J.; Zhai, L.; Hou, W. Implementation of Geographical Conditions Monitoring in Beijing-Tianjin-Hebei, China. Int. J. Geo-Inf. 2016, 5, 89. [Google Scholar] [CrossRef] [Green Version]
  8. Liu, G. Application examples of InSAR and its limitation analysis. Surv. Mapp. Sichuan 2005, 3, 44–48. [Google Scholar] [CrossRef]
  9. Gabriel, A.; Goldstein, R.; Zebker, H. Mapping Small Elevation Changes Over Large Areas: Differential Radar Interferometry. J. Geophys. Res. Solid Earth 1989, 94, 9183–9184. [Google Scholar] [CrossRef]
  10. Hanssen, R.F. Radar Interferometry; Springer: Berlin/Heidelberg, Germany, 2001. [Google Scholar]
  11. Ghiglia, D.C.; Pritt, M.D. Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software. 1998. Available online: https://xueshu.baidu.com/usercenter/paper/show?paperid=f7a53cd9173485b92a3f0756cf0c1d84 (accessed on 10 August 2021).
  12. Zebker, H.; Rosen, P.; Hensley, S. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. J. Geophys. Res. Solid Earth 1997, 102, 7547–7563. [Google Scholar] [CrossRef]
  13. Xu, C.; He, P.; Wen, Y.; Liu, Y. Recent advances InSAR interferometry and its applications. J. Geomat. 2015, 40, 1–9. [Google Scholar] [CrossRef]
  14. China Geological Environmental Monitoring Institute. China Geological Environmental Monitoring Groundwater Level Yearbook; China Dadi Publishing House: Beijing, China, 2007. [Google Scholar]
  15. Li-Tang, H.; Sun, K.; Yin, W. Review on the Application of GRACE Satellite in Regional Groundwater Management. J. Earth Sci. Environ. 2016, 38, 258–266. [Google Scholar]
  16. Strassberg, G.; Scanlon, B.; Rodell, M. Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer (USA). Geophys. Res. Lett. 2007, 34. [Google Scholar] [CrossRef] [Green Version]
  17. Luo, Z.; Li, Q.; Zhong, P. Water storage variations in Heihe river recovered from GRACE temporal gravity field. Acta Geod. Cartogr. Sin. 2012, 41, 676–681. [Google Scholar]
  18. Xie, X.; Caijun, X.; Gong, Z.; Wei, L. Groundwater Storage Changes in Shan-Gan-Jin Plateau Derived from GRACE. Bull. Surv. Mapp. 2018, 1, 133–137. [Google Scholar]
  19. Srivastava, S.; Dikshit, O. Seasonal and trend analysis of TWS for the Indo-Gangetic plain using GRACE data. Geocarto Int. 2019, 35, 1343–1359. [Google Scholar] [CrossRef]
  20. Singh, L.; Subbarayan, S. Satellite-derived GRACE groundwater storage variation in complex aquifer system in India. Sustain. Water Resour. Manag. 2020, 6, 43. [Google Scholar] [CrossRef]
  21. Rodell, M.; Chen, J.; Kato, H.; Famigliettim, J.; Nigro, J.; Wilson, C. Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE. Hydrogeol. J. 2007, 15, 159–166. [Google Scholar] [CrossRef] [Green Version]
  22. Famiglietti, J.; Lo, M.; Ho, S.; Bethune, J.; Anderson, K.; Syed, T.; Swenson, S.; Linage, C.; Rodell, M. Satellites measure recent rates of groundwater depletion in California’s Central Valley. Geophys. Res. Lett. 2011, 38. [Google Scholar] [CrossRef] [Green Version]
  23. Voss, K.; Famiglietti, J.; Lo, M.; Linage, C.; Rodell, M.; Swenson, S. Groundwater depletion in the Middle East from GRACE with implications for transboundary water management in the Tigris-Euphrates-Western Iran region. Water Resour. Res. 2013, 49, 904–914. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Joodaki, G.; Wahr, J.; Swenson, S. Estimating the human contribution to groundwater depletion in the Middle East, from GRACE data, land surface models, and well observations. Water Resour. Res. 2014, 50, 2679–2692. [Google Scholar] [CrossRef]
  25. Velicogna, I.; Tong, J.; Zhang, T.; Kimball, J. Increasing subsurface water storage in discontinuous permafrost areas of the Lena River basin, Eurasia, detected from GRACE. Geophys. Res. Lett. 2012, 39, 1–5. [Google Scholar] [CrossRef] [Green Version]
  26. Suzuki, K.; Matsuo, K.; Yamazaki, D.; Ichii, K.; Lijima, Y.; Paga, F.; Yanagi, Y.; Hiyama, T. Hydrological variability and changes in the Arctic circumpolar tundra and the three largest pan-Arctic river basins from 2002 to 2016. Remote Sens. 2018, 10, 402. [Google Scholar] [CrossRef] [Green Version]
  27. Pan, Y.; Zhang, C.; Gong, H.; Yeh, P.; Shen, Y.; Guo, Y.; Huang, Z.; Li, X. Detection of human--induced evapotranspiration using GRACE satellite observations in the Haihe River basin of China. EGU Gen. Assem. Conf. Abstr. 2017, 44, 190–199. [Google Scholar] [CrossRef]
  28. Gao, M.; Gong, H.; Chen, B.; Zhou, C.; Chen, W. InSAR time-series investigation of long-term ground displacement at Beijing Capital International Airport, China. Tectonophysics 2016, 691, 271–281. [Google Scholar] [CrossRef]
  29. Chen, B.; Gong, H.; Li, X.; Lei, K.; Zhu, L.; Gao, M.; Zhou, C. Characterization and causes of land subsidence in Beijing, China. Int. J. Remote. Sens. 2017, 38, 808–826. [Google Scholar] [CrossRef]
  30. Castellazzi, P.; Martel, R.; Rivera, A.; Huang, J.; Pavlic, G.; Calderhead, A.; Chaussard, E.; Garfias, J.; Salas, J. Groundwater depletion in Central Mexico: Use of GRACE and InSAR to support water resources management. Water Resour. Res. 2016, 52, 5985–6003. [Google Scholar] [CrossRef]
  31. Du, Z.; Ge, L.; Ng, H.; Li, X. Time series interferometry integrated with groundwater depletion measurement from grace. In Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; IEEE: Manhattan, NY, USA, 2016. [Google Scholar] [CrossRef]
  32. Huang, Z.; Pan, Y.; Gong, H.; Yeh, J.; Li, X.; Zhou, D.; Zhao, W. Subregional-scale groundwater depletion detected by GRACE for both shallow and deep aquifers in North China Plain. Geophys. Res. Lett. 2015, 42, 1791–1799. [Google Scholar] [CrossRef]
  33. Gong, H.; Pan, Y.; Zheng, L.; Li, X.; Zhu, L.; Zhang, C.; Huang, Z.; Li, Z.; Wang, H.; Zhou, C. Long-term groundwater storage changes and land subsidence development in the North China Plain (1971–2015). Hydrogeol. J. 2018, 26, 1417–1427. [Google Scholar] [CrossRef] [Green Version]
  34. Zhang, Z.; Shen, Z.; Xue, Y. Evolution of Groundwater Environment in North China Plain; Geological Publishing House: Bath, UK, 2000. [Google Scholar]
  35. Fei, Y.; Zhang, Z.; Zhang, F.; Chen, J.; Chen, Z.; Wang, Z. Factors affecting dynamic variation of groundwater level in North China Plain. J. Hehai Univ. 2005, 5, 538–541. [Google Scholar] [CrossRef]
  36. Chen, Z.; Qi, J.; Xu, J.; Xu, J.; Hao, Y.; Nan, Y. Paleoclimatic interpretation of the past 30 ka from isotopic studies of the deep confined aquifer of the North China Plain. Appl. Geochem. 2003, 18, 997–1009. [Google Scholar]
  37. Guo, H.; Li, W.; Wang, L.; Chen, Y.; Zang, X.; Wang, Y.; Zhu, J.; Bian, Y. Present situation and research prospects of the land subsidence driven by groundwater levels in the North China Plain. Hydrogeol. Eng. Geol. 2021, 48, 1–13. [Google Scholar]
  38. Cui, W.; Lei, K. Some Ideas on Land Subsidence Working from the view of Coordinated Development in Beijing-Tianjin-Hebei Regions. Urban Geol. 2018, 13, 1007–1903. [Google Scholar]
  39. Yi, L.; Fang, Z.; He, X.; Chen, S.; Wei, W.; Qiang, Y. Land subsidence in Tianjin, China. Environ. Earth Sci. 2011, 62, 1151–1161. [Google Scholar]
  40. Xue, X.; Li, J.; Xie, X.; Qian, K.; Wang, Y. Impacts of sediment compaction on iodine enrichment in deep aquifers of the North China Plain. Water Res. 2019, 159, 480–489. [Google Scholar] [CrossRef] [PubMed]
  41. Guo, H.; Bai, J.; Zhang, Y.; Wang, L.; Shi, J.; Li, W.; Zhang, Z.; Wang, Y.; Zhu, J.; Wang, H. The evolution characteristics and mechanism of the land subsidence in typical areas of the North China Plain. Geol. China 2017, 44, 1115–1127. [Google Scholar]
  42. Beijing Institute of Hydrogeology and Engineering Geology. The Annual Report of Beijing Land Subsidence Monitoring; Beijing Institute of Hydrogeology and Engineering Geology: Beijing, China, 2015. [Google Scholar]
  43. Yi, C. The latest progress of land subsidence control in Tianjin. Haihe Water Resour. 2017, 64, 42–43. [Google Scholar]
  44. Liu, X.; Wang, Y.; Yan, S. Ground deformation associated with exploitation of deep groundwater in Cangzhou City measured by multi-sensor synthetic aperture radar images. Environ. Earth. 2017, 76, 6. [Google Scholar] [CrossRef]
  45. Ferretti, A.; Prati, C. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Trans. Geosci. Remote. Sens. 2000, 38, 2202–2212. [Google Scholar] [CrossRef] [Green Version]
  46. Zhang, Y.; Wu, H.; Kang, Y.; Zhu, C. Ground Subsidence in the Beijing-Tianjin-Hebei Region from 1992 to 2014 Revealed by Multiple SAR Stacks. Remote. Sens. 2016, 8, 675. [Google Scholar] [CrossRef] [Green Version]
  47. Tapley, B.; Bettadpur, S.; Ries, J.; Thompson, P.; Watkins, M. GRACE Measurements of Mass Variability in the Earth System. Science 2004, 305, 503–505. [Google Scholar] [CrossRef] [Green Version]
  48. Yeh, J.; Swenson, S.; Famiglietti, J.; Rodell, M. Remote Sensing of Groundwater Storage Changes in Illinois Using the Gravity Recovery and Climate Experiment (GRACE). Water Resour. Res. 2006, 42, 12. [Google Scholar] [CrossRef]
  49. Massoud, E.; Purdy, A.; Miro, M.; Famiglietti, J. Projecting Groundwater Storage Changes in California’s Central Valley. Entific Rep. 2018, 8, 12917. [Google Scholar] [CrossRef]
  50. Zhang, C.; Duan, Q.; Yeh, P.; Pan, Y. Sub-regional groundwater storage recovery in NCP after the South-to-North water diversion project. J. Hydrol. 2021, 597, 126–156. [Google Scholar] [CrossRef]
  51. Wang, J.; Song, C.; Reager, J.; Yao, F.; Famiglietti, J.; Sheng, Y.; MacDonald, G.; Brun, F.; Schmied, H.; Marston, R.; et al. Recent global decline in endorheic basin water storages. Nat. Geosci. 2018, 11, 926–932. [Google Scholar] [CrossRef] [Green Version]
  52. Lorenz, M.O. Methods for measuring the concentration of wealth. Am. Stat. Assoc. 1905, 70, 209–219. [Google Scholar] [CrossRef]
  53. Sangüesa, C.; Pizarro, R.; Ibañez, A.; Ingram, B.; Rivera, D.; Garcia, C.; Ingram, B. Spatial and Temporal Analysis of Rainfall Concentration Using the Gini Index and PCI. Water 2018, 10, 112. [Google Scholar] [CrossRef] [Green Version]
  54. Zou, S.; Abuduwaili, J.; Duan, W.; Ding, J.; Ma, L. Attribution of changes in the trend and temporal non-uniformity of extreme precipitation events in Central Asia. Sci. Rep. 2021, 11, 1–11. [Google Scholar] [CrossRef] [PubMed]
  55. Wanling, X.; Zhu, W.; Zhang, J.; Zheng, X.; Jin, H. Analysis on Temporal Inhomogenenity of Runoff in Tumen River Mainstream Based on Lorenz Curve. Bull. Soil Water Conserv. 2015, 35, 128–132. [Google Scholar]
  56. Save, H.; Bettadpur, S.; Tapley, B. High-resolution CSR GRACE RL05 mascons. J. Geophys. Res. Solid Earth 2016, 121, 7547–7569. [Google Scholar] [CrossRef]
  57. Wiese, D.N.; Landerer, F.W.; Watkins, M.M. Quantifying and reducing leakage errors in the JPL RL05M GRACE mascon solution. Water Resour. Res. 2016, 52, 7490–7502. [Google Scholar] [CrossRef]
  58. Scanlon, B.; Zhang, Z.; Save, H.; Wiese, D.; Landerer, F.; Long, D.; Longuevergne, L.; Chen, J. Global Evaluation of New GRACE Mascon Products for Hydrologic Applications. Water Resour. Res. 2016, 52, 9412–9429. [Google Scholar] [CrossRef]
  59. Zhu, J.; Guo, H.; Peng, L.; Tian, X. Relationship between Land Subsidence and Deep Groundwater Yield in the North China Plain. South--North Water Transf. Water Sci. Technol. 2014, 12, 165–169. [Google Scholar]
  60. Jiang, L.; Bai, L.; Zhao, Y.; Gao, G.; Wang, H.; Sun, Q. Combining InSAR and Hydraulic Head Measurements to Estimate Aquifer Parameters and Storage Variations of Confined Aquifer System in Cangzhou, North China Plain. Water Resour. Res. 2018, 54, 8234–8252. [Google Scholar] [CrossRef]
Figure 1. (a) The location of the NCP in China. (b) The location of BTH in NCP. (c) The location of BTH area.
Figure 1. (a) The location of the NCP in China. (b) The location of BTH in NCP. (c) The location of BTH area.
Remotesensing 13 03773 g001
Figure 2. Schematic diagram of Lorenz curve.
Figure 2. Schematic diagram of Lorenz curve.
Remotesensing 13 03773 g002
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 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).
Remotesensing 13 03773 g003
Figure 4. The average land subsidence rate from 2012 to 2016.
Figure 4. The average land subsidence rate from 2012 to 2016.
Remotesensing 13 03773 g004
Figure 5. The accumulated subsidence rate from 2012 to 2016.
Figure 5. The accumulated subsidence rate from 2012 to 2016.
Remotesensing 13 03773 g005
Figure 6. The time-series subsidence changes and area of subsidence exceeding 50 mm from 2012 to 2016 in BTH.
Figure 6. The time-series subsidence changes and area of subsidence exceeding 50 mm from 2012 to 2016 in BTH.
Remotesensing 13 03773 g006
Figure 7. Comparison between subsidence derived by InSAR and leveling results from 2015 to 2016, whose locations are shown in Figure 1.
Figure 7. Comparison between subsidence derived by InSAR and leveling results from 2015 to 2016, whose locations are shown in Figure 1.
Remotesensing 13 03773 g007
Figure 8. Time series of total water storage anomalies (TWSA) estimated from the CSR and JPL.
Figure 8. Time series of total water storage anomalies (TWSA) estimated from the CSR and JPL.
Remotesensing 13 03773 g008
Figure 9. Time series of the SMSA and SWESA from GLDAS.
Figure 9. Time series of the SMSA and SWESA from GLDAS.
Remotesensing 13 03773 g009
Figure 10. Time series of the groundwater storage anomaly (GWSA) from GRACE and GLDAS.
Figure 10. Time series of the groundwater storage anomaly (GWSA) from GRACE and GLDAS.
Remotesensing 13 03773 g010
Figure 11. Comparison between GWSA derived by GRACE and in situ groundwater level results from 2003 to 2016.
Figure 11. Comparison between GWSA derived by GRACE and in situ groundwater level results from 2003 to 2016.
Remotesensing 13 03773 g011
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 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.
Remotesensing 13 03773 g012
Figure 13. The relationship between groundwater storage consumption and subsidence volume in the cities of BTH.
Figure 13. The relationship between groundwater storage consumption and subsidence volume in the cities of BTH.
Remotesensing 13 03773 g013
Figure 14. Long-term time series trend of abnormal of groundwater storage and accumulated land subsidence.
Figure 14. Long-term time series trend of abnormal of groundwater storage and accumulated land subsidence.
Remotesensing 13 03773 g014
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 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.
Remotesensing 13 03773 g015
Figure 16. Chaoyang District in Beijing.
Figure 16. Chaoyang District in Beijing.
Remotesensing 13 03773 g016
Figure 17. Lorenz curve of Chaoyang District in Beijing.
Figure 17. Lorenz curve of Chaoyang District in Beijing.
Remotesensing 13 03773 g017
Figure 18. Wuqing City in Tianjin and Bazhou City in Langfang.
Figure 18. Wuqing City in Tianjin and Bazhou City in Langfang.
Remotesensing 13 03773 g018
Figure 19. Lorenz curve of Wuqing City in Tianjin and Bazhou City in Langfang.
Figure 19. Lorenz curve of Wuqing City in Tianjin and Bazhou City in Langfang.
Remotesensing 13 03773 g019
Figure 20. Jing County in Hengshui and Dongguang County in Cangzhou.
Figure 20. Jing County in Hengshui and Dongguang County in Cangzhou.
Remotesensing 13 03773 g020
Figure 21. Lorenz curve of Jing County in Hengshui and Dongguang County in Cangzhou.
Figure 21. Lorenz curve of Jing County in Hengshui and Dongguang County in Cangzhou.
Remotesensing 13 03773 g021
Figure 22. Gaoyang County in Baoding.
Figure 22. Gaoyang County in Baoding.
Remotesensing 13 03773 g022
Figure 23. Lorenz curve of Gaoyang County in Baoding.
Figure 23. Lorenz curve of Gaoyang County in Baoding.
Remotesensing 13 03773 g023
Table 1. Total water (TW) supply, groundwater (GW) supply, and MSWDP supply and precipitation (P) in BTH from 2012 to 2016.
Table 1. Total water (TW) supply, groundwater (GW) supply, and MSWDP supply and precipitation (P) in BTH from 2012 to 2016.
YearBeijingTianjinHebei
TW Supply
(km3)
GW Supply (km3)MSWDP Supply (km3)P
(mm)
TW
Supply
(km3)
GW
Supply (km3)
MSWDP Supply (km3)P
(mm)
TW Supply (km3)GW Supply (km3)MSWDP Supply (km3)P
(mm)
20123.592.040.28807.632.010.550820.3219.5315.120692.81
20133.642.010.35516.482.380.570452.0219.1314.460557.83
20143.751.960.08427.32.620.530.0063435.6619.2814.210.00012388.2
20153.821.820.76596.032.570.490.39561.218.7213.360.11526.06
20163.881.750.84661.882.720.470.89623.4418.2612.500.36599.30
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yu, W.; Gong, H.; Chen, B.; Zhou, C.; Zhang, Q. Combined GRACE and MT-InSAR to Assess the Relationship between Groundwater Storage Change and Land Subsidence in the Beijing-Tianjin-Hebei Region. Remote Sens. 2021, 13, 3773. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13183773

AMA Style

Yu W, Gong H, Chen B, Zhou C, Zhang Q. Combined GRACE and MT-InSAR to Assess the Relationship between Groundwater Storage Change and Land Subsidence in the Beijing-Tianjin-Hebei Region. Remote Sensing. 2021; 13(18):3773. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13183773

Chicago/Turabian Style

Yu, Wen, Huili Gong, Beibei Chen, Chaofan Zhou, and Qingquan Zhang. 2021. "Combined GRACE and MT-InSAR to Assess the Relationship between Groundwater Storage Change and Land Subsidence in the Beijing-Tianjin-Hebei Region" Remote Sensing 13, no. 18: 3773. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13183773

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop