Implications of Groundwater Depletion for Aquifer Geomatrix Deformation and Water Availability

Worsening water storage depletion contributes to environmental degradation, land subsidence and earthquake or even disrupts food production/security and social stability. There is also the need for efficient water use strategies in populated regions, especially when such regions also have intensive agricultural and industrial activities. The North China Plain (NCP) is one such region which is not only the seat of power, but also a major agricultural and industrial with a severe water storage depletion. Thus this study integrates satellite, model and field data products to investigate water storage depletion and land subsidence in the Beijing Environs of NCP. In the first step, GRACE (Gravity Recovery and Climate Experiment) mass rates are analyzed for water storage depletion in the region. Next, GRACE total water storage (TWS) is corrected for soil moisture storage (SMS) to derive groundwater storage (GWS) using GLDAS (Global Land Data Assimilation System) data products. The derived GWS is compared with GWS obtained from field-measured groundwater level to show water storage depletion in the study area. Then GPS (Global Positioning System) data of relative land surface change are used to show land subsidence due to water storage depletion. A total of ~96 near-consecutive months (Jan. 2002 through Dec. 2009) of datasets are used in the study. Based on GRACE mass rates, TWS depletion is 36.54±1.74 mm yr -1 or 6.34±0.29 km 3 yr -1 for the 169 000 km 2 study area. Analysis of relative land surface change shows the occurrence of land subsidence at 7.29±0.35 mm yr -1 in the Beijing Environ of NCP. About 7.50% (2.74±0.18 mm yr -1 or 0.46±0.03 km 3 yr -1 ) of the depletion in TWS and 5.25% (1.52±0.07 mm yr -1 or 0.26±0.01 km 3 yr -1 ) of that in GWS are attributed to storage reductions due to the land subsidence. Storage loss in the region justifies the current south-north water transfer efforts in the region. The concurrence of water storage depletion and land subsidence could have adverse implications for the hydrology, ecology, food security and social stability of the region. It is important to devise efficient measures to avert the negative effects of water storage depletion in the study area.


