Detecting seasonal and long-term vertical displacement in the North China Plain using GRACE and GPS

In total, 29 continuous Global Positioning System (GPS) time series data together with data from Gravity Recovery and Climate Experiment (GRACE) are analysed to determine the seasonal displacements of surface loadings in the North China Plain (NCP). Results show significant seasonal variations and a strong correlation between GPS and GRACE results in the vertical displacement component; the average correlation and weighted root-meansquares (WRMS) reduction between GPS and GRACE are 75.6 and 28.9 % respectively, when atmospheric and nontidal ocean effects were removed, but the annual peak-topeak amplitude of GPS (1.2–6.3 mm) is greater than the data (1.0–2.2 mm) derived from GRACE. We also calculate the trend rate as well as the seasonal signal caused by the mass load change from GRACE data; the rate of GRACEderived terrestrial water storage (TWS) loss (after multiplying by the scaling factor) in the NCP was 3.39 cm yr−1 (equivalent to 12.42 km3 yr−1) from 2003 to 2009. For a 10year time span (2003 to 2012), the rate loss of TWS was 2.57 cm yr−1 (equivalent to 9.41 km3 yr−1), which is consistent with the groundwater storage (GWS) depletion rate (the rate losses of GWS were 2.49 and 2.72 cm yr−1 during 2003–2009 and 2003–2012 respectively) estimated from GRACE-derived results after removing simulated soil moisture (SM) data from the Global Land Data Assimilation System (GLDAS)/Noah model. We also found that GRACEderived GWS changes are in disagreement with the groundwater level changes from observations of shallow aquifers from 2003 to 2009, especially between 2010 and 2013. Although the shallow groundwater can be recharged from the annual climate-driven rainfall, the important facts indicate that GWS depletion is more serious in deep aquifers. The GRACE-derived result shows an overall uplift in the whole region at the 0.37–0.95 mm yr−1 level from 2004 to 2009, but the rate of change direction is inconsistent in different GPS stations at the−0.40–0.51 mm yr−1 level from 2010 to 2013. Then we removed the vertical rates, which are induced by TWS from GPS-derived data, to obtain the corrected vertical velocities caused by tectonic movement and human activities. The results show that there are uplift areas and subsidence areas in NCP. Almost the whole central and eastern region of NCP suffers serious ground subsidence caused by the anthropogenic-induced groundwater exploitation in the deep confined aquifers. In addition, the slight ground uplifts in the western region of NCP are mainly controlled by tectonic movement (e.g. Moho uplifting or mantle upwelling).

Abstract.In total, 29 continuous Global Positioning System (GPS) time series data together with data from Gravity Recovery and Climate Experiment (GRACE) are analysed to determine the seasonal displacements of surface loadings in the North China Plain (NCP).Results show significant seasonal variations and a strong correlation between GPS and GRACE results in the vertical displacement component; the average correlation and weighted root-meansquares (WRMS) reduction between GPS and GRACE are 75.6 and 28.9 % respectively, when atmospheric and nontidal ocean effects were removed, but the annual peak-topeak amplitude of GPS (1.2-6.3 mm) is greater than the data (1.0-2.2 mm) derived from GRACE.We also calculate the trend rate as well as the seasonal signal caused by the mass load change from GRACE data; the rate of GRACEderived terrestrial water storage (TWS) loss (after multiplying by the scaling factor) in the NCP was 3.39 cm yr −1 (equivalent to 12.42 km 3 yr −1 ) from 2003 to 2009.For a 10year time span (2003 to 2012), the rate loss of TWS was 2.57 cm yr −1 (equivalent to 9.41 km 3 yr −1 ), which is consistent with the groundwater storage (GWS) depletion rate (the rate losses of GWS were 2.49 and 2.72 cm yr −1 during 2003-2009 and 2003-2012 respectively) estimated from GRACE-derived results after removing simulated soil moisture (SM) data from the Global Land Data Assimilation System (GLDAS)/Noah model.We also found that GRACEderived GWS changes are in disagreement with the groundwater level changes from observations of shallow aquifers from 2003 to 2009, especially between 2010 and 2013.Al-though the shallow groundwater can be recharged from the annual climate-driven rainfall, the important facts indicate that GWS depletion is more serious in deep aquifers.The GRACE-derived result shows an overall uplift in the whole region at the 0.37-0.95mm yr −1 level from 2004 to 2009, but the rate of change direction is inconsistent in different GPS stations at the −0.40-0.51mm yr −1 level from 2010 to 2013.Then we removed the vertical rates, which are induced by TWS from GPS-derived data, to obtain the corrected vertical velocities caused by tectonic movement and human activities.The results show that there are uplift areas and subsidence areas in NCP.Almost the whole central and eastern region of NCP suffers serious ground subsidence caused by the anthropogenic-induced groundwater exploitation in the deep confined aquifers.In addition, the slight ground uplifts in the western region of NCP are mainly controlled by tectonic movement (e.g.Moho uplifting or mantle upwelling).

