Satellite signal shows storage-unloading subsidence in North China

integrates satellite, model and field data products to investigate WSD and land subsidence in North China. In the first step, GRACE (Gravity Recovery and Climate Experiment) mass rates are used to show WSD in the region. Next, GRACE total water storage (TWS) is corrected for soil water storage (SWS) to derive groundwater storage (GWS) using GLDAS (Global Land Data Assimilation System) data products. The 10

derived GWS is compared with GWS obtained from field-measured groundwater level to show land subsidence in the study area. Then GPS (Global Positioning System) data of relative land surface change (LSC) are used to confirm the subsidence due to WSD. A total of ∼ 96 near-consecutive months (January 2002 through December 2009) of datasets are used in the study. Based on GRACE mass rates, TWS depletion 15 is 23.76 ± 1.74 mm yr −1 or 13.73 ± 1.01 km 3 yr −1 in the 578 000 km 2 study area. This is ∼ 31 % of the slated 45 km 3 yr −1 water delivery in 2050 via the South-North Water Diversion Project. Analysis of relative LSC shows subsidence of 7.29 ± 0.35 mm yr −1 in Beijing and 2.74 ± 0.16 mm yr −1 in North China. About 11.53 % (2.74 ± 0.18 mm or 1.58 ± 0.12 km 3 ) of the TWS and 8.37 % (1.52 ± 0.70 mm or 0.88 ± 0.03 km 3 ) of the 20 GWS are attributed to storage reductions accompanying subsidence in the region. Although interpretations of the findings require caution due to the short temporal and large spatial coverage, the concurrence of WSD and land subsidence could have adverse implications for the study area. It is critical that the relevant stakeholders embark on resource-efficient measures to ensure water availability, food security, ecological

6045
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Poland et al. (1975) noted a record groundwater-controlled subsidence of > 8.5 m in San Joaquin Valley of California. There are several other reports of subsidence due to long-term groundwater withdrawals in China and elsewhere across the globe (Ireland et al., 1984;Lofgren, 1991;Zhang et al., 2008;Li et al., 2011).
Monitoring land surface change (LSC) for crustal deformation is a critical first-5 step for effective intervention measures (Galloway et al., 1999;Galloway and Burbey, 2011). 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 spaceborne techniques like GPS (Global Positioning System) and InSAR (Interferometry 10 Synthetic Aperture Radar) are now available. The procedural details of several LSC monitoring techniques are discussed by Lofgren (1991), Abidin et al. (2001), Gambolati et al. (2006), Chen et al. (2007) and Liu et al. (2008).
Here, GPS data products of relative LSC are used in combination with hydrological data products from GRACE (Gravity Recovery and Climate Experiment), GLDAS 15 (Global Land Data Assimilation System) and field-measured groundwater level and soil water storage (SWS) to characterize water storage depletion (WSD) and land subsidence in North China.
Launched in March 2002 (Swenson et al., 2006), GRACE is a twin satellite system that orbits the Earth in tandem at a separation distance and altitude of ∼ 200 and 20 ∼ 450 km, respectively (Syed et al., 2008;Longuevergne et al., 2010). Changes in the satellite separation distance due to variations in gravity pull are measured by an onboard K-band microwave to the nearest nanometer. At any point in time, the satellite orbit position is tracked by GPS receivers and related with the inter-satellite distance for gravity fields (Wahr et al., 2004;Rodell et al., 2009). Then mass redistributions within 25 the Earth (e.g., atmospheric surface and ocean bottom pressures and terrestrial water, snow/ice storage) are inferred from the time-variable gravity fields (Chambers, 2009;Longuevergne et al., 2010). After removing non-hydrologic effects, GRACE gravity field residuals closely track terrestrial water storage changes derived from field measure-6046 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ments or model estimates at various spatial scales (Ramillien et al., 2008;Strassberg et al., 2009;Famiglietti et al., 2011;Moiwo et al., 2011).
GLDAS mission models global land surface energy states (Rodell et al., 2004), including soil water, runoff and evapotranspiration (ET). GLDAS land surface models include Mosaic, Noah, Community Land Model, Variable Infiltration Capacity and Catch-5 ment Model (Hogue et al., 2005). As the accuracy levels of GRACE and GLDAS data products are comparable (Rodell et al., 2009;Longuevergne et al., 2010), their integrative use can improve water storage analysis Immerzeel et al., 2010).
In this study, GRACE total water storage (TWS) is corrected for groundwater storage (GWS) using GLDAS-derived SWS and compared with GWS from field-measurements to show WSD and land subsidence in North China. Then the storage-unloading subsidence is confirmed by analysis of GPS data product of relative LSC in the region. This study extends current applications of GRACE/GLDAS data products well beyond the conventional field of water storage analyses into crustal deformation analyses. This 15 will deepen existing knowledge on the implications of long-term storage depletion for human subsistence. It will also add greater urgency to the need to develop resourceefficient strategies for a truly sustainable world.

