Introduction

Groundwater sustains social and environmental development in many arid and semi-arid regions; however, over-abstraction of groundwater, induced by drought and large water demand for agricultural and domestic use, has become a global issue (Famiglietti 2014; Richey et al. 2015). The Gravity Recovery and Climate Experiment (GRACE) satellites show that many aquifers around the world are being depleted at amazing rates such as in northwest India (Long et al. 2016; Rodell et al. 2009), Central Valley in California, USA (Famiglietti et al. 2011), and the North China Plain (NCP; Feng et al. 2013; Huang et al. 2015; Pan et al. 2017). More information can be found from the review work by Chen et al. (2016a and Wada (2016).

If only focusing on the groundwater depletion rate, the NCP (Fig. 1) does not seem to be the worst among the well-known aquifers worldwide (Chen et al. 2016a; Famiglietti et al. 2011; Rodell et al. 2009); however, groundwater over-abstraction accompanied by the compression of aquifer systems has led to serious land subsidence in NCP (Guo et al. 2015; Ye et al. 2016). The subsidence-affected area in the NCP is around 120,000 km2 (86% of the total area), and an area of 17,500 km2 has been subsiding at a rate of ≥20 mm/yr (Ye et al. 2016). The Interferometric Synthetic Aperture Radar (InSAR) measured a subsidence rate up to 15 cm/yr in Beijing during 2003–2011 (Chen et al.2016b; Zhang et al. 2016). The accumulated displacement of the city of Cangzhou, Hebei province, one of the most serious land subsidence sites in NCP, is up to around 3 m. Land subsidence may cause increased flood risk, structural damage, and seawater intrusion, and result in significant socio-economic distress (Hwang et al. 2016). The highly developed megacities (e.g., Beijing) with lots of large constructions, high-speed railways, and dense population, are under the potential risk of land surface subsiding.

Fig. 1
figure 1

Location map of the study area and the in situ observation gauges, and spatial map of specific yield according to Zhang et al. (2009)

Full size image

One may expect amelioration regarding the situation of groundwater over-abstraction and land subsidence after the implementation of the South-to-North Water Diversion Project in December 2014, which was designed to transfer 9.6 km3 of water per year from the south (Danjiangkou reservoir) to the north (including Beijing, Tianjin, Hebei and Henan). The government plans to reduce groundwater use in Beijing, Tianjin and Hebei province to 65% of current levels by 2023. Such a dramatic change to groundwater consumption may recover the groundwater systems in NCP to some extent.

Numerous studies have explored the groundwater storage (GWS) depletion and land subsidence in NCP; however, previous studies usually focused on a limited time span in this area. For instance, Huang et al. (2015) and Pan et al. (2017) quantified and analyzed the GWS depletion during 2003–2014, while a few other researchers studied the land subsidence during 1992–2014 (e.g. Chen et al. 2016b, 2017; Zhang et al. 2016; Zhu et al. 2015). In order to provide more information for future groundwater resources management, more studies are still needed to investigate the GWS dynamics and the evolution of the aquifer system from a long-term perspective (Castellazzi et al. 2016a, b).

This study aims to explore the long-term changes (1971–2015) of GWS and land subsidence in the NCP, as well as their relationship using in situ observations and satellite measurements from GRACE and InSAR data collected from literature. The spatial and temporal variations of land subsidence will be discussed to provide a historical overview of the groundwater system in NCP and its three sub-regions, i.e., Beijing, Tianjin and Hebei provinces.

Data and methods

In situ groundwater level observations

Historical monthly in situ groundwater level data (Fig. 1) from 1971 to 2013 were collected, including the records of 44 wells in NCP during 1971–1983, obtained from Water Resources Yearbook of Haihe River Basin compiled by the Ministry of Water Resources of China; 5 wells in Beijing during 1971–2013 provided by Beijing Institute of Hydrogeology and Engineering Geology; and 81 wells in NCP during 2005–2013 from the Groundwater Level Almanac compiled by the China Institute of Geological Environment Monitoring. There was a data gap for the entire NCP during 1984–2002. In this study, the continuous groundwater level observations of the five wells in Beijing, accompanied with reported data about the groundwater level changing rate in NCP, were interpolated to obtain a regional averaged estimate of monthly groundwater storage anomaly (GWSA) during 1984–2002. The observed groundwater level changes were converted into GWS changes by multiplying the specific yield according to Zhang et al. (2009; Fig. 1b) as used by Huang et al. (2015). A 20% uncertainty is assumed for the specific yield. The monthly GWSA was then obtained by subtracting mean GWS changes of the whole period (1971–2015).

GRACE-derived groundwater storage anomalies

Monthly GRACE level-3 Release-05 mascon (i.e., mass concentration) products (2003–2015) from the Center for Space Research (CSR) at the University of Texas and Jet Propulsion Laboratory (JPL), USA are utilized. GRACE mascon solutions represent total terrestrial water storage (TWS) anomalies. In arid region like NCP, the TWS values are the aggregation of surface-water storage in lakes and reservoirs, soil moisture and groundwater. In order to estimate GRACE-derived GWS anomalies, the observed reservoir water storage (from 22 primary reservoirs) and the averaged soil moisture storage (SMS) simulated by four land surface models (Noah, Mosaic, Community Land Model and Variable Infiltration Capacity Model) from the Global Land Data Assimilation System (GLDAS; Rodell et al. 2004) were subtracted from the mascon-based TWS anomalies (TWSA). The modeled SMS data from GLDAS are validated against in situ observations from 97 monitoring sites (Fig. 1) provided by China Meteorological Administration.

The GRACE mascon solution has reduced leakage from land to ocean, and no additional leakage correction is required (Scanlon et al. 2016). More details about the mascon solution can be found in Save et al. (2016) and Watkins et al. (2015). The missing values in GRACE data are conservatively estimated through a cubic-spline interpolation. The uncertainties of GWSA are estimated using the propagated errors from GRACE TWSA provided by mascon solutions and the modeled SMS estimated using the standard deviation of the four model simulations.

Land subsidence

The regional averaged annual land subsidence rates in Beijing eastern plain (1971–2014), Cangzhou (1971–2004, 2007, 2009–2010, 2013–2014), and Tianjin urban area (1971–2013) were obtained from previous studies (Wang 2004; Liu et al. 2016; Chen et al. 2010; Pan et al. 2004; Guo 2016). These data were obtained through ground leveling conducted by the institutes of local geologic survey; additionally, the land subsidence in Beijing plain area was estimated from 29 RADARSAT-2 InSAR images from 2011 to 2015. They were processed using the small baseline approaches following Zhou et al. (2016).

Auxiliary data

Groundwater abstraction and precipitation data were used for correlation analysis with the GRACE-derived GWSA. Annual groundwater abstraction data (1971–2015) were provided by Beijing Institute of Hydrogeology and Engineering Geology and collected from Water Resources Yearbook of Haihe River Basin. Monthly gridded (0.5° × 0.5°) precipitation data (version 2.0) are obtained from the China Metrological Administration (CMA). This dataset has been produced through a thin plate spline interpolation using the data from 2,472 rain gauge stations across China. It has been well validated over the study area with 1,684 rain gauge records (2010–2012) in a previous study by Pan et al. (2017).

Results and discussion

Long-term GWS depletion exacerbated by droughts

Figure 2 plots the monthly GWSA of NCP during the period 1971–2015. It is apparent that the GWS in NCP showed prolonged declining trend since the 1970s with a rate of −17.8 ± 0.1 mm/yr (Fig. 2a; Table 1). However, variabilities of GWS depletion associated with the changes of precipitation and abstraction can be found for the sub-periods (Fig. 2b,c). The GWSA slightly decreased (−6.2 ± 1.3 mm/yr) during 1971–1980, followed by a rapid declining (−20.6 ± 0.2 mm/yr) during 1981–2002 with less precipitation and more abstraction. In recent years (2003–2015), the GWSA represented a similar declining trend (−18.2 ± 0.2 mm/yr) with slightly increased mean precipitation and abstraction as compared with that during 1981–2002, despite the fact that the abstraction gradually reduced year by year. The long-term declining trend is a combined result from natural and anthropogenic factors, as drought usually lead to more intensive groundwater abstraction due to increased agricultural water demand (Pan et al. 2017). For continuous dry years (e.g. 1997–2002), the groundwater depletion rate can be as large as −34.7 ± 0.9 mm/yr.

Fig. 2
figure 2

a Monthly time series of the GWSA in NCP from 1971 to 2015. The blue curves are the GWSA estimated from in situ groundwater-level observations, while the red curves are that from GRACE data. The subplot shows the blown-up GWSA time series from GRACE and in situ observations during 2005–2013. The shaded areas are the monthly GWSA errors. Annual rainfall (b) and annual groundwater abstraction data (c) are shown for reference. The straight lines are the linearly fitted trends (a) and long-term mean annual rainfall (b) and groundwater abstraction (c) during the corresponding sub-period

Full size image
Table 1 Fitted seasonal terms (annual and semi-annual amplitudes and phases) and long-term trends of the monthly GWSA time series (Fig. 2) derived from GRACE and in situ observations based on the least-squares harmonic analysis (Jin and Feng 2013; Yan et al. 2006)
Full size table

During 1971–2015, Beijing experienced the largest GWS depletion rate (−34.1 ± 2.4 mm/yr) as compared with Tianjin (−7.0 ± 0.4 mm/yr) and Hebei (−15.2 ± 0.6 mm/yr; Fig. 3). An aggravated depletion (−76.1 ± 6.5 mm/yr) occurred in Beijing since 1999 due to the continuous drought (Chen et al. 2017). In Hebei, the GWS depletion in recent years (2005–2013) is only half of the historical record (−14.7 ± 1.1 mm/yr) of 1971–1990. The GWS in Tianjin experienced accelerated depletion in recent years (2005–2013) at the rate of −20.2 ± 2.0 mm/yr.

Fig. 3
figure 3

The 1971–2013 yearly time series of precipitation, groundwater abstraction, land subsidence, and GWSA in three sub-regions of a Beijing, b Tianjin, and c Hebei. Land subsidence data of typical regions of Beijing eastern plain, Tianjin urban area and Cangzhou was used for Beijing, Tianjin and Hebei, respectively. There are missing records in parts of the years

Full size image

Figure 4 presents temporal relationships between the GWSA, precipitation, and groundwater abstraction for three sub-regions, i.e. Beijing, Tianjin and Hebei. Poor relationships exist between GWSA and precipitation (Fig. 4a–c), while significant correlations can be found between the GWSA and abstraction (Fig. 4d–f). In the early depletion stage (GWSA >0), negative correlations indicate that increased abstraction can significantly deplete the aquifer. In contrast, positive relationships exist in the late depletion stage (GWSA <0), indicating that the declining trend cannot be quickly recovered through reducing abstraction. This is because the reduced abstraction is still unsustainable and it cannot be totally offset by groundwater recharge (Zhang et al. 2009). The response of the GWSA to the changes of abstraction varies among the sub-regions. The GWSA of Beijing is more sensitive to precipitation and abstraction as compared to Tianjin and Hebei (Fig. 4g–i). It means that the GWS of Beijing can be significantly recovered by intensive precipitation (e.g. 1985–1997), or be depleted by continuous droughts (e.g. 1999–2012).

Fig. 4
figure 4

Statistical relationships between the GWSA and a–c precipitation, d–f annual groundwater abstraction, and g–i cumulative groundwater abstraction for Beijing (left panel), Tianjin (middle panel) and Hebei (right panel). The shaded area indicates the 95% confidence interval

Full size image

Table 1 illustrates the least-squares-fitted seasonal terms (i.e., annual and semi-annual amplitudes and phases) and long-term trends of GWSA. In situ GW-level observations from 2005 to 2013 are used in comparison with GRACE-derived GWSA. The two GWSA estimates show a close depletion rate of −17.7 ± 1.1 mm/yr (GRACE) and − 19.2 ± 1.7 mm/yr (in situ) during 2005–2013, with a correlation coefficient of 0.84 and root mean square error (RMSE) of 37.6 mm. Regarding the larger amplitude represented by water wells comparing to that from GRACE (Fig. 2; Table 1), the reasons are probably due to errors associated with specified yield and the GRACE data (Scanlon et al. 2012). Joint use of the results from GRACE and in situ observations lead to caution about the uncertainties.

Sub-regional land subsidence responding to natural and anthropogenic factors

Continuous GWS depletion has led to a large area of land subsidence in the NCP, with the accumulated deformation up to ~3 m in Cangzhou (Fig. 5). The changes of the GWS and land subsidence of sub-periods are given in Table 2. Large subsidence rates (>50 mm/yr) mainly occurred in the late (2003–2015), middle (1980–2000), and early (1971–1987) periods in Beijing, Hebei, and Tianjin, respectively. The land subsidence generally responds to the changes of groundwater abstraction and GWSA with a good correlation, while it shows no significant correlations with annual precipitation (Fig. 6). During the early stage (before 1980s), positive correlations demonstrate that the subsidence was mainly regulated by abstraction and, thus, increased pumping would lead to worse subsidence. Negative correlations existing in the late stage (after the 1980s) indicate that the land subsidence continued with reduced regional abstraction during this period. This has probably resulted from the changed hydrogeological conditions due to the drawdown of groundwater level. Budhu and Adiyaman (2013) found that subsidence continued when pumping is stopped if the groundwater level is below the elevation of the clay zones. In contrast, the GWSA represents consistent negative correlations with subsidence during the whole period, except for Hebei in the late stage because subsidence data of the city of Cangzhou used for analysis fails to represent regional land subsidence of Hebei province during this period.

Fig. 5
figure 5

The graduated color map delineates the spatial distribution of GRACE-derived GWSA trend (2003–2015) and the black contour lines show the accumulated land subsidence from 1971 to 2003 in NCP

Full size image
Table 2 Statistics of annual precipitation (P), groundwater (GW) abstraction, mean subsidence rate, and GWSA trend in the three sub-regions of Beijing, Tianjin, and Hebei during 1971–2013
Full size table
Fig. 6
figure 6

Statistical relationships between annual subsidence and a–c precipitation, d–f groundwater abstraction, and g–i GWSA, for Beijing (left panel), Tianjin (middle panel) and Hebei (right panel). The shaded area indicates the 95% confidence interval

Full size image

The different correlations represented by abstraction and GWSA, as well as their variations in space and time, resulted from the combined effects of natural and anthropogenic factors. Taking Beijing as an example, GWS depletion due to increasing abstraction caused a notable development of subsidence during 1970s, followed by stable GWS changes and subsidence during 1980–1998; however, the 9-year continuous drought that started from 1999 again caused a rapid increase of subsidence, despite reduced abstraction. In this case, the abstraction failed to represent subsidence development with a positive correlation shown in the early stage. In contrast, the GWSA well reflected the variations of subsidence with negative correlations during the whole period as it integrated the information of both discharge and recharge.

It should be noted that the lack of long-term land subsidence data at large scale limits the understanding of the spatio-temporal variabilities of the aquifer system. The relationships derived from this study may vary as the scale changes; besides, the hysteresis of subsidence and lithological differences may also contribute to the spatio-temporal variations (Guo et al. 2015).

Seasonal variations revealed by satellite data

Figure 7a plots the spatial pattern of mean land subsidence rate (2011–2015) from RADARSAT-2 InSAR data and Fig. 7b compares the mean annual cycle (2011–2015) of GWSA from in situ groundwater-level observations and InSAR-derived land deformation in the Beijing plain area. The RADARSAT-2 InSAR data are of high spatial resolution for explaining the distribution and driving factors of land subsidence (Zhou et al. 2016) and provide much temporal information to reveal its seasonal variations. As shown in Fig 7b, the land surface deformation represents a similar seasonal cycle comparing to GWSA. Both GWS depletion and land surface deformation evidence accelerating rates from March to June, and recover in the second half year. It is believed that this seasonal cycle has resulted from intensive groundwater pumping during March–May when little precipitation is available for wheat growth, while precipitation during the rainy season (June–September) sustains the water demand for maize (Pan et al. 2017).

Fig. 7
figure 7

a Spatial distribution of the mean subsidence rate in Beijing for the period 2011–2015. b Seasonal cycle of in situ GWSA and InSAR-derived land deformation during 2011 and 2015. The diagonal-patterned area (a) indicates the cone of depression, and its seasonal deformation is illustrated with the green-diamond curve (b)

Full size image

Conclusions

This paper provides an overview of the long-term evolution of the GWS and land subsidence in the NCP, with a focus on the spatio-temporal variation of the aquifer systems responding to natural and anthropogenic factors. The GWS shows a prolonged declining rate of −17.8 ± 0.1 mm/yr during the whole study period 1971–2015, while more intensive depletion (−76.1 ± 6.5 mm/yr) can be found in Beijing during 1999–2012 due to a precipitation deficit. Large subsidence rates (>50 mm/yr) can generally be found in the late (2003–2015), middle (1980–2000), and early (1971–1987) period for Beijing, Hebei, and Tianjin, respectively. The land subsidence showed varied correlations to GWSA and abstraction in space and time. Contrast correlations were found between subsidence and abstraction during early (positive) and late (negative) stages, while the GWSA represented a well indicator for subsidence due to its integrated representation of discharge and recharge.

This study also demonstrates the potential and the limitations of remote sensing approaches (satellite gravimetry, InSAR) for understanding the spatial and temporal dynamics of the aquifer system. However, the data absence for some periods (mainly 1984–2002) and the lack of overall land surface deformation data at larger scale may lead to uncertainties associated with the correlation analysis. It is expected that the GRACE Follow-On mission as well as succeeding new SAR satellites can enhance the understanding of the regional aquifer system at multiple spatial and temporal scales.