Introduction
Global Positioning System (GPS) was used to monitor crustal motion, especially in the vertical or height component due to its large amplitude, which has been used to study surface loading caused by mass change.Site-position time series recorded by continuous GPS arrays have revealed the vertical displacement variations resulted from trend or seasonal distribution of mass in a region or global changes, e.g. a change of continental water (Bevis et al., 2005; van Dam et L. Wang et al.: Detecting seasonal and long-term vertical displacement al., 2007;Wahr et al., 2013), ice (Sauber et al., 2000;Khan et al., 2010;Nielsen et al., 2013), snow (Heki, 2001;Grapenthin et al., 2006), ocean (van Dam et al., 2012;Wahr et al., 2014), and atmospheric mass (van Dam et al., 1994;Boehm et al., 2007).
On the global scale, terrestrial hydrologic mass exchanges that cause significant large-scale loading occur between the oceans, continents, and atmosphere at seasonal and interannual timescales.On the local scale, for the inter-annual and long-periodic change in the hydrologic cycle, what significantly affects loading is large anthropogenic disturbances on groundwater extraction and artificial reservoir water impoundment and other climate-driven factors, e.g.natural floods and droughts (Chao et al., 2008;Rodell et al., 2009;Feng et al., 2013;Joodaki et al., 2014;Wang et al., 2014a).The global-scale mass variations closely related to changes in terrestrial water storage (TWS) are observed by the Gravity Recovery and Climate Experiment (GRACE) satellite mission, while the surface elastic displacement can be estimated if the load and rheological properties of the Earth were known (Farrell, 1972).The majority of previous loading studies solved the three components of crustal motion by the adopting joint analysis of GRACE time-variable gravity field coefficients and GPS data (Kusche and Schrama, 2005;van Dam et al., 2007;Fu and Freymueller, 2012;Fu et al., 2013).In principle the loading effects caused by the majority of mass redistributions near the Earth's surface are from water, atmospheric, and ocean transports on daily to interannual timescales (Kusche and Schrama, 2005).The contribution to the surface displacement made by the variations of the atmospheric and ocean can be reasonably modelled by using global atmospheric surface pressure data and space geodetic data respectively.Thus, after removing the loading effects of the atmosphere and ocean, GRACE-derived displacement and GPS data allow for the detection of changes in the Earth's larger hydrological storage systems.
In general, changes in TWS capacities depend on precipitation and human consumption.Variations in TWS may be related to precipitation, which is strongly driven by climate, and can be simulated from global water and energy balance models (Syed et al., 2008).This is related to soil texture and root depth in the case of soil water storage; e.g.soil moisture (SM) and vegetation canopy storage can be derived from the Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004), the WaterGAP Global Hydrology Model (WGHM) (Döll et al., 2003), and the Community Land Model (CLM) (Oleson et al., 2013), surface water storage (e.g.water in rivers, lakes, reservoirs, and wetlands can be derived from WGHM while snow or ice can be derived from WGHM and CLM), and naturally occurring (i.e.climate-driven) aquifer storage (e.g.groundwater predicted by WGHM and CLM).Variations in TWS may also be caused by man-made factors, such as water withdrawals for irrigation purposes (Döll, 2009) and dam construction for power generation and navigation (Wang et al., 2011).These changes in TWS can be observed in situ (i.e.groundwater level and impounded water level).The TWS integrated by all the variations can lead to the overall changes in crust displacement.
This study focuses on the crustal deformation of the North China Plain (NCP) (Fig. 1), which is one of the most uniformly and extensively altered areas by human activities in the world (Tang et al., 2013).The NCP is one of the world's largest aquifer systems and supports an enormous exploitation of groundwater.Overexploitation of groundwater has seriously affected agriculture irrigation, industry, public supply, and ecosystems in the NCP.Previous studies used GRACE data, land surface models, and well observations to provide insight on groundwater depletion in the NCP (Su et al., 2011;Zhong et al., 2009;Feng et al., 2013;Moiwo et al., 2013;Tang et al., 2013;Huang et al., 2015).Liu et al. (2014) has discussed loading displacement in the NCP before.Only five GPS stations' (i.e.BJFS, BJSH, JIXN, TAIN, and ZHNZ) data are used in their work.Although they calculated the seasonal amplitudes, phases, and trends of vertical displacement from GRACE and GPS, the atmospheric and non-tidal ocean-loading effects were not removed in the Liu's work (Liu et al., 2014); i.e. they added the Atmosphere and Ocean De-aliasing Level-1B (AOD1B) solution (GAC solution) back to the GRACE spherical harmonic solutions.
Hydrol.Earth Syst.Sci., 21, 2905Sci., 21, -2922Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/2905/2017/ Here, we use GRACE and data from 29 GPS sites to study the seasonal and long-term loading displacement due to dynamic hydrological processes and groundwater-derived land subsidence in the NCP.In contrast to previous focus study (Liu et al., 2014), the most obvious difference between our results and their work is we removed other loading effects (e.g.atmospheric and non-tidal ocean) in order to reflect the seasonal and long-term displacement caused by TWS loads better.Additionally, we discuss long-term trends due to mass changes revealed by GRACE measurements and its impacts on tectonic vertical rates evaluations.

GRACE data
The GRACE mission design makes it particularly useful for time-variable gravity studies.GRACE was jointly launched by NASA and the German Aerospace Center (DLR) in March 2002 (Tapley et al., 2004).The Level-2 gravity products consist of sets of complete spherical harmonic (Stokes) coefficients out to some maximum degree and order (typically l max = 120) averaged over monthly intervals.When considering large-scale mass redistribution in the Earth system on a timescale ranging from weekly to interdecadal, it is reasonable to assume that all relevant processes occur in a thin layer at the Earth's surface (Kusche and Schrama, 2005).In this analysis, we assume that the gravitational and geometrical response of the Earth can be described by the theory of Farrell (1972), where the loads' Love numbers only depend on the spherical harmonic degree.Thus, the elastic displacements due to the surface mass change can easily be represented in terms of spherical harmonic coefficients for the gravity field and load Love numbers, k l , l l , and h l (Wahr et al., 1998;Kusche and Schrama, 2005).Level-2 products are generated at several project-related processing centres, such as the Center for Space Research (CSR) at the University of Texas, GeoForschungsZentrum in Potsdam, Germany, and the Jet Propulsion Laboratory in California.The mass estimates (TWS and sea level) show very good agreement among these products (Fu and Freymueller, 2012;Wahr et al., 2014;Wang et al., 2014b).
This study used monthly sets of spherical harmonic (Stokes) coefficients from GRACE RL05 (i.e.Release 5) gravity field solutions generated from the CSR, spanning from February 2003 to April 2013.Each monthly GRACE field consisted of a set of Stokes coefficients, C l,m and S l,m , up to a degree and order (l and m) of 60.In fact, the GRACE Stokes coefficients ("GSM" coefficients denoted by the GRACE Project) have had modelled estimates of the atmospheric and oceanic mass signals removed.Thus, the GRACE coefficients include the full terrestrial water storage signal with remaining atmospheric and oceanic signals due to errors in the respective models (Swenson et al., 2008).
Generally, using the GRACE AOD1B products can add back the de-aliasing atmospheric and non-tidal oceanic effects to the GRACE data.However, we would like to reduce the environmental loading contributions to the GRACE and GPS observations, if we study on the accurate interpretation of displacement due to TWS loading.Thus, we analysed the effects of non-tidal ocean variations and atmospheric loading on the GRACE model and GPS coordinates; please see Supplement Sect.S1 for details.
We replaced the GRACE C 20 coefficients with C 20 coefficients inferred from satellite laser ranging (Cheng et al., 2013).Due to the fact that the reference frame origin used in the GRACE gravity field determination is the Earth's centre of mass (CM), GRACE cannot determine the degree-1 terms variations in the Earth's gravity field (geocenter motion).Here, we used degree-1 coefficients as calculated by Swenson et al. (2008) to determine the position of the CM relative to the centre of figure (CF) of the Earth's outer surface.We applied the post-processing method described in Swenson and Wahr (2006) to remove north-south stripes.We adopted 250 km as the averaging radius to implement Gaussian smoothing, a technique that suppresses errors at high degrees (Wahr et al., 1998;van Dam et al., 2007).Stokes coefficients resulted from A et al. (2013) were used to remove contributions from glacial isostatic adjustment (GIA).The contribution of GIA is about 0.28-0.33mm yr −1 and nonseasonal in the NCP, which is small and non-seasonal; therefore, their impact on the seasonal results discussed in this paper would be minimal, even if they were not removed.
The spatial pattern of TWS, shown in Fig. 2, was obtained from monthly GRACE mass solutions for NCP and surrounding regions between spring, 2003, and spring, 2013.An obvious negative trend was identified localized over North China, including some of the north-western regions (i.e.Shanxi province) and north-eastern regions (i.e.Liaoning province).The TWS changes derived from the GRACE data show significant loss trends across the entire study area (NCP), specifically in Beijing, Tianjin, Hebei province, and Shanxi province.Previous studies have investigated how much groundwater depletion has caused the GRACE-derived TWS loss in the whole of the NCP (Su et al., 2011;Feng et al., 2013;Moiwo et al., 2013) or in different sub-regions of the NCP (Zhong et al., 2009;Tang et al., 2013;Huang et al., 2015).These investigations, however, did not focus on regional displacement due to seasonal or long-term variations of hydrologic loading.