Materials and method
2.1 Study area 20 The study area, North China, lies between longitudes 107.1-119.2 • E and latitudes 34.6-41.7 • N in an area of ∼ 578 000 km 2 . In addition to Beijing and Tianjin, the area variously covers six provinces ( Fig. 1, top plate) and has a total population exceeding 250 × 10 6 people. The region accounts for over 45 % of China's grain production and 15 % of its gross domestic product (Cao et al., 2013;Feng et al., 2013;Tang et al., 25 2013).

6047
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The average annual ET (985 mm yr −1 ) in cultivated lands is over two times the average precipitation (500 mm yr −1 ), sustained largely by irrigation Moiwo et al., 2011). The thick Quaternary alluvial deposits (> 1000 m) are poorly drained with average hydrologic conductivity of > 100 m d −1 (Foster et al., 2004;Hu et al., 2010;Cao et al., 2013). Thus the several decades of high water use (especially in agriculture and 5 industry) has pushed water exploitation far beyond precipitation recharge in the region (Kendy et al., 2007;Cao et al., 2013). The difficulties in the water balance accounting and the increasing gap between water use and recharge have caused surface water depletion and groundwater level loss of 1-2 m yr −1 (Yang et al., 2006;Hu et al., 2010). This suggests that pumping rates are not balanced by increased recharge and/or de-10 creased discharge due to groundwater irrigation Cao et al., 2013).
The storage-depletion depression cones of 40-50 m in shallow and 60-80 m in deep aquifer systems are increasing concerns for land subsidence (Xue et al., 2005;Li et al., 2011). While groundwater recovery is cited in recent years (Cui et al., 2009), there are no confirmed reports of lessening water use or increasing precipitation/recharge in the 15 region Cao et al., 2013). This suggests uncertainties in the storage dynamics due to difficulties in accounting for water use in the region. There is thus the need to refine water use strategies towards greater efficiency and sustainability in the region. 20 Pumping wells generate stress disturbances which propagate through aquifer systems and cause pressure/head loss at magnitudes dependent on the hydrogeologic conditions (Gambolati et al., 2006). Stress change due to pumping could result in the compaction of hydrogeologic formations and land subsidence (Waltham, 2002). Subsidence, which is the response of compressible hydrogeologic formation to changes in 25 fluid pressure within the formation, is quantified in terms of effective stress as (Poland and Davis, 1969;Poland, 1984;Galloway et al., 1998):

