Vertical Displacements Driven by Groundwater Storage Changes in the North China Plain Detected by GPS Observations

The North China Plain (NCP) has been experiencing the most severe groundwater depletion in China, leading to a broad region of vertical motions of the Earth’s surface. This paper explores the seasonal and linear trend variations of surface vertical displacements caused by the groundwater changes in NCP from 2009 to 2013 using Global Positioning System (GPS) and Gravity Recovery and Climate Experiment (GRACE) techniques. Results show that the peak-to-peak amplitude of GPS-derived annual variation is about 3.7~6.0 mm and is highly correlated (R > 0.6 for most selected GPS stations) with results from GRACE, which would confirm that the vertical displacements of continuous GPS (CGPS) stations are mainly caused by groundwater storage (GWS) changes in NCP, since GWS is the dominant component of total water storage (TWS) anomalies in this area. The linear trends of selected bedrock-located IGS CGPS stations reveal the distinct GWS changes in period of 2009–2010 (decrease) and 2011–2013 (rebound), which are consistent with results from GRACE-derived GWS anomalies and in situ GWS observations. This result implies that the rate of groundwater depletion in NCP has slowed in recent years. The impacts of geological condition (bedrock or sediment) of CGPS stations to their results are also investigated in this study. Contrasted with the slight linear rates (−0.69~1.5 mm/a) of bedrock-located CGPS stations, the linear rates of sediment-located CGPS stations are between −44 mm/a and −17 mm/a. It is due to the opposite vertical displacements induced by the Earth surface’s porous and elastic response to groundwater depletion. Besides, the distinct renewal characteristics of shallow and deep groundwater in NCP are discussed. The GPS-based vertical displacement time series, to some extent, can reflect the quicker recovery of shallow unconfined groundwater than the deep confined groundwater in NCP; through one month earlier to attain the maximum height for CGPS stations nearby shallow groundwater depression cones than those nearby deep groundwater depression cones.


Introduction
Groundwater resources are essential to urban and agricultural needs in North China Plain (NCP), where the precipitation is low, and more than 70-80% of fresh water comes from groundwater [1], and agricultural irrigation almost consumes 70% of total water supply [2].Historical exploitation of groundwater has caused a remarkable decline in groundwater storage (GWS) in NCP, which is referred to as groundwater depletion and has been verified by many previous studies [3][4][5][6].Large-scale GWS changes exert loading and unloading of the Earth's surface, inducing vertically downward and upward displacements of the Earth's crust and uppermost mantle [7][8][9].It should be noted that the ground surface subsidence induced by mining excavation is not discussed by this paper [10,11].
With the advances in space geodetic sensors, modern space geodetic measurements can quantify both total water storage (TWS) variations with the Gravity Recovery and Climate Experiment (GRACE) and Earth surface displacements with the Global Positioning System (GPS) [12][13][14][15].Previous studies have provided significant insights into the seasonal and long-term vertical displacements in NCP by using GPS and GRACE data.For instance, Liu et al. (2014) [16] found that the GPS-and GRACE-derived secular trends and seasonal signals are in good agreement, respectively, with an uplift magnitude of 1-2 mm/year and a correlation of 85.0-98.5%;Wang et al. (2017) [17] identified that the average correlation and weighted root-mean-squares (WRMS) reduction between GPS and GRACE are 75.6% and 28.9%.These previous studies have focused on comparison, cross-validation of GPS, and GRACE results.
According to the distinct hydrogeological settings within the quaternary sediments [18,19], the plain area of NCP can be divided into the piedmont plain (PP) and the east central plain (ECP, alluvial and coastal plains) (Figure 1), where water supply mainly depends on shallow fresh groundwater stored in the shallow unconfined aquifers and deep confined aquifers, respectively.However, the shallow groundwater in PP and deep groundwater in ECP with different characteristics over the heterogeneous aquifer systems have not been studied either from the perspective of GPS displacement monitoring field or from the view of GRACE hydrological applications; this study attempts to explore the different influences of shallow groundwater and deep groundwater on GPS results.
Remote Sens. 2018, 10, x FOR PEER REVIEW 2 of 13 Historical exploitation of groundwater has caused a remarkable decline in groundwater storage (GWS) in NCP, which is referred to as groundwater depletion and has been verified by many previous studies [3][4][5][6].Large-scale GWS changes exert loading and unloading of the Earth's surface, inducing vertically downward and upward displacements of the Earth's crust and uppermost mantle [7][8][9].It should be noted that the ground surface subsidence induced by mining excavation is not discussed by this paper [10,11].With the advances in space geodetic sensors, modern space geodetic measurements can quantify both total water storage (TWS) variations with the Gravity Recovery and Climate Experiment (GRACE) and Earth surface displacements with the Global Positioning System (GPS) [12][13][14][15].Previous studies have provided significant insights into the seasonal and long-term vertical displacements in NCP by using GPS and GRACE data.For instance, Liu et al. (2014) [16] found that the GPS-and GRACE-derived secular trends and seasonal signals are in good agreement, respectively, with an uplift magnitude of 1-2 mm/year and a correlation of 85.0-98.5%;Wang et al. (2017) [17] identified that the average correlation and weighted root-mean-squares (WRMS) reduction between GPS and GRACE are 75.6% and 28.9%.These previous studies have focused on comparison, cross-validation of GPS, and GRACE results.
According to the distinct hydrogeological settings within the quaternary sediments [18,19], the plain area of NCP can be divided into the piedmont plain (PP) and the east central plain (ECP, alluvial and coastal plains) (Figure 1), where water supply mainly depends on shallow fresh groundwater stored in the shallow unconfined aquifers and deep confined aquifers, respectively.However, the shallow groundwater in PP and deep groundwater in ECP with different characteristics over the heterogeneous aquifer systems have not been studied either from the perspective of GPS displacement monitoring field or from the view of GRACE hydrological applications; this study attempts to explore the different influences of shallow groundwater and deep groundwater on GPS results.GRACE grids are drawn with oblique lines, following the regional boundary.All the calculations (e.g., GRACE-derived TWS and GWS, GLDAS-modeled SWS, and precipitation) related to the region are based on the extents of GRACE grids.The background is the region's digital elevation model (DEM) information.
In this paper, in contrast to previous studies, except for seasonal and long-term vertical displacement due to dynamic hydrological processes and groundwater depletion, the capability of GPS to detect the different impacts of shallow and deep groundwater in NCP, as well as the influence of distinct the geological conditions (bedrock or sediment) of continuous GPS (CGPS) stations, are also deeply investigated and discussed.In this study, the assessment will be based on GPS data from the Crustal Movement Observation Network of China (CMONOC, Figure 1 shows the locations of selected CGPS stations), level-2 Release-05 (RL05) GRACE data from the University of Texas Center for Space Research (UTCSR), Global Land Data Assimilation System (GLDAS) hydrologic model simulations, and in situ groundwater-level observations of both shallow unconfined and deep confined aquifers.