GPS data
In total, 24 GPS sites from the Crustal Movement Observation Network of China (CMONOC) and five GPS sites from the International GNSS Service (IGS) (Table 1) were analysed in this study (Fig. 1 shows the locations of the GPS stations).From the GPS sites, eight of them were located in the surrounding area of the NCP.Daily values for the upward, eastward, and northward coordinates were determined by GPS data of IGS stations between 2003 and 2013, which is consistent with GRACE time span.The 24 GPS sites of CMONOC provided data from 2010 to 2013.GIPSY/OASIS-II (version 5.0) software was used in point positioning mode to obtain the daily coordinates and covariances; these were used to transform the daily values into ITRF2008 (Altamimi et al., 2011).We estimated this daily frame alignment transformation using a set of reliable ITRF stations (∼ 10 stations each day).In the GPS processing, corrections for solid Earth tides were undertaken, and ocean tide loading effects were corrected using ocean tide model FES2004 with Greens Functions modelled in the reference frame of CM (centre of the mass of the total Earth system) to maintain theoretical consistency and adherence to current IERS conventions (Hao et al., 2016;Fu et al., 2012), but atmospheric pressure loading or any other loading variations (non-tidal ocean loading) with periods > 1 day were not removed.In order to focus on the seasonal and trend feature over the entire observation time, we first smoothed the data to reduce large scatter before using a 3-month-wide moving window filter to remove the short-period terms.
Due to the coseismic displacement of the 2011 M w 9.0 Tohoku earthquake, we estimated and removed offsets (i.e. using the differences of the average values of seven days between before and after earthquake to obtain the coseismic displacement) at those times for the vertical time series of GPS stations at Eastern China.Wang et al. (2011) study re-sults reveal that the coseismic horizontal displacements induced by the earthquake are at the level of millimetres to centimetres in North and Northeast China, with a maximum of 35 mm, but the vertical coseismic and post-seismic displacements are too small to be detected.In order to maintain consistency with the GIA effects presented in the GRACE solutions, we remove GIA effects for all GPS, stations by using Stokes coefficients (l = 100) results computed by A et al. ( 2013) which used the ICE5G ice history and VM2 viscosity profile (Peltier et al., 2004).
Figure 3a shows the time series (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) of daily solutions for IGS GPS sites BJFS, ZHNZ, BJSH, JIXN, and TAIN.The long-term linear trends are mainly dominated by surface mass loading and tectonic processes, and the GPS time series show significant seasonal variations.The peak-topeak seasonal amplitude can be seen to be more than 20 mm, which reflects the strong seasonal mass changes in the NCP.The GRACE data from CSR use model output to remove the gravitational effects of atmospheric and oceanic mass variability from the satellite data before constructing monthly gravity field solutions.In order to compare the displacement from GPS with GRACE, the effects of atmospheric and nontidal oceanic loading on the GPS coordinates needed to be removed.Displacements due to atmospheric loading were calculated using data and programs developed by the GGFC (Global Geophysical Fluid Center) (van Dam, 2010; NCEPderived 6 hourly, global surface displacements at 2.5 • × 2.5 • spacing; http://geophy.uni.lu/ncep-loading.html).These utilized the NCEP (National Center of Environmental Protection) reanalysis surface pressure data set.The 12 h sampling model ECCO (Estimating the Circulation & Climate of the Ocean; http://www.ecco-group.org/) is used to compute the surface displacement driven by non-tidal ocean effects and its spatial resolution is 1 • × 0.3-1.0• , i.e. 1 • longitude (zonal) interval and 0.3 to 1.0 • in latitude (meridian) intervals from Equator to high latitude.An example of the effects of the non-tidal ocean and atmospheric loading in the GPS and GRACE data is provided in the Supplement (Fig. S1).
The displacements caused by atmospheric pressure and non-tidal ocean loading mainly show seasonal fluctuations and no obvious long-term trend during GPS observations (e.g.time series of height from atmospheric and non-tidal ocean loading at IGS sites in Fig. 3b).The annual amplitude is 4.0-4.6 and 0.24-0.42mm for the atmospheric and non-tidal ocean-loading effects respectively, while the semiannual amplitude is about 0.3 and 0.03 mm respectively.But the phases between the atmospheric and non-tidal oceanloading effects have more apparent difference.The results of the seasonal amplitudes and phase fits of vertical displacements, derived by GRACE and GPS for IGS stations between before and after correcting atmospheric and non-tidal ocean, are summarized in Table S1 in the Supplement.

Elastic displacements due to mass loads
GRACE Stokes coefficients (Wahr et al., 1998) and load Love numbers (Farrell, 1972) can be used to estimate the displacement effects in three components (up, north and east) caused by mass load changes.The mathematical relationships (Kusche and Schrama, 2005;van Dam et al., 2007) between the radial surface displacement (Up or Height) and the Stokes coefficients of mass is • (C l,m cos(mϕ) + S l,m sin(mϕ)) where h is the displacement of the Earth's surface in the radial direction at latitude θ (theta) and eastward longitude φ (phi), N max = 60, R is the Earth's radius, P l,m is fully normalized Legendre functions for degree l and order m, C l,m and S l,m are time-variable components of the (l, m) Stokes coefficients for some months, and h l , k l , and l l are the three degree-dependent load Love numbers, which are functions of Earth's elastic property.In this equation we adopted the load Love numbers provided by Han and Wahr (1995).
Similarly, horizontal displacements (North and East) can be calculated using the following equations: where n and e are north and east components of the displacement respectively, both having positive values when the crust moves towards the north and east respectively.As is mentioned in Sect.2.1 above, in order to be consistently comparable to the GPS time series, we corrected the degree-1 components to GRACE-derived mass variations, using Stokes coefficients derived by Swenson et al. (2008).Corresponding to degree-1 contribution to vertical displacement, the value of load Love numbers of the degree-1 in the CF frame should be computed by using Eq. ( 23) in Blewitt (2003).Figure 4 shows an example (site BJFS, JIXN, TAIN, and ZHNZ) of the GRACE-derived vertical (Fig. 4c) and horizontal displacements (Fig. 4a and b) before and after destriping.It can be clearly seen that the maximal amplitude of vertical displacement is several order of magnitude higher than horizontal displacements.In addition, the calculated results using the monthly GRACE model data after destriping show that the effects of TWS (e.g.soil moisture) on surface displacements are seasonal variations and long-term changes on vertical and horizontal components.As most of the stations are located in areas of TWS loss in the NCP (see sites location in Fig. 2), the fact is that the motion is upward (see the positive trend of GRACE-derived vertical in Table 1) during this event (if a load is removed, the site uplifts and moves away from the load, e.g.Wahr et al., 2013).Identified horizontal displacements are important as they constrain the location of load changes (Wahr et al., 2013;Wang et al., 2014c).The displacement of the ZHNZ site is upwards and to the south (see the negative trend of the ZHNZ north component in Fig. 4a) due to the mass loss almost due north of the site.Correspondingly, the displacement of the TAIN site is upwards and to the south-east (see the negative trend of the north component and the positive trend of the east compo-nent of the TAIN site in Fig. 4a and b) caused by the mass loss located to the north-west of the site, based on the use of GPS horizontals for loading studies from Wahr et al. (2013).

GRACE-derived seasonal variations and comparison with GPS measurements
Using Eq. ( 1) and GRACE-derived Stokes coefficients, the vertical displacements at the GPS sites in the NCP and its surrounding region can be calculated.To focus on these changes, GRACE-derived vertical displacements were computed by fitting a model with a linear trend and annual periodic terms using least-squares method over the entire 11year time span, for a comparison to the seasonal variations observed by GPS (Table 1).Figure 5 shows time series of vertical displacements for GPS sites of IGS stations (BJFS, BJSH, JIXN, TAIN, and ZHNZ).The fitting results show the GRACE-derived (without ADO1B) peak-to-peak annual amplitudes can be more than 2 mm, and the semi-annual amplitude are also visible at these five GPS sites.This reflects the climate-derived seasonal hydrological fluctuations in the NCP.
Compared with GRACE results mainly due to the mass change in seasonal and long-term linear periods, all GPS time series show significant seasonal and long-term trends,  which are mainly dominated by tectonic and hydrological processes.The fitting results (after least-squares fitting) show the peak-to-peak vertical seasonal displacements from GPS time series are to be larger than GRACE-derived results at those GPS sites, and the peak-to-peak seasonal amplitude changes between 5 and 6 mm (Table 1).The results of the comparison between GPS and GRACE-derived seasonal height variations at 24 GPS sites from CMONOC can be seen in Fig. S3.For all the selected GPS sites, the annual component is more dominant than the semi-annual one.The peak-to-peak annual amplitude is 1.2-6.3 and 1.0-2.2mm for the GPS and GRACE solutions respectively, while the semi-annual amplitude is about 1/2-1/3 times of that in annual amplitude.These more consistent seasonal variations of GRACE and GPS height time series reflect the climatederived seasonal hydrologic process; i.e. heavy monsoonal precipitation in the late summer months result in mass loads increase (the maximum negative of vertical amplitude) and largely pumping for agricultural usage in late spring months cause mass loads decrease (the maximum positive of vertical amplitude).The facts show that GPS data relatively has a larger amplitude than GRACE-derived displacements, which does not merely exist in the IGS stations but also almost in all CMONOC stations except SXGX (Table 1).This indicates that GPS has a strong sensitivity for local surface loading.By contrast, because the spatial resolution of GRACE data is limited to approximately 300 km (N max = 60), GRACEderived results are mainly constrained by large-scale areas.This means that GRACE-derived vertical displacements show a small difference between stations due to the results are averaged over scales of several hundred kilometres or more.
Next, we compare GPS observed and GRACE-derived seasonal height variations.The estimated annual amplitudes and initial phases derived from GPS (grey vector) and GRACE (red vector) are shown in Fig. 6.We find that there are many sites where the signals disagree in both amplitude and phase.The annual amplitudes and phases from GRACEderived results are much more spatially coherent than those determined from the GPS heights because GRACE solutions truncate to l max = 60 and the Gaussian filtering was used to lead to smooth out concentrated loads.The phase results of the GRACE-derived displacements show that the annual signal peaks basically appear between September and October, which indicates a maximum load occurs in these 2 months (summer monsoon).However, there are several differences between the results of GPS and GRACE in some sites, more specifically, the signal peaks sometimes appear between August and September according to GPS data.The five signals of sites in the north-western foothills region of NCP agree in phase, while annual amplitudes from GRACE are significantly less than GPS, e.g.NMTK, NMZL, HEYY, HEZJ, and HECC.The cause of most phase inconsistency may be the different spatial resolution of GRACE compared to GPS.That is, GPS measurements can sense the difference between loads very near the site, and loads a bit further away, but GRACE with wavelengths on the order of 300 km reflects this variation at a monthly scale.Another important reason is that a 1-month sampling of GRACE means a phase sampling of 30 • , while a 1-day sampling of GPS means a phase sampling of ∼ 1 • , the different temporal sampling rate caused the inconsistent phase between GRACE and GPS.
With the purpose of quantitatively evaluating the consistency between GPS and GRACE, we remove GRACEderived seasonal displacement from GPS-observed detrended height time series to compute the reductions of weighted root-mean-squares (WRMS) based on the Eq.(2) in van Dam et al. (2007).Correlation between GPS and GRACE-derived seasonal variations and WRMS reduction ratio-of-remove GRACE-derived seasonal displacement from GPS-observed detrended height time series; please see the Table S3 for details.All the selected sites show high correlation (85-99 %, without TJBH site) when atmospheric and non-tidal ocean effects were not removed.By contrast, the seasonal amplitudes and phases from GRACE results are much more spatially coherent than those determined from the GPS heights, caused by the different spatial resolution between them.In addition, we also attempt to calculate GRACE-derived horizontal displacement using Eqs.( 2) and ( 3), and compare it with the GPS measurements.An example (five IGS sites) of the comparison between GRACEderived and GPS-observed horizontal displacements was presented to demonstrate the correlation of seasonal horizontal variation caused by surface hydrological load.Please see the Fig. S4.
4 Long-term uplift caused by TWS loss

Compare groundwater storage (GWS) variations with in situ measurements and previous results
To estimate TWS changes averaged over the NCP, an averaging kernel based on a calculation that using weighted Gaussian convolution to construct monthly time series from GRACE Stokes coefficients described by Eq. ( 5) of Wahr et al. (2014) was used.This method extends the averaging kernel convolution approach (Swenson and Wahr, 2002) by allowing for non-uniform weighting during the convolution.We took the NCP "basin function" from the China provincial boundary grid points and we convolved with a 250 km Gaussian smoothing.We then applied this averaging kernel  2).The estimated results for the time series analysis also include some contributions outside the NCP due to the finite number of harmonic degrees in the GRACE solution (e.g.l max = 60 for CSR solutions).The average kernel in our study is also not an exact unity cover for the entire NCP area; these two factors result in under-or overestimation of the true TWS time series signal.To estimate this "leakage in" signal, a scaling factor method was used to restore the amplitudedamped TWS time series.This method, as described by Wahr et al. (2014), requires the construction of a set of simulated Stokes coefficients, which represents the signal from a uniformly distributed 1 cm water depth change over the NCP.This estimates a water volume = 3.6626 km 3 based on the overall area of basin function (i.e.366260 km 2 ).By applying our GRACE analysis procedure to these simulated Stokes coefficients, we can infer an average water thickness change equal to 0.47 cm for the NCP.Each monthly GRACE estimate of NCP water thickness is then multiplied by a scaling factor = 1 (cm)/0.47(cm) to obtain variations in the total water thickness per area of the NCP.Multiplying the monthly GRACE estimates of NCP water thickness by a scaling factor = 3.6626 (km 3 )/0.47(cm) provides a mass change of the NCP.Table 2 shows the rate of GRACE-derived TWS loss (after multiplying by the scaling factor) in the NCP was 3.39 cm yr −1 from 2003 to 2009; this is equivalent to a volume of 12.42 km 3 yr −1 .For a 10-year time span, the rate was 2.57 cm yr −1 , which is equivalent to a volume of 9.41 km 3 yr −1 .
In our study, the GRACE-based TWS time series covers water change at all depths: surface water storage, soil moisture, snow and groundwater, including anthropogenic effects (i.e.groundwater withdrawal, inter-basin diversion, reservoir and coal transport).To isolate the groundwater contributions, the Noah version of GLDAS, which possesses monthly intervals and spatial resolution of 1.0 • (Rodell et al., 2004), was used to subtract monthly water storage estimates predicted by land surface models.GLDAS generates a series of land surface forcing (e.g.precipitation, surface meteorology and radiation), state (e.g.soil moisture and temperature, and snow), and flux (e.g.evaporation and sensible heat flux) data simulated by land surface models.The GLDAS/Noah model can provide values of snow, vegetation, and all soil moisture layers, but it does not include anthropogenic and climatedriven groundwater depletion.So we isolated GWS variations by subtracting simulated SM data from GLDAS/Noah model from GRACE-derived total TWS.
In order to confirm validation of the results in this study, our GRACE-based estimate was compared with field measurement data of groundwater level (e.g. in situ water table observations) and the results from previous studies, e.g. the reported TWS loss from Zhong et al. (2009), Su et al. (2011), Moiwo et al. (2013), and the reported GWS loss from Huang et al. (2015), Feng et al. (2013), Tang et al. (2013).We have acquired in situ groundwater level measurements (most of groundwater table depth in the shallow unconfined aquifers, available from 2002 to 2013), which are mainly located in the central and eastern plain of the NCP (including the Beijing, Tianjin, and some cities of Hebei, Henan, and Shandong provinces).The data series are obtained from Ministry of Water Resources of China (MWR) (the website of data integration and services in the groundwater resources is available at http://www.groundwater.cn/default.html).We get the area-weighted mean groundwater level change series in the NCP from time series of monthly groundwater table depth changes of 20 cities in our study region.We also collected the daily precipitation data (rainfall amount) for weather stations during the period of 2003-2012 from China Meteorological Data Sharing Service System (CMDSSS) (available at http: //data.cma.cn/site/index.html).Figure 7 shows our GRACEbased estimate is generally consistent with those monthly groundwater level changes observed by monitoring wells after multiplying by mean value of specific yields in the NCP during 2002-2012.
Note that our GRACE-based TWS time series covers all depth of water mass changes and most of the groundwater level changes from observations of shallow aquifers, which showed the long-term mass loss in the NCP from 2003 to 2009, but the rate of this decrease slowed towards the end of 2009 and then increases again after 2010.This difference between the two sub-periods (2004-2009 and 2010-2012) is mainly induced by climate-driven precipitation recharge in NCP (Fig. 7).In addition, comparison between monthly GWS variations estimated from GRACE minus GLDAS/Noah model and in situ groundwater level measurements also confirmed the difference of trend changes of GWS in these two sub-periods (Table 2).The rate of GRACE-derived GWS loss (after multiplying by the scaling factor) in the NCP was 2.49 cm yr −1 from 2003 to 2009 and 2.72 cm yr −1 from 2003 to 2009; this is equivalent to a volume of 9.12 and 9.96 km 3 yr −1 respectively.Our GRACEbased depletion of groundwater was significantly higher than ground-based measurements.Although the shallow GWS increases from 2010 to 2012 due to precipitation recharge in NCP, but the increase of rainfall is difficult to recharge the current more serious depletion of groundwater in the deep aquifers.
In addition, we compare the depletion in TWS or GWS between our results and previous studied results.Table 2 shows our results and compare them with the earlier analysis in different zones where TWS or GWS loss surveys have been published.We found that the trend rate of our GRACE-based GWS in the whole NCP region is in good agreement with that reported by Huang et al. (2015) during 2003-2012and Feng et al. (2013) during 2003-2010, which are estimated from the  2003-2009  2003-2012  2003-2009  2003-2012   a This study −1.62 ± 0.39 −1.23 ± 0.23 −1.17  level 2 Release-05 GRACE data and multiplied by the scaling factor.However, other previous results showed obvious difference between the loss rate of TWS and GWS because these studies used the early versions of the GRACE data, different defined area of NCP, or do not use scaling factor compared with Huang's, Feng's and our study.For instance, Zhong et al. (2009) found a rate of 2.4 cm yr −1 from 2003 to 2007 based on level 2 Release-04 GRACE data in Beijing, Hebei, and Tianjin; Su et al. (2011) calculated TWS and GWS declining at a rate of 1.1 and 0.5 cm yr −1 from 2002 to 2010 respectively, based on level 2 Release-04 GRACE data in Beijing, Tianjin, Hebei, Shandong, and Henan, and Moiwo et al. (2013) estimated a TWS loss rate of 1.68 cm yr −1 from 2002 to 2009 in the vast northern China (i.e. in addition to Beijing and Tianjin, the study area is comprised of 12 other provinces).Although Tang et al. (2013) did not applied scaling factor to restore the amplitude-damped GRACE signal, but they used the latest GRACE products (RL05) and same region of NCP with our study.Thus, Tang et al. (2013) estimated a GWS depletion rate of 0.84 to 1.4 cm yr −1 (2003-2011) is also in good agreement with our estimated result before being multiplied by scaling factor (i.e. the rate of TWS loss was 1.23 cm yr −1 from 2003 to 2012).

Groundwater depletion contributions to long-term uplift
Loading or unloading of the crust from surface mass changes will cause the crust to subside or uplift with different amplitudes.These displacements depend on the amplitude of the load and the distance between the load and the observation point (Farrell, 1972).On this basis, we used GRACE- derived vertical displacements (the method of elastic displacements due to mass loads described by Sect.2.3) to evaluate TWS loss contributions for the evident crustal uplift in the GPS measurements.Time series of monthly predicted vertical surface displacements from GRACE for 25 GPS sites in the NCP were plotted (Fig. 8a).The fitting results (after least-squares fitting) show the trend rate of GRACEderived vertical displacements for the whole region at the 0.37-0.95mm yr −1 level from 2004 to 2009, but the rate of change direction is inconsistent in different GPS stations at −0.40-0.51mm yr −1 level from 2010 to 2013 (Table 1).The smoothed results indicate a rising trend from 2004 to 2012 (Fig. 8a), which represented the TWS loss in the observation time span.Figure 8a also clearly shows mass anomaly due to TWS changes in the vertical component, e.g. a notable negative peak from 2003 to 2005 and subsidence in 2012 (grey background in Fig. 8a).These GRACE-derived long-term height fluctuations mainly include variations in the storage of natural surface water: high storage in wet years and low storage in dry years (Tang et al., 2013), which can be modelled using land surface model output such as those provided by the GLDAS (Rodell et al., 2004).
Figure 8b shows the GRACE-derived height amplitudes after removing the continental water storage signal, which uses the output from the GLDAS/Noah hydrology model.The calculated results show that the contributions of other types of TWS effects (except groundwater) on the surface are small relative to groundwater depletion, and those main loading effects on the amplitudes of seasonal displacement with no obvious long-term trend.Meanwhile, we clearly see that these fluctuations (grey background in Fig. 8a) are almost erased from the GRACE-derived minus GLDAS/Noah vertical displacements, and the obvious continuous uplift are presented in the grey background of Fig. 8b, which is mainly because the contributions from groundwater depletion in the NCP.Contrasts between the seasonal amplitudes phases and trends fit-of-vertical displacement derived by GRACE after removing GLDAS/Noah effects and the original ones.Please see the Table S4.
For the results described above, after the subtraction of the GLDAS/Noah contributions, GRACE-derived heights largely reflect loading effects from the groundwater (natural and anthropogenic factors) and anthropogenic contributions.The anthropogenic impact on mass change was investigated by Tang et al. (2013) for the effect of inter-basin diversion, reservoir, and coal transport distribution on the GRACEderived estimates of groundwater depletion in the NCP.Results from their investigation showed that the trend of anthropogenic contributions was equivalent to 4.83 mm yr −1 water thickness (described by equations in Table 2 of Tang et al., 2013) during 2003-2011 for the whole NCP.This means that there was a large groundwater depletion contribution for the GRACE-derived vertical displacements in long-term uplift.Investigating groundwater withdrawal due to anthropogenic activities (drinking water extraction, agricultural irrigation, and industrial manufacturing) should be of high importance because precipitation data for this area (shown in Fig. 7) indicated no long-term droughts during the GRACE observation period of 2003-2011.

The loading effects of non-tidal ocean and atmospheric variations
As part of the processing performed by the GRACE Project, the GRACE Stokes coefficients (denoted by the GRACE Project as "GSM" coefficients) have had modelled estimates of the atmospheric and oceanic mass signals removed.Thus, the GRACE coefficients include the full effects of terrestrial water storage.The GRACE Project provides the modelled atmospheric and oceanic contributions to the Stokes coefficients in two forms: "GAC" files, which include the global atmospheric and oceanic effects, and "GAD" files, which have had the atmospheric signals over land set to zero.The coefficients in the GAD file therefore represent ocean bottom pressure variations."GAA" files are those CSR GAC files added to the GSM files, and subtracted from the GAD files.
The coefficients in the GAA file represent atmospheric pressure variations.
As is mentioned in Sect. 1, the cause of the difference between our results and Liu's work (Liu et al., 2014) is we removed atmospheric and non-tidal ocean-loading effects while they did not.However, we found that the amplitude of GPS after removing atmospheric and non-tidal ocean-loading effects is still greater than the GRACE while we added the AOD1B (GAC solution) de-aliasing model to the GRACE solutions (i.e.no atmospheric and non-tidal ocean corrected, please see the Table S1).The most obvious difference between our results and Liu's work (Liu et al., 2014) is that they adopt the load Love numbers from Guo et al. (2004) to transform these coefficients into vertical surface displacement estimates.We check the two results of Love numbers (ocean load and atmospheric pressure load) from Guo et al. (2004), there are significant differences between ocean-load and atmospheric pressure-load Love numbers.Meanwhile, we compared k n Love numbers from Guo et al. (2004) (work of Liu et al., 2014) and k n Love numbers from Han and Wahr (1995) (our work) with the k n Love numbers used in ADO1B products (Farrell, 1972) respectively.The different Love numbers have caused the amplitude of the same station from Liu's GRACE-derived vertical displacements much more than GPS and our GRACE results, due to k n from atmospheric pressure-load Love numbers (Guo et al., 2004) significantly larger than Love numbers from Han and Wahr (1995) and Farrell (1972).The detailed analysis of the different Love numbers from Guo et al. (2004), Han andWahr (1995), andFarrell (1972), please see the Sect.S1.Moreover, Table S3 indicates that our correlation results of IGS stations (BJFS, BJSH, JIXN, TAIN, and ZHNZ) are consistent with Liu's work (Liu et al., 2014), indicating that the seasonal variations might come from the same geophysical process.The WRMS residual reduction ratio for all the stations ranged from 19 to 85 %, which is better than Liu's work (Liu et al., 2014).However, the correlation and WRMS reduction between GPS and GRACE are weak when atmospheric and non-tidal ocean effects were removed, with the average correlation and WRMS reduction reduce to 75.6 and 28.9 % respectively.This is mainly because the seasonal hydrologic process is a major contributor to seasonal changes, and different stations are greatly influenced by the surrounding hydrological process.

Removing hydrological loading displacement from GPS using GRACE data
Coordinate variations measured by GPS stations, principally for the vertical component, have been used to investigate global (Dong et al., 2002) and local (Grapenthin et al., 2006) tectonic activity, as well as seasonal displacement modes for constraining estimates of continental, atmospheric, and ocean water storage.Some previous studies (e.g.Fu et al., 2012)   tion with relying on the accurate interpretation of GPS motion in terms of surface stress or tectonic movement.Thus, the displacement signal from surface mass loading is a source of noise (van Dam et al., 2007).For these applications, they would like to obtain reliable loading models or even surface mass observations, which can be used to reduce the environmental loading contributions to the GPS observations.In this study, we also attempt to separate tectonic and hydrological effects using GRACE-derived hydrological vertical rates.As mentioned in the Sect.3, the good seasonal correlation between GRACE and GPS signals indicates that the long-term uplifts revealed by GRACE detections are probably true and mixed in the GPS measurements.
Figure 9 shows results for individual GPS time series.Crustal subsidence or uplift due to vertical tectonic motion and TWS changes in the studied period are clearly evident in the vertical component shown in most of the GPS stations; in fact, the analysis in the five IGS GPS stations (BJFS, BJSH, JIXN, TAIN, and ZHNZ) suggests that the GPS vertical time series can be described by two different rates around 2010, due to a continuous decrease in TWS from 2004 to 2009; towards the end of 2009 the rate of this decrease slowed and rate started to rise since 2010 (please see Fig. 6).Thus, we divided it into two sub-periods when fitting GPS and GRACE trend for these five stations (Fig. 9a).GPS trend changes indicate an overall uplift for the whole region at the 0.04-1.47mm yr −1 level from 2004 to 2009, but the rate of change direction is inconsistent in different GPS stations at −0.94-2.55mm yr −1 level from 2010 to 2013 (Table 1).
In addition, the long-term trend rate is different in different areas from 2010 to 2013 (Fig. 9b).For example, the trend rate from GPS measurements shows the uplift in western NCP (Shanxi-SX, and some of Hebei-HE stations), but opposite trends in the central and eastern plain of NCP (Beijing-BJ, Tianjin-TJ, and some of Hebei-HE stations).The groundwater depletion, which occurs in the shallow unconfined aquifers in Piedmont Plain, leads to the loading uplift effect from mass loss.But groundwater depletion occurs in the deep confined aquifers in the central and eastern plain of NCP, which causes serious ground subsidence, rather than ground uplift caused by groundwater loss.
It is also possible that some GPS signals could be a result of loading from changes in the distribution of water stored in the surface and ground around the GPS surrounding region.To remove those contributions, Stokes coefficients output from the GRACE model were used to compute crustal motion at the NCP (Fig. 8a), and then transform the monthly results into daily resolution data using a spline interpolation.In addition, GRACE solutions are corrected for GIA while GPS ones are not.Here Stokes coefficients results from A et al. (2013) were used to remove GIA effects from GPS measurements, which is about 0.2-0.4mm yr −1 in the land areas of China and 0.28-0.33mm yr −1 in the NCP (please see the Fig. S5 and Table S4).
We compute the GRACE-derived long-term uplift for all continuous GPS sites used in this paper.The results indicate an overall uplift for the NCP region.Then we remove this TWS-induced uplift from GPS actual observed vertical rates to derive the corrected vertical velocities.Figure 10 divided the time into two sub-periods (2004-2009 and 2010-2013) to indicate an overall long-term trend before (grey arrow) and after (red arrow) removing hydrological loading displacement for the whole region.Secular displacement results between 2004 and 2009 show that loading displacement due to the TWS loss reduce the uplift rate of GPS to some extent, and groundwater exploitation was the main contributor to crustal uplift caused by TWS loss in the NCP (BJFS, BJSH, and JIXN in the Fig. 10a).However, studies indicate that groundwater withdrawal produces localized subsidence, which can be largely relative to tectonic displacement (Bawden et al., 2001).Therefore, in this study, more attention was paid to land subsidence due to groundwater loss.

Land subsidence in the central and eastern of NCP
Land subsidence has been commonly observed in the NCP, and has become the main factor that impacts regional sustainable economic and social development (Guo et al., 2015).Over the past years, the scope and magnitude of land subsidence has expanded.In this study, we used GPS sites to obtain time series of land subsidence evolution characteristics.The trend rates from GPS sites, after removing the rates from GRACE-derived long-term uplift and GIA effects, can be seen in Fig. 10b (the grey background areas in the dashed white box) to reflect the rate of land subsidence from 2010 to 2013, which is because the groundwater exploitation in the deep confined aquifers has a more serious impact on land subsidence (Guo et al., 2015).The results show that Tianjin becomes the most serious subsidence area, e.g. in the Tanggu and Hangu districts (TJBH), with an average subsidence of ∼ 14 mm yr −1 after 2010; in Wuqing district (TJWQ), recent subsidence-averaged ∼ 43 mm yr −1 .Because the Cangxian district (HECX) is close in proximity to the Jinghai region of Tijian, the sedimentation rate of ∼ 20 mm yr −1 can represent the subsidence trend of south-western Tianjin.However, the difference of spatial distribution of land subsidence is large in Tianjin, and uneven settlement characteristics are obvious.For example, the south-western and western areas of Tianjin are the most serious areas, and the trend of land subsidence exists in the northward but the amplitude is small relative to the south-western and western areas, i.e.JIXN site shows a small negative trend (∼ −0.6 mm yr −1 ).The cause for subsidence in the Tianjin area is linked to over exploitation of groundwater, an issue that has not been effectively controlled resulting in rapidly developing land subsidence in the suburbs in recent years (Yi et al., 2011).
In the central and eastern region of NCP, where disastrous land subsidence has also occurred in Beijing and cities in central Hebei province and the north-east of Shandong province, for instance, large subsidence zone in Hebei province, has formed from north to south, which start from the western region of Beijing (BJFS, BJSH, and BJYQ station), via the eastern region of Xingtai and Handan (HELY station), extend to northern Hebi (HAHB station belong to Henan province).
However, results from our investigation show that the centre of land subsidence does not completely overlap the TWS loss contributions (see the secular trend maps of the TWS changes of NCP in Fig. 2).The uplift still exists even when we removed the rates from GRACE-derived and GIA effects in the piedmont of Taihang Mountains and the western part of NCP (Shanxi province), where the groundwater depletion occurs in the shallow unconfined aquifers have not led to a large area of subsidence.The reason for this difference with the western region of NCP is that crustal uplift is mainly controlled by tectonic movement, which is the orogenic belt and plateau area in the western Taihang Mountains, which is basic in the uplift.In our results, most of the corrected vertical velocities at GPS stations, especially in the central and eastern region of NCP, agree with the previous study results, i.e. combining with mobile and continuous GPS observation (Zhao et al., 2014) and using GPS stations from GNSS and levelling data (MLR, 2015).For results of vertical crust movement in the NCP from the previous study, please see the Fig. S6.

Conclusions
Temporal variations in the geographic distribution of surface mass (continental water, ocean mass, and atmospheric mass) can lead to displacement of the Earth's surface.Due to excessive exploitation of groundwater resources the NCP area has become susceptible to land subsidence, and it has become one of the most affected areas in the world.Calculating the loading displacement can explain the natural displacement phenomenon, and it presents new insight into the dynamics of land subsidence.
Traditional displacement observation has space limitations.Based on the elastic displacement of the Earth's crust by surface loadings, this study combined GRACE and GPS data to investigate vertical displacements in the NCP area.GRACE data were used to model vertical displacements due to changes in hydrological loads.The results showed both GPS and GRACE data to observe strong seasonal variations.Comparisons between the observed GPS seasonal vertical displacement and GRACE-derived seasonal displacement demonstrated that a consistent physical mechanism is responsible for TWS changes; i.e. the seasonal hydrospheric Hydrol.Earth Syst.Sci., 21, 2905Sci., 21, -2922Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/2905/2017/ mass movements due to climate variability cause periodic displacements of the lithosphere.
As well as the significant seasonal characteristics, GRACE also exhibited a long-term mass loss in this region; the rate of GRACE-derived TWS loss (after multiplying by the scaling factor) in the NCP was 3.39 cm yr −1 from 2003 to 2009, which is equivalent to a volume of 12.42 km 3 yr −1 .The rate was 2.57 cm yr −1 from 2003 to 2012, equivalent to a volume of 9.41 km 3 yr −1 .The TWS loss was principally due to groundwater depletion in the NCP.We calculated that the consequent trend rate caused by the load mass change using GRACE data and removed this hydrological effect from observed GPS vertical rates.Secular displacement results showed that TWS losses reduced loading displacement to some extent, but the trend rates disagree due to the difference of spatial distribution with anthropogenic depletion of TWS in the NCP.
Particularly, land subsidence has been affecting the central and eastern region of NCP, especially in Tianjin, for the past years.Overpumping of groundwater is the main cause of land subsidence, which has led to comprehensive detrimental effects on the society, the economy, and the natural environment.The impact of groundwater exploitation in different aquifer systems and active faults in the different regions on land subsidence needs to be analysed in future investigations.For example, using GRACE to remove mass-loading signals from a GPS record requires either confidence that there is no concentrated load signal very near the site, or a scaling factor based on a reliable model of the mass change (the groundwater depletion rate estimated from monitoring well stations) pattern around the site.
Data availability.The GPS data of CMONOC and IGS were made by First Crust Monitoring and Application Center, China Earth-quake Administration (http://www.eqdsc.com/data/pgv-sjxl.htm).The groundwater level data series were obtained from MWR (the website of data integration and services in the groundwater resources is available at http://www.groundwater.cn/default.html),and the daily precipitation data for weather stations were collected from CMDSSS (http://data.cma.cn/site/index.html).The GRACE solutions are available at ftp://podaac.jpl.nasa.gov/allData/grace/L2/CSR/RL05/ and the GLDAS/Noah model data provided by the NASA Goddard Earth Sciences Data and Information Services Center (http://disc.sci.gsfc.nasa.gov/).

Figure 1 .
Figure 1.Study region of the North Plain China (NCP) showing locations of continuous GPS stations.Blue dots represent continuous GPS sites in the Crustal Movement Observation Network of China (CMONOC) and red stars represent the International GNSS Service (IGS) sites).Cities and provinces are labelled as follows: Beijing (BJ), Tianjin (TJ), Hebei province (HE), and Shanxi province (SX).

Figure 2 .
Figure 2. The 2003-2012 secular trend maps (cm yr −1 ) of the terrestrial water storage (TWS) changes in the North Plain China (NCP) and surrounding regions derived from GRACE data.Results have been destriped and smoothed with a 250 km Gaussian smoothing function.

Figure 3 .
Figure 3. (a) Daily values of the vertical (positive upward) components of position, as measured at IGS GPS sites BJFS, ZHNZ, BJSH, JIXN, and TAIN.The example of displacement due to atmospheric and non-tidal ocean loading at BJFS IGS sites are shown in (b).

Figure 4 .
Figure 4. Surface horizontal (north and east components) and vertical deformation modelled by GRACE in four IGS sites.(a) and (b) show the time series and trend rates of north and east components in BJFS, JIXN, TAIN, and ZHNZ respectively, (c) shows the time series of vertical displacements.

Figure 5 .
Figure 5.Time series showing daily values (a) and fitting results (b) of the vertical (positive upward) components from GPS and GRACEderived at five IGS GPS sites.

Figure 6 .
Figure 6.Comparison of annual amplitudes and initial phases between GPS (grey) and GRACE (red).The initial phases are anticlockwise from the east (reference time is 2004.0).

Figure 7 .
Figure 7. Time series showing total terrestrial water storage (TWS) changes in the spatially averaged area (kernel) of the NCP estimated from CSR GRACE data, monthly groundwater level changes observed by monitoring wells after multiplying by mean value of specific yields in the NCP during 2002-2012 and the daily precipitation data (rainfall amount) for weather stations during the period of 2003-2012 from CMDSSS.The black dashed curve is the temporal smoothing GRACE-based result, the red and blue dashed curve are the long trend of GRACE-based result during 2003-2009 and 2003-2012 respectively.

Figure 8 .
Figure 8. GRACE-derived smoothed (dash curves) and long-term (solid curves) vertical displacement time series due to load changes (a), the groundwater depletion contributions estimated from GRACE minus GLDAS data for smoothed (dash curves) and long-term (solid curves) vertical displacements (b), as measured at five IGS stations and 20 CMONOC stations in NCP and its surrounding region.The grey background shows inflexion effects due to TWS changes in the vertical component.

Figure 9 .
Figure 9. Smoothed (solid curves) and long-term (dash curves) versions of daily values of the vertical (positive upward) component of position, as measured at 29 GPS sites in NCP and its surrounding region, (a) five IGS stations and (b) 24 CMONOC stations.

Figure 10 .
Figure 10.GPS (grey arrow, positive upward) and corrected GPS (red arrow, positive upward) vertical trend rate after subtracting the GRACE-derived long-term uplift rate due to load changes and GIA effect between 2004 and 2009 (a), and between 2010 and 2013 (b).

Table 1 .
GPS station information.
1 IGS sites: the observation time between 2003 and 2013. 2 CMONOC sites: the observation time between 2010 and 2013.

Table 2 .
Trends of GRACE-derived TWS and GWS (GRACE minus GLDAS/Noah) in situ measurements (shallow aquifers) and compared with the previous studies during 2003-2012.