Basic concept
where σ T is total stress (M L −1 T −2 ), σ e is effective stress (M L −1 T −2 ), and p w is porewater pressure (M L −1 T −2 ). The three common compaction stresses are gravitational stress, hydrostatic stress and dynamic seepage stress (Lofgren et al., 1968). While hydrostatic stress is largely 5 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 ; M L −1 T −2 ), of sediments and water at a reference plane below the water table is the unit weight (γ m ; M L −2 T −2 ) times the thickness (z 1 ; L) of moist sediments above the water table plus the unit weight 10 (γ b ; M L −2 T −2 ) times the thickness (z 2 ; L) of buoyant saturated sediments below the water table (Poland and Davis, 1969;Poland, 1984): in which g is average specific gravity of the grain deposits (-), r s is average specific retention of the moist grains (-), 15 n is average porosity of the aquifer deposits (-), and γ w is specific weight of water (M L −2 T −2 ). The functional connections of the above equations with 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 detailed by Poland and Davis (1969) and Galloway et al. (1998).
Since storage largely occurs as groundwater and soil water (especially in semiarid regions), most of GRACE analyses have focused on these storage components (Syed et al., 2008;Strassberg et al., 2009;Moiwo et al., 2011). Thus the time-variant change in soil water storage (∂Sws/∂t; L T −1 ) and in groundwater storage (∂Gws/∂t; L T −1 ) is related to the time-variant change in GRACE-derived total water storage (∂Tws/∂t; Note that the terms in Eq. (3) are running averages for the given region of interest. Generally, ∂Tws/∂t can be computed from GRACE month-to-month gravity fields inverted for TWS anomaly (TWSA) (Rodell et al., 2009;Longuevergne et al., 2010). Storage anomaly is the residual storage content at a given time with respect to that 5 at a reference epoch. Then storage change is the difference in storage anomalies between successive time steps (Moiwo et al., 2011). Similarly, GLDAS global SWS or field-measured SWS can be inverted for ∂Sws/∂t (Immerzeel et al., 2010;Moiwo et al., 2011). Then groundwater levels can be inverted for ∂Gws/∂t using storage coefficients (Feng et al., 2013;Tang et al., 2013); specific yield (S y ) for unconfined aquifer 10 and specific storage (S s ) for confined aquifer (Syed et al., 2008). Generally, TWS well corrected for SWS can track field-monitored GWS especially in semiarid hydrologic conditions (Strassberg et al., 2009). Sustained groundwater overdraft in North China since the 1978 land reforms (Yang and Tian, 2009) can potentially induce aquifer compaction and land subsidence. The change in the position of mass 15 along the axis of gravity occurs when compaction/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.
Thus the trends in SWS, GWS and TWS are likely to suggest subsidence or uplift, which either way needs verification by LSC analysis. In this study, GRACE TWS is 20 corrected for SWS to derive GWS. The GRACE-corrected GWS anomaly (CGWSA) is then compared with field-measured GWS anomaly (GWSA). While negative trends in GWSA and CGWSA could suggest compaction and subsidence, positive trends suggest expansion and uplift of the crustal formations. To verify the LSC due to WSD in the study area, GPS data products of relative LSC are analyzed for limited region of

Data acquisition and processing
The GRACE, GLDAS and GPS data products are used in combination with fieldmeasured SWS and GWS. Except for GRACE (which spans from April 2002 to December 2009), the datasets cover January 2002 through December 2009. The GRACE release-05 data products from CSR (Center for Space Research), JPL (jet Propulsion 5 Laboratory) and GFZ (German Research Center for Geosciences) used are available at http://geoid.colorado.edu/grace/dataportal.html.
The GRACE monthly solutions are filtered for leakages (using an averaging kernel for the study area), destriped for north-south trending errors and corrected for glacial isostatic adjustments and other non-hydrological effects (Swenson et al., 2006;Paul-10 son et al., 2007;Immerzeel et al., 2010). The data are smoothed at 200 km Gaussian half-width, fitted with the Wahr et al. (2004) error bar and truncated to within one degree of the study area. The fields are afterwards spatially averaged (based on the averaging kernel) to create time series of TWS anomaly (TWSA), from which ∂Tws/∂t is computed. The spatial average of the GRCAE-estimated monthly TWSA from the 3 centers 15 for 2002-2009 is plotted in the bottom right plate of Fig. 1.
The average of the GLDAS-estimated optimal fields of SWS from the 5 model is used in this study (Rodell et al., 2009). The GLDAS data product is verified with fieldmeasured SWS data in the study area. Like GRACE, the GLDAS data field is truncated to within one degree of the study area and spatially averaged (based on the averaging 20 kernel) to create time series of SWS anomaly (SWSA) and then ∂Sws/∂t. The GLDAS data products used are also available at http://gdata1.sci.gsfc.nasa.gov/daac-bin/G3/ gui.cgi?instance_id=GLDAS10_M.
The groundwater data are from 205 fairly distributed monitoring wells in the study area ( Fig. 1, top plate). About one-third of the wells monitor groundwater level in the Storage coefficients of the unconfined (S y ) and confined (S s ) aquifers are 0.012-5 0.264 and 0.5639-1.0668 × 10 −3 mm −1 , respectively (Feng et al., 2013;Tang et al., 2013). The mean S s is 0.07829 × 10 −3 mm −1 , which is ∼ 33 % (0.235) that of S y -see http://en.cgs.gov.cn/ or http://www.cgs.gov.cn/ for details. GWSA is derived by multiplying the separately-interpolated GWS with storage coefficient for the respective aquifer systems. The GWSA for the separate aquifers are then aggregated and interpolated 10 for ∂Gws/∂t 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 15 groundwater level and storage. However, these effects are not critical in this study as the analyses are long-term averages.
The GPS data of relative LSC (frame IGS08) are 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 20 GAMIT/GLOBK for daily loosely-constrained solutions of long-term surface deformation Palano et al., 2012) 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 for analyzing GPS measurements, primarily to study crustal deformations.
Generally, GRACE accuracy increases at higher spatial/temporal scales (Swenson et al., 2006), so the analysis is scaled up from monthly to yearly cycles. Here, monthly cycles are the month-to-month time series for 2002-2009 spatially averaged over the study area. Average monthly cycles are the time series averaged temporally along the 12 months of the years in 2002-2009 and spatially over the study area. The seasonal and average seasonal cycle definitions follow the definitions for the monthly and average monthly cycles. In the study area, 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
Because GLDAS generally underestimates SWS, the standard deviation of the fivecontributing products is used as a measure of uncertainty (Rodell et al., 2009). The contributing model-estimated standard deviation is 2.16 mm month −1 , < 7 % of the average value and therefore acceptable for this study. The GLDAS SWS is in good agree-10 ment with field-observed SWS (adjusted for irrigation) with R 2 = 0.899 (Fig. 1, bottom right plate). The SWS data from 80 monitoring sites in the study area are spatially averaged on monthly basis (for January 2002 through December 2009) and plotted against the GLDAS-estimated SWS in Fig. 1.
The GRACE dataset is processed based on the one-to-zero (from center-to-15 peripheral) averaging kernel method (Swenson et al., 2006) to minimize leakage. The kernel is spatially smoothed and truncation as a function of the Stokes coefficient, noise filtering and Gaussian convolution. Amplitude damping due to the kernel averaging is corrected following the method used by Strassberg et al. (2009). Least squares analysis suggests that average random errors in the GRACE monthly data are < 5.32 mm, 20 which is < 5 % of the average error. Also the standard deviation of the three-contributing GRACE data products is used as a measure of uncertainty. The estimated standard deviation is 1.64 mm, < 6 % of the average value. Because errors generally cancel out with higher spatial/temporal averaging and storage change is derived from storage anomaly but also for clarity, error bars are only plotted for monthly storage anomalies. 25 Also linear trends and regression equations are shown on the plots for a quick grasp of the storage dynamics.