Introduction
Groundwater is a vital source of water supply for humans, agriculture, industries and high-value ecosystems [1,2]. Current global groundwater withdrawal is estimated at 750−800 km 3 yr -1 [3], providing some 50% of drinking water needs, 40% of industrial water demand and 20% of irrigation water supply [4]. Along with the benefits (e.g., increased food production and socio-economic growth) are also the consequences (e.g., storage depletion, quality deterioration, ecological degradation, land subsidence, seawater intrusion and water conflicts) of groundwater exploitation [1,5]. Excessive groundwater exploitation over long periods could induce land subsidence [6][7][8] or even earthquake [8,9].
Groundwater withdrawal reduces pore pressure and also concurrently increases effective stress; a condition that causes aquifer compaction, land subsidence or earthquake [9][10][11]. Subsidence is the elastic/inelastic compaction of crustal materials due to changes in stress dynamics [12,13]. Subsidence due to groundwater unloading is the compaction Deformation and Water Availability of clay-silt/sand-gravel deposits [1,14]. Aquifer recharge restores lost pore-water pressure that in turn induces elastic rebound of aquifer materials [14].
Despite its significance or even because of it, current global rates of groundwater development are unsustainable [1,14,16]. Under certain conditions, physical changes within aquifer/crustal systems due to storage loss can be permanent, leading to unrecoverable loss of storativity. Groundwatercontrolled crustal deformations like groundwater depression cones and land subsidence have been reported in Arizona/Texas [6], Mexico City [17], Jakarta [18], Ravenna [19], Taiwan [20] and Lorca/Spain [9]. Other crustal deformations and the related effects are detailed by González et al. [9] for the Los Angeles Basin in USA and by Rodolfo and Siringan [21] for East Asia. Poland et al. [22] detected a record groundwater-controlled subsidence of > 8.5 m in San Joaquin Valley of California. There are several such reports on land subsidence due to long-term groundwater withdrawals in China and elsewhere across the globe [11,[24][25].
Monitoring land surface change (LSC) for crustal deformation is the critical first step towards developing effective subsidence intervention measures [12,26]. Several techniques exist for LSC monitoring, including spirit leveling, piezometric/reservoir pressure measurement, geodimeter measurement and extensometer measurement. Also radioactive marker techniques, modeling techniques and space-borne techniques like GPS (Global Positioning System) and InSAR (Interferometry Synthetic Aperture Radar) are now available. The procedural details of several LSC monitoring techniques are discussed by [8,10,18,23,27]. Here, GPS data product of relative LSC is used in combination with hydrological data products from GRACE (Gravity Recovery and Climate Experiment), GLDAS (Global Land Data Assimilation System) and field-observed groundwater level and soil moisture storage (SMS) to characterize water storage depletion and land subsidence in the Beijing Environs of North China Plain (NCP).
GRACE is a twin-satellite system launched in March 2002 to monitor the dynamics in global water storage [28]. The twin satellites orbit the Earth in tandem at a separation distance of ~200 km and initial altitude of ~450 km [29,30]. Changes in the separation distance between the twin satellites due to variations in gravity pull are measured by an onboard K-band microwave to the nearest nanometer. At any given point in time, the orbit positions of the twin satellites are tracked by GPS receivers and related with the inter-satellite distance driven by gravity fields [31,32]. Then mass redistributions within the Earth (e.g., due to atmospheric surface pressure, ocean bottom pressure, and changes in storage of terrestrial water and snow/ice) are inferred from the time-variable gravity fields [30,33]. After removing nonhydrologic effects, the GRACE gravity field residuals are closely related with terrestrial water storage changes derived from field measurements or model estimates at spatial scales ≥ 160 000 km 2 [34][35][36][37][38][39][40][41].
GLDAS uses land surface models (e.g., Mosaic, Noah, Community Land Model, Variable Infiltration Capacity and Catchment Model) to estimate global land surface energy states, including soil moisture, runoff and potential evapotranspiration (ET) [38,39]. Due to the high complementarity of the data types in terms of spatial/temporal coverage and accuracy, GRACE and GLDAS data products are integratively used to analyze terrestrial water storage dynamics and to develop sustainable management measures [30,37,40]. In this study, GRACE total water storage (TWS) is corrected for groundwater storage (GWS) using GLDASderived SMS and then compared with GWS from fieldmeasurements to show water storage depletion and land subsidence in the Beijing Environs of NCP. Next, the storage-unloading subsidence is verified by GPS data products of relative LSC in the region. The study depicts a novel demonstration of the application of GRACE/GLDAS data products in crustal deformation analysis, going beyond the conventional water storage analysis for such data products. Thus the study not only deepens our current understandings of the implications of storage depletion, but also strengthens the need for developing sustainable water resources use strategies. This is significant not only for research, planning and policy judgement, but also for food production and social stability in the region and beyond.

Study Area
The study area, the Beijing Environs of NCP, lies between longitudes 113.8-119.8 ºE and latitudes 37.5-42.6 ºN and covers an area of ~169 000 km 2 (Fig. 1, top plate). With a population of over 112 million, NCP accounts for over 45% of China's grain production and 15% of its gross domestic product [2,5,41,42].
Potential ET (985 mm yr -1 ) in cultivated lands is about two times the average annual precipitation (500 mm yr -1 ), which is sustained by intensive irrigation [15,37]. The thick Quaternary alluvial deposits (> 1000 m) are a little poorly drained with average hydraulic conductivity of > 100 m d -1 [2,20,47]. Because of intensive water use (especially in agriculture and industry), water exploitation in the region far exceeds precipitation recharge [2,45]. There is sustained surface water depletion and water table head loss (of about 1-2 m yr -1 ) driven by the ever-widening gap between water use and recharge [15,20]. The pumping rates (due to groundwater irrigation) are neither balanced by increased recharge nor decreased discharge, worsening water shortage in the region [2,15].
Depression cones of 40-50 m and 60-80 m have developed respectively in the shallow and deep aquifer systems, raising concern for storage depletion effects beyond water shortage for it various uses [7,25]. While groundwater recovery is cited in recent years [46], there is no confirmed report of decreasing water use or increasing precipitation/recharge in the region [2,15]. The difficulties in the water budget accounting deepens the uncertainties in the storage dynamics, hence the need for efficient and sustainable use of the water resources in the study area.

Basic Concept
Pumping wells generate stress disturbances that propagate through aquifer systems and cause pressure/head loss at magnitudes dependent on the hydrogeologic conditions. Stress change due to pumping could result in the compaction of hydrogeologic formations and land subsidence [13]. Subsidence, which is the response of compressible hydrogeologic formation to changes in fluid pressure within the formation, is quantified in terms of effective stress as [6,47,48]: where T σ is total stress (ML-1T-2), e σ is effective stress (ML-1T-2), and w p is pore-water pressure (ML-1T-2).
The three common compaction stresses are gravitational stress, hydrostatic stress and dynamic seepage stress [23]. While hydrostatic stress is largely neutral, gravitational stress and dynamic seepage stress are additive in effect and together change the void ratio and mechanical properties of aquifer deposits. The combined effect, known as geostatic pressure or total stress ( T σ ; ML-1T-2), of sediments and water at a reference plane below the water table is the unit weight ( m = − + ; in which g is average specific gravity of the grain deposits (-), s r is average specific retention of the moist grains (-), n is average porosity of the aquifer deposits (-), and w γ is specific weight of water (ML-2T-2). Details on the links between subsidence and elastic/inelastic storage coefficients (e.g., specific storage/yield, aquifer thickness, etc.) and other aquifer properties (e.g., pore-water pressure, hydraulic conductivity, etc.) are documented by Poland [47] and Galloway et al. [12].
The degree 2, order 0 (C20) of GRACE spherical harmonic (SH) coefficients can be estimated from satellite laser ranging [28]. Then the other GRACE SH coefficients approximated to degree and order 60 for variations in surface mass density ( 0 S ∆ ) in an area ( 0 R ) are calculated as [5]: After glacial isostatic adjustment [49], north-south stripe filtering, high-frequency noise smoothing and local kernel averaging [28], the surface mass variations are inverted for TWS anomalies comparable to field-measured or modelestimated values [5]. Then piezometric head measurements are inverted for groundwater storage anomaly via storage coefficients [5,42] -specific yield ( y S ) for unconfined aquifer and specific storage ( s S ) for confined aquifer [29]. Also as storage (especially in semiarid regions) mainly occurs as groundwater and soil moisture [29,35,37], GRACE TWS analyses are largely limited to these storage components.
A number of studies have pointed out to aquifer compaction and land subsidence in NCP due to sustained groundwater overdraft since the 1978 land reforms [7,50,51]. The trends in SMS, GWS and TWS can serve as a vital signal for potential land subsidence or uplift in a region, either of which can be verified through LSC analysis. The change in the position of mass along the axis of gravity occurs when compaction or subsidence moves mass towards the center of mass of the Earth, causing change in gravity. However, this effect is small with respect to the mass changes addressed in this study and therefore not considered.
In this study, GRACE TWS is corrected for SMS to derive GWS and then the GRACE-corrected GWS compared with field-measured GWS. While negative trends in the water storage anomalies suggest compaction and subsidence, positive trends on the other hand suggest expansion and uplift of the aquifer systems. To verify LSC due to water storage depletion in the Beijing Environs study area, GPS data products of relative LSC are analyzed in relation to the derived TWS anomaly.

Data Acquisition and Processing
The GRACE, GLDAS and GPS data products are used in combination with field-measured SMS and GWS. Except for GRACE (which dataset spans from Apr. 2002 to Dec. 2009), the datasets cover the period from Jan. 2002 through Dec. 2009. The GRACE release-05 data products from CSR (Center for Space Research), JPL (Jet Propulsion Laboratory) and GFZ (German Research Center for Geosciences) used in this study are available at http://geoid.colorado.edu/grace/dataportal.html.
The concept of GRACE data generation and application in hydrology is today quite well documented. GRACE signal detects water storage dynamics across regions with spatial resolutions of 160 000 km 2 or more [32], which makes it applicable in the 169 000 km 2 Beijing Environs study area. Here, GRACE monthly solutions are corrected for glacial isostatic adjustments and other non-hydrological effects, destriped for north-south errors, smoothened at 200 km Gaussian half-width for high-frequency noise and truncated within 1° (one degree) of the study area over a local averaging kernel [5,28,40,49]. The GRACE-derived hydrological fields are afterwards spatially averaged to create time series of TWS anomaly. A linear trend of the spatially averaged GRCAE-estimated monthly TWS anomalies from the 3 geo-centers for the period 2002-2009 is plotted in the bottom left plate of Fig. 1.
The average of the GLDAS-estimated optimal fields of SMS from the 5 surface energy models is used in this study [32]. The GLDAS data product is verified using fieldmeasured SMS data in the study area. Like GRACE, the GLDAS data field is truncated to within 1° of the study area and spatially averaged (based on the local averaging kernel) to create time series of SMS anomaly. Long-term GLDAS data products are available at http://gdata1.sci.gsfc.nasa.gov/daacbin/G3/gui.cgi?instance_i d=GLDAS10_M.
The groundwater data are from 138 fairly distributed monitoring wells in the study area (Fig. 1, top plate). About one-third [46] of the monitoring wells are for groundwater level in the confined aquifer and the rest [92] for the unconfined aquifer. The data, also for the period from Jan. 2002 through Dec. 2009, are the averages of three observations (8 th , 18 th and 28 th ) per month. The data are converted into GWS anomaly using storage coefficients. To do this, the monthly groundwater data are interpolated (Inverse Distance Weighted Interpolation) in ArcMap separately for the unconfined and confined aquifers.
For the period 2002-2009, groundwater level range across the study area is 66 m to -26 m for the unconfined aquifer and 14 m to -87 m for the confined aquifer. Storage coefficients of the unconfined ( Sy ) and confined ( Ss ) aquifers are 0.012−0.264 and 0.5639-1.0668 × 10 -3 mm -1 , respectively [5,42]. The mean Ss is 0.07829 × 10 -3 mm -1 , which is ~33% (0.235) that of Sy (see http://water.cgs.gov.cn for details). The GWS anomaly derived from field measurements is the product of the interpolated value and the storage coefficient of the respective aquifer systems. This is then summed up and processed in the same way as described for GRACE and GLDAS data products.
Note that the concept of storage coefficient is a simplified assumption of instantaneous drainage of aquifer systems with dropping pressure heads as it ignores delayed drainage effects. While storage coefficients could vary spatially, the change as a function of depth to groundwater can result in non-linear relationship between in-situ groundwater level and storage. However, these effects are not critical in this study as the analysis focuses on long-term averages.
The (frame IGS08) GPS data product of relative LSC is part of the Global Navigation Satellite System (GNSS) or International GNSS Service (IGS) available at http://www.unavco.org/crosscutting/cc-data.html. The data product is processed in GAMIT/GLOBK for daily, looselyconstrained solutions of long-term surface deformation [9,52] in the region. GAMIT/GLOBK, developed by Massachusetts Institute of Technology (MIT) and sponsored by National Science Foundation (NSF), is a comprehensive package of programs used for analyzing GPS measurements, primarily to study crustal deformations.
Generally, GRACE signal accuracy increases at higher spatial/temporal resolutions [28], thus the analysis focuses on seasonal and annual cycles. However, both average monthly and seasonal cycles are included to show GRACE limitation at short temporal scale. Here, storage anomaly is the storage variation relative to the mean (seasonal or annual) storage for the period 2002-2009. Average monthly cycle is the temporal average for the 12 months in the year and for the 9 years of 2002-2009, along with the spatial average for the study area. Then seasonal cycle is the season-to-season time series for 2002-2009 spatially averaged over the study area. The average seasonal and the yearly (annual) cycles are respectively defined after the average monthly and seasonal cycles. In the Beijing Environs, the period from December to February is winter, March to May is spring, June to August is summer and September to November is autumn.

Uncertainty/Bias
In this study, uncertainty in the GLDAS-derived SMS is measured as the standard deviation of the five-contributing products [32]. The analysis shows that the 5-contributing model-estimated standard deviation is 2.16 mm/month, which is < 7% of the average value. Because that GLDAS generally underestimates SMS (especially for irrigated regions), the above level of accuracy is sufficient for a study of this nature. The SMS data from the 60 monitoring sites are spatially averaged (for Jan. 2002 through Dec. 2009) and plotted against the GLDAS-estimated SMS in Fig. 1 (bottom  right plot). The GLDAS-derived SMS agrees well (R 2 = 0.899) with the field-observed SMS adjusted for irrigation (Fig. 1, bottom right plot).
The GRACE data processing is based on the one-to-zero (center-to-peripheral) averaging kernel function [28] for the study area, which reduces leakage. Induced bias, leakage and amplitude damping as a result of the SH coefficient truncation, glacier isostatic adjustment, destripping, kernel averaging and Gaussian convolution are corrected to reduce uncertainty and bias [35]. Least squares analysis suggests that random errors in the GRACE monthly data are < 5.32 mm, which is < 5% of the average error. A further measure of uncertainty in the GRACE data is based on the standard deviation of the 3-contributing GRACE data products. The estimated standard deviation is 1.64 mm, < 6% of the average value. Although error bars are not shown on the plots, the depicted linear trends and regression equations on the timeseries plots sufficiently highlight the accuracy of the storage anomalies.

Seasonal Cycles
Time series of the seasonal storage anomalies for GWS, SWS and TWS are plotted in Fig. 2 (also see details in Table  1). The anomaly range (denoted as "Dif" in Table 1) is the difference between the maximum and minimum storage amplitudes. It is smallest for SMS and highest for GWS, suggesting that storage in the study area mainly occurs as GWS [5,42]. The seasonal amplitudes are similar to the seasonal phases because up-scaling (which is upward averaging) limits outlier effects on trend curves. Based on the trend lines, GWS (Fig. 2, top plot), SMS (Fig. 2, middle plot) and TWS (Fig. 2, bottom plot) decline in the study area. The storage trends are not only seasonally driven, but are also significant at p < 0.05 and α = 0.05. The seasonality of the storage dynamics is attributed to the prevailing hydroclimatic and agronomic conditions in the region [37]. Note that Min is the minimum amplitude; Max is the maximum amplitude; Dif is the difference between Max and Min amplitudes; STD is the amplitude standard deviation; Avgmon is the average monthly storage anomaly; and Avgseas is the average seasonal storage anomaly.

Yearly Cycles
The amplitudes of the yearly storage dynamics in Fig. 3 are much less (due to the up-scale averaging), also with the smallest range for SMS and the highest for GWS. Also the phases of the storage anomalies ( Fig. 2; top, middle and bottom plots) are more similar than at the seasonal scale. Like the seasonal cycles, the trends in TWS, SMS and GWS are negative [5,42]. The low SMS variability in the study area is attributed to irrigation during non-raining farming seasons. In fact several other studies have noted a generally low SMS anomaly in irrigated farmland regions [32,37].
The annual storage trends (which are significant at p < 0.01 and α = 0.05) suggest that storage in the region is highest in 2004. Although the trends suggest an apparent recovery in 2006−2008, the overall trend is negative; contradicting recent reports of groundwater recovery in the region [46]. The continuous water storage depletion raises concern for crustal deformation and land subsidence in the study area [53].  Figure 4 depicts the dynamics of the average monthly storage anomalies in the study area. GWS increases from January through February and then decreases through March. It is lowest in June and rebounds through December (Fig. 4,  top plot). The storage dynamics is strongly influenced by the rates of precipitation (storage input), irrigation (storage input/output) and ET (storage output) in the study area [2,20,42].

Average Monthly Cycles
The SMS and TWS anomalies differ in amplitude, but similar in phase. The phases generally decline from January through May/June, increase to peak values in August and again decline through December (Fig. 4, top and middle plots). Since GRACE is less sensitive to short-term storage anomalies, TWS anomaly (Fig. 4, bottom plot) is a bit out of phase with measured GWS anomaly (Fig. 4, top plot) in the study area. The above discrepancies could also be due to errors in the GRACE, GLDAS and/or observation data or in the data processing method used.  Figure 5 plots the dynamics of the storage anomalies averaged seasonally for the period 2002-2009. Although the GWS and SMS anomalies are similar, both are different from TWS anomaly. This further confirms that GRACE signal is less sensitive to short-term storage variations [28]. The GWS anomaly declines from winter through summer, before rebounding in autumn (Fig. 5, top plot). That of SMS declines from winter to spring and rebounds through autumn (Fig. 5, middle plot).

Average Seasonal Cycles
The trends in the GRACE-derived average seasonal TWS anomaly (Fig. 5, bottom plot) depict horizontal S-shaped (or a breakthrough) curve. It is lowest in spring and highest in summer, which depicts the overall storage characteristics in the region. The storage anomaly depicts the combined effects of water use, precipitation and ET in the study area [2,15].

Groundwater Storage Anomaly
The GWS anomaly derived from observation data is plotted side-by-side that derived from GRACE/GLDAS data products in Fig. 6. The trend in GWS anomaly reflects the rate of groundwater depletion in the study area. The magnitude of the trend (given by the trend-lines and regression equations) in the field-measured GWS anomaly (Fig. 6, top and middle plots on the left) and that in GRACE-GLDAS-derived GWS anomaly (Fig. 6, top and middle plots on the right) are negative at both seasonal and annual scales, suggesting an overall storage loss in the study area.
Water budget analysis based on GRACE data shows water storage depletion at the rate of 36.54±1.74 mm yr -1 , the equivalent of 6.34±0.29 km 3 yr -1 for the 169 000 km 2 Beijing Environs study area. Also analysis of GLDAS data products shows that 20.75% (7.58±0.36 mm yr -1 ; 1.28±0.06 km 3 yr -1 ) of that loss is in SMS. Groundwater and other storage sources account for the remaining 79.25% (28.96±1.38 mm yr -1 or 4.89±0.23 km 3 yr -1 ) of the storage loss. Feng et al. [5] noted that about 67% of storage loss in North China is from aquifer storage. The estimated contribution of GWS to TWS loss in the region is 76.84% (28.28±1.02 mm yr -1 or 4.79±0.23 km 3 yr -1 ). The cumulative storage loss for the 8- Thus as noted by Strassberg et al. [35] for a semiarid highland in the US, the overall storage dynamics in the study area is little affected (< 3%) by surface water use. Sustained water storage depletion over a long period could lead to land subsidence, an event which limits water storage and supply [9]. To determine the occurrence of land subsidence in the Beijing Environs due to water storage depletion, GPS data products of relative LSC are analyzed.

Land Subsidence
Because of the dense population, high agricultural/industrial activities and severe semiarid conditions, water use in NCP exceeds natural recharge [5,20,55]. As there are no earthquakes or large faults in the region, land deformation could only be due to the abstraction of groundwater, hydrocarbons or coal [27]. Thus GPS data of relative LSC can be used to determine land subsidence due to water storage loss in the region.
The GPS data analysis suggests the occurrence of subsidence at an estimated rate of 7.29±0.35 mm yr -1 (Fig. 6, bottom plot) in the vertical component of the IGS08 station in Beijing (Fig. 1, top plate). Note that Beijing has one of the largest groundwater depression cones in China. Since land reforms in the 1970s, there is sustained subsidence at the rate of 2.97 mm yr -1 in NCP [50]. The vertical dip of hydraulic gradient from the moderate-to-worst water storage depletion zone in Beijing is 4−30% [54]. Groundwater accounts for > 61% of water supply in Beijing, 26% in Tianjin, 80% in Hebei province and 58% in Shanxi province [2], all of which form the NCP.
The magnitude of subsidence is loosely related to the degree of groundwater depletion [9], which is severe for the Beijing Environs,, including Tianjin and the contiguous prefectures (Fig. 1, top plate). For the spatial averaging, it is assumed that the hydro-geologic conditions are uniform across the Beijing Environs study area. This makes the above estimate conservative since it ignores the complex processes of land subsidence. Irrespectively, the analysis lays the basis for future applications of GRACE/GLDAS data products in the study of not only water storage loss, but also other related processes such as subsidence and earthquakes.
There is the possibility of skeletal matrix rearrangement under increased effective stress exceeding pre-consolidation stress in fine grain deposits, resulting in reduced porosity that is permanent and non-recoverable. Skeletal matrix generally compresses under decreased effective stress (not exceeding pre-consolidation stress), also resulting in reduced porosity. However, this deformation is recoverable under stresses in the elastic range. Crushing of granular components of the matrix is possible in sand and especially diatomaceous deposits, but this is rare and tends to occur under very large increases in effective stress [48]. The high storage depletion [2,54] in NCP could result in the above subsidence conditions.

Discussions
Groundwater accounts for some 80% of irrigation [15,20] and 20% of industrial and domestic [5] water use in North China. The year-on-year groundwater level decline (1-2 m) suggests that withdrawal exceeds recharge in the region [2]. The over reliance on groundwater is also reflected in the widespread water storage depletion in the region [54,55,57].
Based on GPS data analysis, subsidence due to water storage depletion is 7.29±0.35 mm yr -1 in the Beijing Environs of NCP. Assuming that the subsidence is all drainable water, it is the equivalent of 1.23±0.06 km 3 for the 169 000 km 2 study area. There is also a considerable water unloading via mineral/coal mining operations in the region [5]. As water storage depletion will continue into the foreseeable future due to the lack of reliable alternatives, the environmental conditions could further deteriorate due to pumping-induced effects such as subsidence [54].
Droughts, degraded water/wetland ecosystems and earthquakes are variously associated with long-term water storage depletion [54]. Pumping-drained groundwater level loss of 100 m in Las Vegas Valley (USA) in the 1950s resulted in ~2 m subsidence, severely damaging engineered infrastructures in the region [57]. Intense pumping in the 1940s in Antelope Valley caused some 2 m of subsidence in Lancaster and over 1 m of subsidence in Rogers Lake region, USA [12,58]. Groundwater drawdown of > 250 m during 1960-2010 triggered the 2011 earthquake of 5.1 moment magnitude in Lorca, Spain [10].
Numerous reports have pointed to water storage depletion and/or subsidence in North China [5,25,56,59]. Consistent with other studies (42)(43)(44)(45), this study shows water storage depletion in the Beijing Environs of NCP. The storage loss due to subsidence is ~20% (1.27±0.06 km 3 yr -1 ) of TWS loss (6.34±0.29 km 3 yr -1 ) in the 169 000 km 2 study area. This underscores the need for water-efficient strategies that increase recharge and limit discharge in the region [2,46].
There is SMS, GWS and TWS loss in the Beijing Environs study area of NCP. After correction for storage loss due to subsidence, annual water storage depletion in the region is 21.02±1.58 mm yr -1 (3.55±0.27 km 3 yr -1 ) for TWS and 15.52±0.76 mm yr -1 (2.62±0.13 km 3 yr -1 ) for GWS. Of course root-zone SMS is largely insensitive to subsidence. The low SMS anomaly [37] is due to intensive farm irrigation in the region [15,60]. Thus storage-unloading subsidence exists mainly in deep aquifer systems in the region [2,7,11,25].
The storage trends are significant at the seasonal and annual cycles and with clear seasonality. Storage recovery measures (which can reduce groundwater use or increase recharge) detectable by GRACE signal [5] are lacking in the region [46]. GWS is lowest and SMS highest during summer precipitation months (Fig. 4), suggesting that precipitation is only effective for root-zone SMS in the region [60]. The high SMS during the other seasons is due to intensive farm irrigation [15]. Over-reliance on groundwater causes a steady decline in storage [55]. GRACE captures the overall hydrologic conditions as storage loss (Fig. 4, bottom plot). While the average monthly trends suggest that storage is lowest in summer, the average seasonal trends clearly reflect prevailing hydro-climatic conditions in the study area.
The monthly, seasonal and annual storage anomalies for 2002−2009 are negative for GWS derived from well data (Fig. 6, top and middle plots on the left) and that from GRACE/GLDAS data (Fig. 6, top and middle plots on the right). As in Eq. (3) and the subsequent discussions, GWS derived from GRACE/GLDAS should track that derived from field observations. The storage trends at the various cycles are negative, again confirming storage loss. Even historical records of water use in the region suggest water storage depletion [11,43,44]. Excessive storage depletion can induce land subsidence [9,50,54].
The potential for subsidence in the study area is verified by GPS data of relative LSC. The analysis suggests the occurrence of subsidence at the 7.29±0.35 mm yr -1 in the Beijing Environs study area (Fig. 5, bottom plot). Subsidence due to water storage depletion could have disastrous implications for ecological sustainability, water availability, food security and social stability in the region and beyond. This calls for strong water-use policies directed towards greater water use efficiency in the semiarid region.

Conclusions
The current understanding of water storage depletion and the approaches towards sustainable water use is limited. This study deepens our existing knowledge on hydrological events in to sustainable of water resources. The study integrates satellite, model and field data products to show water storage depletion and land subsidence in the Beijing Environs of NCP. Error analysis suggests that the results of the study are not the effects of data bias or artifacts of the analysis.
TWS loss in the 169 000 km 2 study area is 36.54±1.74 mm yr -1 (6.34±0.29 km 3 yr -1 ), justifying the need for the southnorth water transfer project. Cumulative water storage depletion for 2002−2009 is 190.08±13.92 mm (49.40±2.35 km 3 ), which is 109.77% of the slated annual water transfer (45 km 3 ) in 2050. Subsidence due to water storage depletion is 7.29±0.35 mm yr -1 in the Beijing Environs of NCP. About 7.50% (2.74±0.18 mm; 0.46±0.03 km 3 ) of TWS and 5.25% (1.52±0.07 mm; 0.26±0.01 km 3 ) of GWS are attributed to storage reductions accompanying land subsidence in the region. Excessive water use over the long-term can induce storage depletion, which could in turn trigger subsidence. Subsidence could further limit GWS by crushing down available pore spaces in the crustal matrix. As GRACE signal detects mass change, it is used in combination with GPS data to show subsidence in this study. The concurrence of groundwater depletion and land subsidence increases the risk of water crisis in the study area.
The NCP is a region with intensive agricultural, industrial and political activities. This suggests that water storage depletion with subsidence could be a disastrous combination and with a destabilizing effect on the region. It is therefore important that the relevant stakeholders embark on efficient water use strategies in the region. The use of alternative water resources, aquifer recharge and inter-basin water transfer could enhance storage, food security and social stability in the region.