GPS Data
GPS can measure the crustal deformation at millimeter-level accuracy, due to the recent improvements of GPS processing strategies [20,21].Ten former and current CGPS stations in NCP (Figure 1) are analyzed in this study.Three of them (BJFS, BJSH, JIXN) are from the International GNSS Service (IGS), and others are from CMONOC [22].All the CGPS stations are equipped with multi-frequency receivers, with a sampling rate of 30 s.The specific information of selected CGPS stations are listed in Table 1, including their topographic features (PP or ECP) and basement type (bedrock or sediment).In this study, the rinex files of the three IGS stations, from 2009 to 2013, are used, and that of other CGPS stations are from 2010 to 2013, since these stations began to collect data from 2010.GAMIT software 10.4 is employed in "baseline mode" to process the GPS data [23], to obtain GPS daily, loosely constrained coordinates, and covariance.In this procedure, Solid-Earth tides and pole tides are corrected by models recommended by the International Earth Rotation Service (IERS) [24].In addition, Finite Element Solutions 2004 (FES2004) ocean tidal loading corrections and absolute phase-center corrections for satellites and antenna, as issued by International GNSS Service (IGS) [25], are used to correct corresponding errors.The tropospheric delays are modeled using the Global Mapping Function (GMF) [26] every two hours and the first order ionospheric delays are eliminated by means of the ionosphere-free linear combination.Gross error rejection is implemented to the baseline resolutions, and a sequential Kalman filter method is employed to combine regional, loosely constrained solutions with global solutions obtained from SOPAC (http://sopac.ucsd.edu)using GLOBK software 10.4.Finally, the datum transformation is performed to obtain a station coordinate time series with respect to the ITRF2008 using 45 global distributed IGS stations [27].According to the above GPS data processing strategy, the daily vertical coordinate time series of the selected CGPS stations is obtained.

TWS and GWS Variations Derived from GRACE
Spherical harmonic coefficients of the Earth's gravity field estimated from GRACE data [28,29] and load Love numbers can be used to model the TWS variations [30] and vertical surface displacements [31] caused by changing water mass loading.The monthly GRACE Level-2 RL05 solution GSM gravity data product, in the form of spherical harmonic coefficients (SHCs) provided from the UTCSR, whose period is from January 2009 to December 2013, up to degree 96, are used.C20 terms are replaced with the results from observations of Satellite Laser Ranging [32].Due to the fact that the reference frame origin used in the GRACE gravity field determination is the Earth's center of mass, GRACE gravity field cannot determine the degree-1 terms variations (geocenter motion), whereas GPS data contain degree-1 terms; therefore, in order to be consistently comparable to the GPS time series, we replace the degree-1 components in GRACE with the results derived by Swenson et al. (2008) [33].
First, the average mean gravity field model between 2009 and 2013 is removed from the data product.Then, the destriping method P3M6 of Chen et al. (2007) [34] and Gaussian smoothing with averaging radius of 300 KM are employed to the SHCs of the gravity fields, in the purpose of reducing the north-south stripes and errors at high degrees.Finally, the TWS variations, in the form of equivalent water height, and vertical surface displacements due to changes in mass loading, can be, respectively, expressed in terms of residual SHCs of gravity field and load Love numbers as follows [30,31]: in which ∆H(θ, λ) and dr(θ, λ) are the variation of equivalent water height and surface displacement in the radial direction at geocentric colatitude θ and geocentric longitude λ, respectively; R 0 and ρ E are the mean radius and density of the Earth, respectively; ρ w is water density; P lm is fully normalized Legendre function of degree l and order m; ∆C lm and ∆S lm represent the residual spherical harmonic coefficients of the destriped and smoothed gravity field from which the average gravity field between 2009 and 2013 have been removed; k l and h l are Load Love numbers at degree l.Here, the load Love numbers from Farrell [7] are adopted.
Since GRACE observations at small scale regions tend to be attenuated in the post-processing of smoothing, the modeled TWS should be adjusted by scaling factors.Four land surface models (LSMs: Community Land Model, Mosaic, Noah, and Variable Infiltration Capacity model) in GLDAS version 1 [35] are used to calculate scaling factors, by calculating the ratio between the original true TWS and the truncated and filtered TWS [36].In this study, the average scaling factor in NCP is about 1.22.
At monthly or longer time scales, soil water storage (SWS) and GWS changes account for the majority of TWS changes in NCP; other components such as the snow water equivalent are negligible [4].Therefore, GWS anomalies are obtained by subtracting the SWS anomalies simulated by hydrologic models from the adjusted GRACE TWS anomalies.The SWS anomalies are also estimated by using the average of simulated soil moisture data (2009-2013) from the four LSMs in GLDAS.

In Situ Groundwater LEVEL MEASUREMENTS
The in situ monthly groundwater level measurements are obtained from daily observations from 107 monitoring well stations in NCP from 2009 to 2013, including 47 shallow unconfined wells and 60 deep confined wells (provided by Institute of China Geological Environment Monitoring).The groundwater data adopted in this study can provide more spatial details of shallow and deep GWS variations compared with that used in previous similar study of Wang et al. (2017) [17], which contained the monthly groundwater table depth of 20 cities, most of them from a shallow unconfined well.As shown in Figure 1, shallow unconfined wells are mainly distributed in PP, while deep confined wells are mainly located in ECP.We obtain in situ GWS variations by multiplying in situ-measured groundwater level anomalies with the specific yield of each well as obtained from the spatial distribution map of specific yields in NCP provided by the China Geological Survey Bureau [37].The distribution of the specific yield in NCP ranges widely from <0.05 in the coastal plain around the margin of the Bohai Sea to >0.2 in the piedmont region of the Taihang Mountains.Based on the isopleth map of specific yields in the NCP and the locations of monitoring well, we conclude that a reasonable range of average values for specific yields was 0.05~0.07,which we use to estimate GWS variations.

Seasonal and Secular Signals Revealed by GPS and GRACE
The GPS-measured and GRACE-modeled [31] vertical displacements are fitted with an empirical model including epoch position, speed, and the amplitude and phase of annual and semiannual harmonic components (to model seasonal effects) using least-squares method [38] over the time span from 2009 to 2013.Moreover, the offsets caused by identified equipment changes and earthquakes are removed in GPS analysis, and only stations with, at least, two-and-a-half years of data are included.As shown in Figure 2, both GPS and GRACE results show apparent seasonal variations.GPS and GRACE record coherent vertical crustal seasonal variations from 2009 to 2013 with a high correlation for most selected CGPS stations (R > 0.6, except for BJSH and TJBH) (Table 2).For BJSH and TJBH stations, it can be inferred from Figure 2 (left) that the low correlations are mainly caused by the disagreements in the phases of GPS-and GRACE-derived results.These inconsistences of phases reflect the different local sensitivity of the two techniques.The good correlation between these two geodetic measurements confirms that the seasonal crustal deformation in NCP is mainly caused by the seasonal hydrological mass loading variations, which can be resolved by geodetic measurements of both GPS and GRACE.For all the selected CGPS stations, the annual component is more dominant than the semi-annual one (Figure 2).The peak-to-peak amplitude of GPS-and GRACE-derived annual variation is about 3.7~6.0mm and 2.5~3.3 mm, respectively (Table 2), which is about two and three times of semi-annual amplitude.Therefore, the following discussions are focused on annual variations.
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 13 contained the monthly groundwater table depth of 20 cities, most of them from a shallow unconfined well.As shown in Figure 1, shallow unconfined wells are mainly distributed in PP, while deep confined wells are mainly located in ECP.We obtain in situ GWS variations by multiplying in situ-measured groundwater level anomalies with the specific yield of each well as obtained from the spatial distribution map of specific yields in NCP provided by the China Geological Survey Bureau [37].The distribution of the specific yield in NCP ranges widely from ˂0.05 in the coastal plain around the margin of the Bohai Sea to ˃0.2 in the piedmont region of the Taihang Mountains.Based on the isopleth map of specific yields in the NCP and the locations of monitoring well, we conclude that a reasonable range of average values for specific yields was 0.05~0.07,which we use to estimate GWS variations.

Seasonal and Secular Signals Revealed by GPS and GRACE
The GPS-measured and GRACE-modeled [31] vertical displacements are fitted with an empirical model including epoch position, speed, and the amplitude and phase of annual and semiannual harmonic components (to model seasonal effects) using least-squares method [38] over the time span from 2009 to 2013.Moreover, the offsets caused by identified equipment changes and earthquakes are removed in GPS analysis, and only stations with, at least, two-and-a-half years of data are included.As shown in Figure 2, both GPS and GRACE results show apparent seasonal variations.GPS and GRACE record coherent vertical crustal seasonal variations from 2009 to 2013 with a high correlation for most selected CGPS stations (R > 0.6, except for BJSH and TJBH) (Table 2).For BJSH and TJBH stations, it can be inferred from Figure 2 (left) that the low correlations are mainly caused by the disagreements in the phases of GPS-and GRACE-derived results.These inconsistences of phases reflect the different local sensitivity of the two techniques.The good correlation between these two geodetic measurements confirms that the seasonal crustal deformation in NCP is mainly caused by the seasonal hydrological mass loading variations, which can be resolved by geodetic measurements of both GPS and GRACE.For all the selected CGPS stations, the annual component is more dominant than the semi-annual one (Figure 2).The peak-to-peak amplitude of GPS-and GRACE-derived annual variation is about 3.7~6.0mm and 2.5~3.3 mm, respectively (Table 2), which is about two and three times of semi-annual amplitude.Therefore, the following discussions are focused on annual variations.Table 2 shows the peak-to-peak amplitudes and initial phases of annual vertical displacements derived by GPS and GRACE.The amplitude of annual variation derived from GPS (3.7~6.0 mm) is slightly higher than that from GRACE (2.5~3.3 mm) for all the selected CGPS stations.This is because GPS observes all the earth surface deformation induced by tectonic and non-tectonic Table 2 shows the peak-to-peak amplitudes and initial phases of annual vertical displacements derived by GPS and GRACE.The amplitude of annual variation derived from GPS (3.7~6.0 mm) is slightly higher than that from GRACE (2.5~3.3 mm) for all the selected CGPS stations.This is because GPS observes all the earth surface deformation induced by tectonic and non-tectonic (hydrology, atmosphere, ocean et al.) factors, while GRACE just models the hydrological component.Compared to GPS solutions, the spatial coherence of phases from GRACE solutions is more apparent, which indicates that GPS has a strong sensitivity for local surface loading, while the spatial resolution for GRACE is several hundred kilometers.This is also the reason for the small difference of GRACE-derived vertical displacements among selected stations.Besides, the fitting results in Table 2 show that the linear trend rate of GRACE-derived vertical displacements for the whole region is 0.5~0.7 mm/a from 2009 to 2013, which reveals the TWS loss for the entire area of NCP in observation time span.However, the rate of change direction is inconsistent for different GPS stations ranging from −44.0 mm/a to 1.5 mm/a, which should be mainly due to Earth surface's elastic and porous response to hydrological mass variations and is deeply discussed in Section 3.2.As mentioned above, the good seasonal correlation between GPS and GRACE signals indicates that the long-term uplifts revealed by GRACE solutions are probably true and mixed in the GPS measurements.Figure 3  (hydrology, atmosphere, ocean et al.) factors, while GRACE just models the hydrological component.Compared to GPS solutions, the spatial coherence of phases from GRACE solutions is more apparent, which indicates that GPS has a strong sensitivity for local surface loading, while the spatial resolution for GRACE is several hundred kilometers.This is also the reason for the small difference of GRACE-derived vertical displacements among selected stations.Besides, the fitting results in Table 2 show that the linear trend rate of GRACE-derived vertical displacements for the whole region is 0.5~0.7 mm/a from 2009 to 2013, which reveals the TWS loss for the entire area of NCP in observation time span.However, the rate of change direction is inconsistent for different GPS stations ranging from −44.0 mm/a to 1.5 mm/a, which should be mainly due to Earth surface's elastic and porous response to hydrological mass variations and is deeply discussed in Section 3.2.As mentioned above, the good seasonal correlation between GPS and GRACE signals indicates that the long-term uplifts revealed by GRACE solutions are probably true and mixed in the GPS measurements.Figure 3 shows the daily observations of three bedrock-located IGS CGPS stations (BJFS, BJSH, JIXN), since other CGPS stations' data are only available after July 2010.The long-term linear trends of crustal uplift or subsidence are mainly dominated by surface mass loading and tectonic processes.The analysis of the three IGS CGPS stations (BJFS, BJSH, JIXN) suggests the slight change of linear trend around 2010.From 2009 to 2010, the linear trend of vertical displacements is uplift, which is consistent with the reported results of Wang et al. (2017) [17] and might imply the unloading of hydrological mass in NCP.However, the trend changes to lightly downward or relatively stable after the end of 2010, indicating the increase of hydrological mass loading in this area.These facts reveal that the severe water condition in NCP might has been alleviated in recent years.In order to verify and explain the linear trend change of three IGS CGPS results, the GRACE-derived and in situ GWS measurements are also explored.Focusing on Figure 4a, compared to GWS variation, SWS component contributes quite insignificantly to GRACE-derived TWS changes, proving that the total hydrological mass loss in NCP is mainly caused by depletion of GWS.The GRACE-derived GWS monthly anomalies in NCP are correlated with in situ-measured GWS anomalies at a coefficient of R ≈ 63 from 2009 to 2013.However, the relatively high correlation is mainly due to the agreement in the phase of variation.As presented in Figure 4a, from 2009 to 2013, both the GRACE and in situ results show a distinct feature with two distinguishable stages: Stage 1: from 2009 to 2010, with continuously decreasing trend at a rate of about −45 mm/a for both GRACE-derived GWS results and in situ-GWS observations, leading to upward vertical movement of CGPS stations (bedrock); Stage 2, from 2011 to 2013, with slightly rebounded and relatively stable for both GRACE-derived and in situ results, resulting in downward or stable vertical movement of CGPS stations (bedrock).This difference between the two sub-periods (2009-2010 and 2011-2013) could be mainly induced by climate-driven precipitation annual variations in NCP (Figure 4c) and partially owing to the impacts of the implementation of the South-to-North Water Diversion (SNWD) project since 2003, which reflect the influence of the SNWD project on reducing GWS depletion in NCP.Some previous literatures also reported positive impacts on GWS recovery in NCP [39].
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 13 In order to verify and explain the linear trend change of three IGS CGPS results, the GRACE-derived and in situ GWS measurements are also explored.Focusing on Figure 4a, compared to GWS variation, SWS component contributes quite insignificantly to GRACE-derived TWS changes, proving that the total hydrological mass loss in NCP is mainly caused by depletion of GWS.The GRACE-derived GWS monthly anomalies in NCP are correlated with in situ-measured GWS anomalies at a coefficient of R ≈ 63 from 2009 to 2013.However, the relatively high correlation is mainly due to the agreement in the phase of variation.As presented in Figure 4a, from 2009 to 2013, both the GRACE and in situ results show a distinct feature with two distinguishable stages: Stage 1: from 2009 to 2010, with continuously decreasing trend at a rate of about −45mm/a for both GRACE-derived GWS results and in situ-GWS observations, leading to upward vertical movement of CGPS stations (bedrock); Stage 2, from 2011 to 2013, with slightly rebounded and relatively stable for both GRACE-derived and in situ results, resulting in downward or stable vertical movement of CGPS stations (bedrock).This difference between the two sub-periods (2009-2010 and 2011-2013) could be mainly induced by climate-driven precipitation annual variations in NCP (Figure 4c) and partially owing to the impacts of the implementation of the South-to-North Water Diversion (SNWD) project since 2003, which reflect the influence of the SNWD project on reducing GWS depletion in NCP.Some previous literatures also reported positive impacts on GWS recovery in NCP [39].

Impacts of Bedrock and Sediment Basement to GPS Results
As presented in Table 1, except for HAHB, HECX, TJBH, and TJWQ stations, which are located on soft sediment and mainly distributed over the ECP area, the other CGPS stations are located on bedrock throughout the PP area and north part of ECP.Their monthly vertical displacement variations are displayed in Figure 5.For most bedrock-located CGPS stations (red bars, except for BJSH), whose vertical displacements mainly reflect elastic response of solid Earth's surface to the loading of hydrological masses in this area, the maximum height occurs around early summer (June), and a minimum height takes place around early winter(December).This temporal pattern of vertical displacement should be mainly determined by the composite impact from agricultural groundwater extraction and climate factors.On one hand, over the whole region of NCP, especially the PP area, water-consumptive winter wheat is widely cultivated and agriculture is highly supported by artificial irrigation, which contribute greatly to groundwater depletion, as Hu et al. ( 2016) [40] reported that the correlation coefficients between irrigation amount and groundwater level anomalies is about −0.7.In addition, Tan et al. (2012) [41] showed that the maximum agricultural irrigation is concentrated in February to May in this region.On the other hand, precipitation is mainly focused on July and August (Figure 4c), which is stored in aquifers for next spring due to the lower irrigation requirement in those months.Hence, the vertical displacement due to groundwater loading displays a clear seasonal pattern of upward trend from January to June, when the groundwater is continuously extracted for irrigation, and downward trend from July to December, when the groundwater is recharged and accumulated from precipitation.

Impacts of Bedrock and Sediment Basement to GPS Results
As presented in Table 1, except for HAHB, HECX, TJBH, and TJWQ stations, which are located on soft sediment and mainly distributed over the ECP area, the other CGPS stations are located on bedrock throughout the PP area and north part of ECP.Their monthly vertical displacement variations are displayed in Figure 5.For most bedrock-located CGPS stations (red bars, except for BJSH), whose vertical displacements mainly reflect elastic response of solid Earth's surface to the loading of hydrological masses in this area, the maximum height occurs around early summer (June), and a minimum height takes place around early winter(December).This temporal pattern of vertical displacement should be mainly determined by the composite impact from agricultural groundwater extraction and climate factors.On one hand, over the whole region of NCP, especially the PP area, water-consumptive winter wheat is widely cultivated and agriculture is highly supported by artificial irrigation, which contribute greatly to groundwater depletion, as Hu et al. ( 2016) [40] reported that the correlation coefficients between irrigation amount and groundwater level anomalies is about −0.7.In addition, Tan et al. (2012) [41] showed that the maximum agricultural irrigation is concentrated in February to May in this region.On the other hand, precipitation is mainly focused on July and August (Figure 4c), which is stored in aquifers for next spring due to the lower irrigation requirement in those months.Hence, the vertical displacement due to groundwater loading displays a clear seasonal pattern of upward trend from January to June, when the groundwater is continuously extracted for irrigation, and downward trend from July to December, when the groundwater is recharged and accumulated from precipitation.By contrast, sediment-located CGPS stations (blue bars) respond differently.Those stations attain their maximum height from April to May and minimum height from October to November, about one to two months earlier compared to bedrock-located CGPS stations, which are likely to be attributed to the conjunctive impacts of elastic and porous response of Earth's surface to groundwater variation.Except for elastic response to groundwater loading and unloading, due to soil pore compression, the Earth's surface rises when the aquifer is recharged with water and subsides when the groundwater is withdrawn from the aquifer [42], which is opposite to Earth's elastic response to the loading of groundwater and could obfuscate the effects of hydrological mass variation on position [43,44].For the CGPS stations located on soft sediment, their vertical positions record not only the Earth's elastic deformation to groundwater mass loading and unloading, but also the porous response associated with groundwater filling and emptying the aquifers.Thus, for these stations, their vertical displacement variations reflect both elastic and porous responses of Earth's surface to groundwater changes, while bedrock-located CGPS stations mainly reveal the elastic response.In our study, the five of six bedrock-located CGPS stations in NCP attain their maximum heights at the time of early summer when groundwater is minimum due to large By contrast, sediment-located CGPS stations (blue bars) respond differently.Those stations attain their maximum height from April to May and minimum height from October to November, about one to two months earlier compared to bedrock-located CGPS stations, which are likely to be attributed to the conjunctive impacts of elastic and porous response of Earth's surface to groundwater variation.Except for elastic response to groundwater loading and unloading, due to soil pore compression, the Earth's surface rises when the aquifer is recharged with water and subsides when the groundwater is withdrawn from the aquifer [42], which is opposite to Earth's elastic response to the loading of groundwater and could obfuscate the effects of hydrological mass variation on position [43,44].For the CGPS stations located on soft sediment, their vertical positions record not only the Earth's elastic deformation to groundwater mass loading and unloading, but also the porous response associated with groundwater filling and emptying the aquifers.Thus, for these stations, their vertical displacement variations reflect both elastic and porous responses of Earth's surface to groundwater changes, Remote Sens. 2018, 10, 259 9 of 13 while bedrock-located CGPS stations mainly reveal the elastic response.In our study, the five of six bedrock-located CGPS stations in NCP attain their maximum heights at the time of early summer when groundwater is minimum due to large pumping for agricultural usage and low recharge from rainfall, while all the sediment-located CGPS stations attain their maximum heights in late spring, due to opposite influence from Earth's surface elastic and porous response to groundwater changes.
Figure 6 shows distribution of groundwater depression cones across NCP.Shallow groundwater depression cones are mostly distributed in the PP area near Shijiazhuang, while the deep groundwater depression cones are primarily located in ECP, near Cangzhou, Hengshui, and Tianjin.The linear trend of vertical displacements at individual CGPS stations is also shown in Figure 6, the uplift or subsidence of which is distinct in different areas.Most selected bedrock-located CGPS stations in PP and north part of ECP show slight secular trend of −0.69~1.5 mm/a from 2009/2010 to 2013.However, sediment-located CGPS stations in ECP (HECX, TJBH, and TJWQ) show significant subsidence vertical rate at −44~−17 mm/a.The three CGPS stations are located on soft sediments and along the margins of deep groundwater depression cones, near Tianjin or Cangzhou, which have been suffering the most severe groundwater depletion in these areas [45] (Hu et al., 2004).For these stations, as discussed above, their GPS vertical position record both earth's elastic deformation associated with changes in GWS and pore compaction of groundwater depletion.Hence, it could be concluded that, for the sediment-located stations nearby deep groundwater depression cones, the subsidence caused by pore compaction are much larger than the elastic uplift induced by groundwater depletion.This indicates that the groundwater depletion in deep confined aquifers is the major triggering factor for inducing land subsidence in middle-east plain of NCP, which is consistent with the reported results of Guo et al. 2015 [46].
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 13 pumping for agricultural usage and low recharge from rainfall, while all the sediment-located CGPS stations attain their maximum heights in late spring, due to opposite influence from Earth's surface elastic and porous response to groundwater changes.Figure 6 shows distribution of groundwater depression cones across NCP.Shallow groundwater depression cones are mostly distributed in the PP area near Shijiazhuang, while the deep groundwater depression cones are primarily located in ECP, near Cangzhou, Hengshui, and Tianjin.The linear trend of vertical displacements at individual CGPS stations is also shown in Figure 6, the uplift or subsidence of which is distinct in different areas.Most selected bedrock-located CGPS stations in PP and north part of ECP show slight secular trend of −0.69~1.5 mm/a from 2009/2010 to 2013.However, sediment-located CGPS stations in ECP (HECX, TJBH, and TJWQ) show significant subsidence vertical rate at −44~−17 mm/a.The three CGPS stations are located on soft sediments and along the margins of deep groundwater depression cones, near Tianjin or Cangzhou, which have been suffering the most severe groundwater depletion in these areas [45] (Hu et al., 2004).For these stations, as discussed above, their GPS vertical position record both earth's elastic deformation associated with changes in GWS and pore compaction of groundwater depletion.Hence, it could be concluded that, for the sediment-located stations nearby deep groundwater depression cones, the subsidence caused by pore compaction are much larger than the elastic uplift induced by groundwater depletion.This indicates that the groundwater depletion in deep confined aquifers is the major triggering factor for inducing land subsidence in middle-east plain of NCP, which is consistent with the reported results of Guo et al. 2015 [46].Long-term uplift or decrease rates of selected CGPS stations; note that the length of arrow is only serves as an index.The background is long-term rate of groundwater level anomalies throughout NCP, which is inferred by fitting a continuous, curved surface to the in situ groundwater level observations using GMT program Surface.Long-term uplift or decrease rates of selected CGPS stations; note that the length of arrow is only serves as an index.The background is long-term rate of groundwater level anomalies throughout NCP, which is inferred by fitting a continuous, curved surface to the in situ groundwater level observations using GMT program Surface.

Characteristics of Shallow and Deep Groundwater Changes in NCP
Figure 4b shows the monthly variation of deep confined groundwater level and shallow unconfined groundwater level.The rate of shallow groundwater decreases at a rate of −0.48 m/a till the end of 2010 and then rebounds, while deep groundwater declines at a rate of −0.73 m/a up to the end of 2012, indicating that deep groundwater suffers more severe groundwater depletion than shallow groundwater in NCP, which is in agreement with reported results of Hu et al. (2016) [40].For a better investigation of the shallow and deep groundwater variability, we further compare the in situ groundwater level observations with precipitation.As shown in Figure 4c, the long-term mean precipitation in NCP is 578 mm.The shallow groundwater level decreases when the annual precipitation in NCP is less than the long-term mean (e.g., 2009-2010), whereas it increases when precipitation exceeds the long-term mean (e.g., 2011-2012).The fluctuation of shallow groundwater level consists with the variation of precipitation.For deep groundwater, the situation is quite different; the deep groundwater level continued declining from 2009 to 2012.Overall, it implies that the shallow GWS variability is largely regulated by precipitation, while the deep GWS variability in NCP is less influenced.
The corresponding response of GPS results is the difference between the times when their displacements reach peak value.In order to avoid the influence of geological condition (bedrock or sediment) of CGPS stations, here the discussion is only focused on the bedrock-located CGPS stations.As shown in Figure 5, some CGPS (e.g., BJSH) stations located in PP attain their peaks about one month earlier than the ones (e.g., HETS) located in ECP, which is most likely due to the different depletion situation and renewal characteristics over the heterogeneous aquifers in PP and ECP.For piedmont plain, which is at the feet of Taihang Mountains, the water supply is mainly from the shallow fresh groundwater stored in the upper unconfined aquifers and is largely influenced by precipitation, as discussed before.Besides, the stations in this area can receive more lateral recharge from the Taihang Mountains than the stations in east central plain due to the decrease of hydraulic gradients and transmissivity from the west to east in NCP [39].We can conclude that, compared with ECP areas, the PP areas are easier to recover from the groundwater extraction because they are more regulated by precipitation and receive more lateral recharge from Taihang Mountains.Therefore, the minimum and maximum of hydrologic mass loading in PP occur earlier than those of ECP, leading to the peak values of some CGPS stations in PP arriving earlier than the ones in ECP.

Conclusions
This study demonstrates the potential of GPS to detect vertical displacements induced by groundwater depletion over the heterogeneous aquifers in North China Plain (NCP).Comparisons between GPS-measured and GRACE-modeled seasonal vertical deformation show good correlation after correcting the bias error of GRACE data using the multiplicative factor, which confirms that the seasonal vertical displacements in NCP are mainly caused by groundwater mass loading variation, since the TWS anomalies are mainly caused by GWS variations in this area.
The seasonal vertical displacements of individual CGPS stations reflect the groundwater pumping activities for agriculture irrigation in spring, leading to less water storage in summer than in winter.The extraction of groundwater causes two processes: one is the land subsidence due to compaction of sediments in the aquifer, while the other is the uplift due to the reduction of the surface loading.In our study, the first process probably affects sediment stations rather than bedrock stations.Thus, the linear trend rate of the sediment-located CGPS stations is decreasing and mostly higher (−44~−17 mm/a) than the slight uplift or subsidence rate (−0.69~1.5 mm/a) of the other bedrock-located CGPS stations.Besides, the geological condition (bedrock or sediment) in CGPS stations, as well as their geographical location (PP or ECP) with distinct characteristics of shallow and deep groundwater renewal, differ with regard to the time when CGPS stations attain their maximum height, which has been discussed in Section 3.Both of the GRACE-derived and in situ-observed GWS anomalies show apparent variation from "obviously declining" (2009-2010) to "relatively stable" (2011-2013) in NCP, which also could be inferred from the linear trend changes of vertical displacements of three IGS CGPS stations.These results imply the slowdown of groundwater depletion rate in recent years.Overall, GPS-derived vertical displacements have the ability to detect groundwater variations in NCP.Moreover, to some extent, their result can reveal the distinguished characteristics of groundwater changes in different sub-regions.

Figure 1 .
Figure 1.Study region of NCP with two sub-regions: the piedmont plain and the east central plain.The distribution of CGPS stations (black and magenta diamonds for bedrock and sediment CGPS stations, respectively) and groundwater level monitoring well stations of different types (blue and yellow dots for shallow unconfined and deep confined well, respectively) in the NCP are shown here.1° × 1° GRACE grids are drawn with oblique lines, following the regional boundary.All the calculations (e.g., GRACE-derived TWS and GWS, GLDAS-modeled SWS, and precipitation) related to the region are based on the extents of GRACE grids.The background is the region's digital elevation model (DEM) information.

Figure 1 .
Figure 1.Study region of NCP with two sub-regions: the piedmont plain and the east central plain.The distribution of CGPS stations (black and magenta diamonds for bedrock and sediment CGPS stations, respectively) and groundwater level monitoring well stations of different types (blue and yellow dots for shallow unconfined and deep confined well, respectively) in the NCP are shown here.1 • × 1 • GRACE grids are drawn with oblique lines, following the regional boundary.All the calculations (e.g., GRACE-derived TWS and GWS, GLDAS-modeled SWS, and precipitation) related to the region are based on the extents of GRACE grids.The background is the region's digital elevation model (DEM) information.

Figure 2 .
Figure 2. Annual (left) and semi-annual (right) signals (right) of vertical displacements measured by GPS and modeled by GRACE in NCP.

Figure 2 .
Figure 2. Annual (left) and semi-annual (right) signals (right) of vertical displacements measured by GPS and modeled by GRACE in NCP.
shows the daily observations of three bedrock-located IGS CGPS stations (BJFS, BJSH, JIXN), since other CGPS stations' data are only available after July 2010.The long-term linear trends of crustal uplift or subsidence are mainly dominated by surface mass loading and tectonic processes.The analysis of the three IGS CGPS stations (BJFS, BJSH, JIXN) suggests the slight change of linear trend around 2010.From 2009 to 2010, the linear trend of vertical displacements is uplift, which is consistent with the reported results of Wang et al. (2017)[17] and might imply the unloading of hydrological mass in NCP.However, the trend changes to lightly downward or relatively stable after the end of 2010, indicating the increase of hydrological mass loading in this area.These facts reveal that the severe water condition in NCP might has been alleviated in recent years.Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 13

Figure 3 .
Figure 3. Daily observations and their linear trends (red dashed line) of the vertical displacements at three IGS CGPS stations.Figure 3. Daily observations and their linear trends (red dashed line) of the vertical displacements at three IGS CGPS stations.

Figure 3 .
Figure 3. Daily observations and their linear trends (red dashed line) of the vertical displacements at three IGS CGPS stations.Figure 3. Daily observations and their linear trends (red dashed line) of the vertical displacements at three IGS CGPS stations.

Figure 4 .
Figure 4. (a) The monthly time series of modeled SMS (magenta line) and GRACE-derived TWS (black line) and GWS (blue line) anomalies in NCP from 2009 to 2013.The GRACE-derived GWS anomalies are compared with that from the in situ measurements (red line); (b) observations of groundwater level (GWL) monthly anomalies of deep confined well (blue line) and shallow unconfined well (red line) during 2009 to 2013, and the blue and red dashed line are the linear trend of in situ shallow and deep GWL observations, respectively, with trend rates shown above.(c) Precipitation is shown as the annual total (Figure 4c above) and the monthly observations (Figure 4c below).The precipitation data are provided by the China Meteorological Administration (http://cdc.cma.gov.cn/).

Figure 4 .
Figure 4. (a) The monthly time series of modeled SMS (magenta line) and GRACE-derived TWS (black line) and GWS (blue line) anomalies in NCP from 2009 to 2013.The GRACE-derived GWS anomalies are compared with that from the in situ measurements (red line); (b) observations of groundwater level (GWL) monthly anomalies of deep confined well (blue line) and shallow unconfined well (red line) during 2009 to 2013, and the blue and red dashed line are the linear trend of in situ shallow and deep GWL observations, respectively, with trend rates shown above.(c) Precipitation is shown as the annual total (Figure 4c above) and the monthly observations (Figure 4c below).The precipitation data are provided by the China Meteorological Administration (http://cdc.cma.gov.cn/).

Figure 5 .
Figure 5. Histograms of peak uplift phase (binned by month) for selected CGPS stations in NCP, red bars for CGPS stations located on bedrock and blue bars for ones located on sediment.

Figure 5 .
Figure 5. Histograms of peak uplift phase (binned by month) for selected CGPS stations in NCP, red bars for CGPS stations located on bedrock and blue bars for ones located on sediment.

Figure 6 .
Figure 6.Long-term uplift or decrease rates of selected CGPS stations; note that the length of arrow is only serves as an index.The background is long-term rate of groundwater level anomalies throughout NCP, which is inferred by fitting a continuous, curved surface to the in situ groundwater level observations using GMT program Surface.

Figure
Figure 4b shows the monthly variation of deep confined groundwater level and shallow unconfined groundwater level.The rate of shallow groundwater decreases at a rate of −0.48 m/a till the end of 2010 and then rebounds, while deep groundwater declines at a rate of −0.73 m/a up to the end of 2012, indicating that deep groundwater suffers more severe groundwater depletion than shallow groundwater in NCP, which is in agreement with reported results of Hu et al. (2016) [40].For a better investigation of the shallow and deep groundwater variability, we further compare the

Figure 6 .
Figure 6.Long-term uplift or decrease rates of selected CGPS stations; note that the length of arrow is only serves as an index.The background is long-term rate of groundwater level anomalies throughout NCP, which is inferred by fitting a continuous, curved surface to the in situ groundwater level observations using GMT program Surface.

Table 1 .
Specific information of selected CGPS stations.

Code Topography and Geomorphology Station Geologic Characteristic (Basement Type)
Mt. is the abbreviation of Mountains.2HAHB is located on the boundary between PP and ECP.

Table 2 .
Annual amplitudes and phases, linear trend fit of vertical displacements derived by GPS and GRACE for selected CGPS stations, and their correlation coefficients of annual components.

Table 2 .
Annual amplitudes and phases, linear trend fit of vertical displacements derived by GPS and GRACE for selected CGPS stations, and their correlation coefficients of annual components.