6053
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 2 plots the time series of the monthly GWS, SWS and TWS anomalies and changes in the study area (see details in Table 1). The range of storage change ("Dif" in Table 1) is the difference between the maximum and minimum amplitudes. It is smallest 5 for SWS and highest for GWS, suggesting that the storage is mainly driven by change in GWS (Feng et al., 2013;Tang et al., 2013). The time series (especially storage anomalies in Fig. 2a, c and e) are similar in phase and the amplitude range smallest for SWS. The trend-lines and regression functions for all the storage terms suggest declining storage in the region (Liu et al., 2006;10 Zheng et al., 2010). Because up-scale averaging reduces errors by limiting outlier effects (Swenson et al., 2006), the datasets are represented at seasonal and annual cycles to more clearly show storage dynamics in the study area.

Seasonal cycles
The time series of the seasonal storage anomalies and changes are plotted in Fig. 3 15 (see details in Table 1). Due to suppressed outlier effects by the up-scale averaging, the amplitudes are smaller and phases more similar at the seasonal than at the monthly cycles. Based on the trend lines, GWS ( Fig. 3a and b), SWS ( Fig. 3c and d) and TWS ( Fig. 3e and f) generally decline in the study area. Like in the monthly cycles, the amplitude change is smallest for SWS. The storage trends show not only seasonality, 20 but are also significant at p < 0.05 and α = 0.05. The seasonality is attributed to the hydro-climatic and agronomic conditions in the study area (Moiwo et al., 2011).

Yearly cycles
The amplitudes of the yearly storage dynamics in Fig. 4 are further reduced (due to the up-scale averaging), with the smallest amplitude range still for SWS. The phases of the storage anomalies (Fig. 4a, c and e) and changes (Fig. 4b, d and f) are a lot more similar at this scale. Like in the other cycles, all the trends in TWS, SWS and GWS terms are negative (Feng et al., 2013;Tang et al., 2013). The low SWS variability is attributed to irrigation during non-raining cropping periods. In fact, SWS change is generally minimal in irrigated regions (Rodell et al., 2009). 5 The annual storage trends are significant at p < 0.01 and α = 0.05. The storage anomalies (Fig. 4a, c and e) and changes (Fig. 4b, c and f) show that the storage is highest in 2004. Although the trends suggest an apparent recovery in 2006-2008, the overall storage remains negative. This contradicts recent reports of groundwater recovery in the region (Cui et al., 2009). The continuous WSD could cause deformation 10 of the geo-matrix and land subsidence in the region (Zheng et al., 2010). Figure 5 depicts the dynamics of the average monthly storage terms in the study area. GWS increases from January through February and then decreasing through March. It is lowest in June and again rebounds through December (Fig. 5a). The corresponding 15 GWS change decreases from January through April, recovers through September and again decreases through December (Fig. 5b). The rates of precipitation (storage input), irrigation (storage input/output) and ET (storage output) greatly influence storage dynamics in the study area (Hu et al., 2010;Cao et al., 2013a;Tang et al., 2013).

Average monthly cycles
The SWS and TWS terms have different amplitudes but similar phases. The phases 20 generally decline from January through May/June, increase to peak values in August and again decline through December ( Fig. 5c and e). The corresponding trends in SWS and TWS change track a trough-shaped curve for the period from May through September. The trends in storage change are not so pronounced before May and after September ( Fig. 5d and f). As GRACE is relatively less sensitive to short-term storage 25 changes, TWS ( Fig. 5b and f) is a bit out of phase with measured GWS (Fig. 5a and e) in the study area.

6055
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 6 presents the dynamics of the storage terms averaged seasonally for [2002][2003][2004][2005][2006][2007][2008][2009]. The dynamics of GWS and SWS are similar but altogether different from that in TWS; further pointing out the low sensitivity of GRACE to short-term storage changes (Swenson et al., 2006). GWS anomaly declines from winter through summer, before 5 rebounding in autumn (Fig. 6a). That of SWS declines from winter to spring and rebounds through autumn (Fig. 6c). Also GWS as well as SWS change ( Fig. 6b and d) decline from winter through spring and rebound through autumn. The trends in GRACE average seasonal TWS anomaly (Fig. 6e) and change (Fig. 6f) track a forward-leaning S-shaped-like curve. It is lowest in spring and highest in sum-10 mer; depicting the overall storage characteristics in the region. The storage dynamics is the combined effects of the water use, precipitation and ET in the region Cao et al., 2013).

Groundwater storage anomaly
GWS anomaly derived from field observations (GWSA) is plotted along with that de-15 rived from GRACE/GLDAS datasets (CGWSA) in Fig. 7. Here, storage anomaly is the storage variation relative to the mean (monthly, seasonal or annual) storage for 2002-2009 and then storage change is the difference in two successive time periods. Thus the trend in GWS anomaly in effect reflects the rate of groundwater depletion in the study area. The degree of the trends (given by the trend-lines and regression func-20 tions) in GWSA (Fig. 7, left plots) and in CGWSA (Fig. 7, right plots) are negative at the monthly, seasonal and yearly scales, suggesting that there is storage loss in the study area (Yang et al., 2006;Hu et al., 2010).
Water budget analysis based on the GRACE data shows WSD at 23.76 ± 1.74 mm yr −1 , the equivalent of 13.73 ± 1.01 km 3 yr −1 for the 578 000 km 2 study area.
Sustained WSD over long periods of time could cause land subsidence and further limit the storage and water supply (Gambolati et al., 2006;. To 10 determine that there could be subsidence due to sustained WSD in the study area, GPS data products is used to analyze LSC in the region.

Land subsidence
Because of the dense population, high agro-industrial expansion and semiarid climatic conditions, water use in North China exceeds natural recharge rate (Hu et al., 2010;15 Wang et al., 2011;Feng et al., 2013). Since there are hardly any occurrences of earthquakes or large-scale earth-faults in the region, land surface deformation could only be caused by abstractions of groundwater, hydrocarbons or coal (Liu et al., 2008). Thus GPS data product of relative LSC is used to analyze for land subsidence due to loss of water storage in the region. 20 The GPS data analysis suggests the occurrence of subsidence at an estimated rate of 7.29±0.35 mm yr −1 in the vertical component of the IGS08 station in Beijing (Fig. 7g). Note that Beijing is one of the zones with severe groundwater depression cones in China. The vertical dip of hydraulic gradient from moderate-to-worse WSD zones is 4-30 % in Beijing (Wang et al., 2011). Groundwater accounts for > 61 % of water supply 25 in Beijing, 26 % in Tianjin, 80 % in Hebei province and 58 % in Shanxi province (Cao et al., 2013).

6057
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | As the degree of subsidence is loosely related to the magnitude of groundwater depression cone , subsidence in Beijing (where depression cones are severe) could be worse than the average for the rest of North China (where there are many areas with little/no depression cones). Analysis of well data shows that storage loss is 62.41 % higher in Beijing than the average for the 578 000 km 2 North China 5 study area. Given the spatial averaging, a similar hydrogeologic conditions is assumed across the study area. Then a relational analysis suggests the occurrence of subsidence at an average rate of 2.74 ± 0.16 mm yr −1 in the entire study area. Because the assumption here ignores the basis for complex subsidence processes, the estimated subsidence only conservative and must be treated with caution. Irrespectively, it could guide future subsidence analyses at specific strategic locations of the study area (e.g., Tianjin, Shijiazhuang and Baoding) with WSD comparable to that in Beijing. Shi et al. (2006) noted that since land reforms in the 1970s, subsidence has been occurring at the rate of 2.97 mm yr −1 in North China Plain. Skeletal matrix can be rearranged under increased effective stress exceeding pre-consolidation stress in fine 15 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 20 and tends to occur under very large increases in effective stress (Poland, 1994). Thus high storage depletion (Wang et al., 2011;Cao et al., 2013) in North China could cause any of the above subsidence conditions.

Discussions
Groundwater accounts for some 80 % of irrigation Hu et al., 2010) 25 and 20 % of industrial and domestic (Feng et al., 2013) water use in North China. The year-on-year groundwater level decline (of 1-2 m) suggests that withdrawal exceeds recharge (Cao et al., 2013). This unsustainable reliance on groundwater is also reflected in the widespread WSD in the region (Yang et al., 2006;Wang et al., 2011).
Based on GPS data analysis, subsidence due to WSD is 7.29 ± 0.35 mm yr −1 in Beijing and 2.74 ± 0.16 mm yr −1 in the rest of North China. Assuming that the subsidence 5 is all drainable water, it is the equivalent of 0.12 ± 0.01 km 3 for the 16 080 km 2 Beijing area and 1.58 ± 0.09 km 3 for the 578 000 km 2 North China region. There is also a considerable water unloading via mineral/coal mining operations in the region (Feng et al., 2013). As WSD will continue into the foreseeable future due to the lack of reliable alternatives, pumping-induced subsidence could worsen environmental conditions in the 10 region in the long run (Wang et al., 2011).
Droughts, degenerated water/wetland ecosystems and earthquakes are variously associated with long-term WSD. Pumping-drained groundwater level loss of 100 m in Las Vegas Valley (USA) in the 1950s resulted in ∼ 2 m subsidence, severely damaging infrastructures in the region (Bell, 1981). Intense pumping in the 1940s in Antelope 15 Valley caused ∼ 2 m subsidence in Lancaster and > 1 m subsidence in Rogers Lake region, USA (Galloway et al., 1998;Hoffmann et al., 2003). Groundwater drawdown of > 250 m during 1960-2010 triggered the 2011 earthquake of 5.1 moment magnitude in Lorca, Spain (González et al., 2012).
Numerous reports have pointed to WSD and/or subsidence in a number pf regions 20 North China (Li, 2003;Liu et al., 2006;Wang et al., 2011;Feng et al., 2013). Consistent with other studies (Foster et al., 2004;Kendy et al., 2007;Tang et al., 2013), this study shows WSD in North China. The storage loss (1.58 ± 0.09 km 3 ) due to subsidence is ∼ 12 % of TWS loss (13.73 ± 1.01 km 3 ) in the (578 000 km 2 ) study area. This underscores the need for water-wise strategies which will enhance recharge 25 and limit discharge in the region (Cui et al., 2009;Cao et al., 2013). There is SWS, GWS and TWS loss in the study area. After correcting for storage loss due to subsidence, annual WSD in the region is 21.02 ± 1.58 mm (12.15 ± 0.91 km 3 ) for TWS and 15.52 ± 0.76 mm (8.97 ± 0.44 km 3 ) for GWS. Of course root-zone SWS is 6059 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | largely insensitive to subsidence. The small SWS change (Moiwo et al., 2011) is due to intensive irrigation in the region (Han et al., 2008;. The storageunloading subsidence mainly occurs in the deep aquifer systems (Xue et al., 2005;Zhang et al., 2008;Li et al., 2011;Cao et al., 2013).
The storage trends are significant at the seasonal and yearly scales and have clear 5 seasonality. Currently, there are no significant storage recovery measures such as to reduce groundwater use or increase recharge in the region (Cui et al., 2009). In fact, GRACE signal can detect storage recoveries driven by such measures (Feng et al., 2013). GWS is lowest and SWS highest during summer precipitation months (Fig. 5), suggesting that precipitation is most effective for root-zone SWS in the region (Han 10 et al., 2008). High SWS in the other seasons is due to irrigation . Over-reliance on groundwater causes a steady decline in GWS (Yang et al., 2006). GRACE captures the overall conditions as storage loss (Fig. 5e). While the average monthly trends suggest that the lowest storage is in summer, the average seasonal trends clearly reflect the prevailing hydro-climatic conditions. 15 The monthly, seasonal and yearly storage anomalies are negative for GWS derived from well data (GWSA) and that from GRACE/GLDAS data (CGWSA) for 2002-2009 (Fig. 7). As argued in Eq. (2) and the subsequent discussions, GWS from GRACE/GLDAS should track that from well data. The storage trends at the various scales are negative, suggesting storage loss. Even historical records of water use in 20 the region suggest WSD (Foster et al., 2004;Zhang et al., 2008). Excessive WSD can cause land subsidence (Shi et al., 2006;Wang et al., 2011;. The condition for subsidence is verified using GPS data of relative LSC. The analysis suggests the occurrence of subsidence at 7.29 ± 0.35 mm yr −1 in Beijing and 2.74 ± 0.16 mm yr −1 in the North China study area (Fig. 7g). Subsidence driven by WSD could 25 have disastrous implications for ecological sustainability, water availability, food security and social stability in this pivotal region. Thus there is the need to reorient water-use policies towards greater efficiency and sustainability in this semiarid region.
Our current understanding of WSD and the approaches towards sustainable water use is limited. This study contributes to the broader efforts of deepening this knowledge for greater efficiency and sustainability in water use. It integrates satellite, model and field data products to show WSD and subsidence in North China. Error analysis shows that 5 the results are not the effects of data biases or artifacts of the analysis.
The rate of TWS loss in the study area is 23.76±1.74 mm yr −1 (13.73±1.01 km 3 yr −1 ), which is 30.52 % of the planned south-north water transfer (45 km 3 yr −1 ) in 2050. Cumulative WSD over 2002-2009.08 ± 13.92 mm (109.87 ± 8.05 km 3 ), which is 244.12 % of the slated water transfer. Subsidence due to WSD is 7.29 ± 0.35 mm yr −1 10 in Beijing and 2.74 ± 0.16 mm yr −1 in North China. About 11.53 % (2.74 ± 0.18 mm; 1.58 ± 0.12 km 3 ) of TWS and 8.37 % (1.52 ± 0.70 mm; 0.88 ± 0.03 km 3 ) of GWS are attributed to storage reductions accompanying subsidence in the region. Excessive groundwater use causes storage depletion and subsidence in the long-term. Subsidence further limits GWS by reducing crustal pore spaces. Because GRACE merely detects mass change, it is used in combination with LSC data to show subsidence in this study. The concurrence of groundwater depletion and subsidence could grow into water crisis in the region.
North China is a pivotal agricultural, industrial and political base where WSD and subsidence could destabilize millions in the region. It is critical for relevant stakeholders 20 to embark on greater efficiency in water use. Use of alternative water resources could mitigate the WSD and subsidence. Aquifer recharge, saltwater use and inter-basin water transfer could enhance GWS, food security and social stability in the region.