Evaluation of Groundwater Storage Depletion Using GRACE/GRACE Follow-On Data with Land Surface Models and Its Driving Factors in Haihe River Basin, China

Groundwater storage (GWS) in the Haihe River Basin (HRB), which is one of the most densely populated and largest agricultural areas in China, is of great importance for the ecosystem environment and socio-economic development. In recent years, large-scale overexploitation of groundwater in HRB has made it one of the global hotspots of GWS depletion. In this study, monthly GWS variations in HRB from 2003 to 2020 were estimated using the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) data in combination with three land surface models (LSMs) from the Global Land Data Assimilation System (GLDAS). The results show the following: (1) HRB suffered extensive GWS depletion from 2003 to 2020, which has been aggravated since 2014, with a mean rate of 1.88 cm·yr−1, which is equivalent to a volume of 6 billion m3·yr−1. The GWS depletion is more serious in the plain zone (−2.36 cm·yr−1) than in the mountainous zone (−1.63 cm·yr−1). (2) Climate changes are excluded from the reasons for GWS depletion due to annual precipitation and evaporation being close to normal throughout the period. In addition, GWS changes show a low correlation with meteorological factors. (3) The consumption of groundwater for irrigation and land use/cover changes have been confirmed to be the dominant factors for GWS depletion in HRB. (4) The effects of inter-basin water transfer projects cannot be obviously observed using the GRACE and GRACE-FO; more inter-basin water transfers are needed for recovering the GWS in HRB. Therefore, it is imperative to control groundwater exploitation and develop a more economical agricultural irrigation structure for the sustainability of groundwater resources in HRB.


Introduction
Groundwater, the third largest water storage after the ocean and the cryosphere, is the world's largest freshwater resource for agriculture, industry, public supply, and ecosystems [1,2]. However, some regions consume more groundwater than is recharged, causing groundwater storage (GWS) depletion [3]. In the past several decades, GWS depletion has occurred in many areas throughout the world, such as in north China [4], northwest India [5], and the Central Valley in the United States [6]. Globally, the magnitude of GWS depletion may be so large that it has become a measurable contributor to rising sea levels. GWS depletion has already caused some negative geological impacts on the environment, such as groundwater funnels, land subsidence, soil salinization, seawater intrusion, groundwater quality deterioration, and, in some cases, ecosystem degradation [7][8][9][10]. However,

Study Area
Haihe River Basin (HRB) is located between 112 • -120 • E and 35 • -43 • N with an area of about 31.8 × 10 4 km 2 . HRB is a political, economic, and cultural center, adjacent to the Bohai Sea in the east, Yellow River in the south, Taihang Mountain in the west, and the Inner Mongolia Plateau in the north. According to the administrative division, HRB covers Beijing City, Tianjin City, most of Hebei Province, part of Shandong Province, Henan Province, Shanxi Province, Inner Mongolia Autonomous Region (Neimeng), and Liaoning Province (Figure 1). HRB has a typically semi-humid and semi-arid continental monsoon climate. Affected by the continental air mass of Mongolia, the temperature of spring rises quickly with high winds and strong evaporation. Influenced by the marine air mass, it is humid in summer, with high temperatures and precipitation. Autumn brings less precipitation, and winter is controlled by the Siberian continental air mass with less snow. The annual average temperature is 1.5-14°C, and the annual average relative humidity is 50-70%. The annual average precipitation is 535 mm, with over 70% of precipitation falling in the summer months of June to August, and annual average surface evaporation of 850-1300 mm [34].
There are the following three river systems in HRB: Haihe River System, Luanhe River System, and Tuhaimajiahe River System. As the biggest, the Haihe River system contains Jiyun River, Chaobai River, Beiyun River, and Yongding River in the northern part, and Daqing River, Ziya River, and Zhangwei River in the southern part. The Luanhe River system and Tuhaimajiahe River system are in the northern and southern parts of HRB, respectively. As the largest land use type in the region, agriculture accounts for over 90% of the arable lands. In addition, over 70% of the arable lands are under continuous wheat-maize rotation farming. HRB grain production accounts for 10% of China's total grain production, including 30% of wheat production and 20% of corn production [35,36].
The bedrock of HRB consists of Archean gneiss and Proterozoic carbonate rocks overlain by thick Tertiary and Quaternary deposits. In a plain area, groundwater mainly occurs in Quaternary unconsolidated strata. The Quaternary aquifers in HRB are recognized as one of the world's 37 large aquifer systems in the Worldwide Hydrogeological Mapping and Assessment Program (WHYMAP) [37]. The Quaternary aquifers in HRB consist of four major aquifer units: Holocene series (Q4), late Pleistocene series (Q3), middle Pleistocene series (Q2), and early Pleistocene series (Q1). Q4 is an unconfined aquifer, 40-170 m in depth, consisting of coarse-grained sand in the piedmont plain and fine-grained sand in the littoral plain. Q3 is a shallow confined aquifer, 120-70 m in depth, consisting of sandy gravel, and medium to fine sand. Q2 is a confined aquifer, 250-350 m in depth, consisting of sandy gravel in the piedmont plain, and medium to fine sand in the central and the littoral plains. Q1 is 400-600 m in depth and consists of cemented sandy gravel and a thin layer of weathered sand. Both aquifers Q1 and Q2 are the targets of exploitation.
Since the 1970s, rapid agricultural and industrial development has resulted in a great demand for water resources [33]. In 2000-2019, groundwater accounted for 81% of the total water supply in HRB, and 79%, 44%, and 83% in Beijing City, Tianjin City, and Hebei Province, respectively [38]. Groundwater in HRB has been in a state of unsustainable development. In the year 2000, 7.7 × 10 9 m 3 of groundwater was overexploited, of which 3.9 × 10 9 m 3 was from shallow groundwater and 3.8 × 10 9 m 3 was from deep groundwater, forming 60 × 10 3 km 2 shallow groundwater overexploitation areas and 56 × 10 3 km 2 deep groundwater overexploitation areas.

GRACE/FO Data
There are two primary GRACE/FO datasets, the spherical harmonic coefficients (SHCs) and the mass concentration solutions (mascon). In this study, two versions of GRACE/FO mascon products, CSR-M (from the Center for Space Research) and JPL-M (from the Jet Propulsion Laboratory) were used to obtain the gridded TWSA. These mascon data are now developed to the 6th generation (RL06), which was recently released to minimize the common north-south striping errors and noises. CSR-M provides monthly TWSA data with a grid of 0.25 • × 0.25 • that are available from April 2002 to June 2017 for GRACE and from June 2018 to December 2020 for GRACE-FO [39]. TWSA are relative to the 2004-to-2009 mean baseline. CSR-M could be used without any post-processing. JPL-M provides monthly TWSA data with a grid of 0.5 • × 0.5 • with the same time span as CSR-M. JPL-M needs to multiply the scale factor to recover the leakage signals. In both products, the degree 2 order 0 coefficients have been substituted with more accurate estimations from satellite laser ranging (SLR) [40,41], and the degree 1 coefficient has been added [42]. Glacial isostatic adjustment (GIA) correction has been applied based on the model from Geruo et al. [43]. The leakage of continental hydrology is corrected following the method proposed by Chen et al. [44]. The missing GRACE/FO data are filled by using the cubic-spline interpolation.

Land Surface Models
The National Aeronautics and Space Administration (NASA) and National Oceanic and Atmospheric Administration (NOAA) jointly developed the Global Land Data Assimilation System (GLDAS). The land surface models (LSMs) in the GLDAS adopt advanced data assimilation technology to integrate satellite and ground observation data and can simulate the global distributions of land surface states/fluxes in real time. In this study, the monthly soil moisture content (SM), snow water equivalent (SWE), and total canopy water storage (CWS) from the following three LSMs were used as auxiliary data to estimate the GWSA: the Catchment Model (CLSM), NOAH Model (NOAH), and Variable Infiltration Capacity Model (VIC) from the GLDAS (v2.1). The simulated SM in the CLSM, NOAH, and VIC only hold representatives for the 343 mm, 200 mm, and 190 mm soils in depth. The time span ranges from January 2000 to December 2020 with a spatial resolution of 1 • × 1 • . In addition, evapotranspiration and precipitation flux were also used to calculate the regional evapotranspiration and precipitation.

GWSA Based on GRACE/FO and LSMs Data
GWSA can be isolated from TWSA derived by the GRACE/FO and other components estimated by auxiliary data using the following formula: According to the statistical data from Haihe River Basin Water Resources Bulletins [38], the water storage in reservoirs and lakes are very small in HRB. Therefore, the surface water storage anomalies (SWSA) consist of SM anomalies (SMA), SWE anomalies (SWEA), and CWS anomalies (CWSA) in HRB. SMA, SWEA, and CWSA are obtained by subtracting the mean value of 2003-2020 from the original data. Due to the inconsistent spatial resolution between LSMs data and GRACE/FO data, LSMs data and GRACE/FO data are downscaled at a spatial resolution of 0.025 • .

Time Series Decomposition
The time series is composed of seasonal signal, trend signal, and residuals based on the non-parametric STL (Seasonal-Trend decomposition procedure based on Loess) method.
The STL method is a robust and computationally efficient approach to detect nonlinear patterns in trend estimates, based on a locally weighted regression [45]. The Mann-Kendall trend test with Sens's slope estimator is a widely used non-parametric trend test for analyzing climatic, hydrological, and environmental time series data, used to assess the statistical significance of the trend signal in GWSA at the pixel scale.

Groundwater Storage Deficit Index
The monthly groundwater storage deficit (GWSD) is used to indicate the hydrological characteristics as follows: where GWSD i,j and GWSA i,j represent the GWSD and GWSA for the jth month in the year I, and GWSA j is the long term mean of GWSA for the jth month. A negative GWSD indicates that GWS is less than the normal value. GWSDI is the normalized GWSD by the zero mean normalization method as shown in the following Equation (4): where µ and σ are the mean and standard deviation of GWSD time series, respectively.

Uncertainty Assessments
The TWSA uncertainty is based on the method proposed by Landerer and Swenson [46]. First, the linear trends and seasonal components were removed from TWSA to obtain the residuals. Then, the inter-annual signals were removed by fitting a 13-month moving average to the residuals. The root mean square (RMS) of the residuals approximated the uncertainty in TWSA. For SWSA, the uncertainty was estimated from the monthly standard deviation (STD) among the three GLDAS models. Considering the error propaga-tion law during the least squares fit and uncertainties in TWSA and SWSA, the uncertainty in the GWSA was estimated by using the following equation: Compared with the TWSA, the SWSA obtained from LSMs show a non-significant increase or decrease trend from 2003 to 2020 with smaller amplitudes (Figure 2b). In addition, the monthly SWSA from the CLSM, NOAH, and VIC show obvious differences. The CLSM-estimated SWSA is larger than that estimated from NOAH and VIC ( Figure 2b). The differences are mainly caused by the differences in the soil depths. The uncertainties of the SWSA are 3.19, 3.23, and 2.25 cm for the CLSM, NOAH, and VIC, respectively.     Figure 5 shows the spatial patterns of GWSA in HRB in the period of 2003-2020. On one hand, the groundwater depletion in HRB is serious year by year. On the other hand, the distribution of GWSA is uneven in HRB, and the southern area is more depleted than that in the northern area. There are two GWS depletion centers: one is on the border of Hebei Province, Shanxi Province, and Henan Province, including Xingtai City, Handan City, Changzhi City, and Hebi City. The other is in the east of Hebei Province, including Hengshui City, Cangzhou City, and DeZhou City.   Figure 6). In fact, the previous studies conducted in HRB also found inconformity between the GWSA based on the GRACE/FO and LSMs with the in situ groundwater level [47,48]. There are four reasons for these differences.

The Spatial Patterns of GWSA in HRB
(1) The GRACE/FO examines a large spatial region, whereas monitoring wells are for point scale [5]. (2) Groundwater monitoring wells only monitor water level changes in a certain aquifer, while the GRACE/FO monitors the comprehensive performance of all aquifers.
(3) The unevenly distributed groundwater monitoring wells do not have adequate spatial resolutions to reflect the holistic situation of Beijing City. (4) The coarse spatial resolution of the GRACE/FO is more suitable for a region with an area larger than 20 × 10 4 km 2 .

Comparisons with the Data from Water Resources Bulletins
Haihe River Basin Water Resources Bulletins give the shallow groundwater storage variations in the plain zone [38]. Hebei Water Resources Bulletins give the shallow groundwater storage variations and groundwater depth variations in the deep aquifers [49]. Figure 7a shows the trend of GWSA in HRB based on the GRACE/FO, which is similar to the trend of GWSA in the shallow aquifer from Haihe River Basin Water Resources Bulletins, with a correlation coefficient of 0.90. Figure 7b shows the trend of GWSA in Hebei Province based on the GRACE/FO, which is similar with the trend of GWSA in the shallow aquifer from Hebei Water Resources Bulletins, with a correlation coefficient of 0.85. This phenomenon indicates that the shallow aquifer is the main force of groundwater depletion. However, GWS depletion derived from the GRACE/FO has been aggravated since 2014, while this phenomenon is not shown in the shallow groundwater storage changes in the Water Resources Bulletins. This indicates that the groundwater depletion has been deteriorative in the deep aquifer since 2014. This behavior is consistent with the previous geological survey results; the shallow groundwater funnel is over 20,000 km 2 , and the deep groundwater funnel is over 70,000 km 2 in the plain area of the north. In addition, the GWSA in Hebei Province based on the GRACE/FO show a linear correlation with the deep groundwater depth changes in Cangzou City, Hengshui City, and Xingtai City (Figure 7c-e). Since the age of deep groundwater in the North China Plain is up to 25 ka B.P [50], it is very slow to recover from overexploitation.  Table 2, our results show TWS/GWS change trends that are consistent with prior studies. However, there are still some subtle differences in the change rate among different studies. These are mainly due to the differences in periods of data, kinds of data version, data combinations, post-processing methods, and so on. The superiority of our results lies in their deriving from the latest data version and multi-data combinations. On one hand, studies have proved that mascon solutions are more accurate compared with spherical harmonics solutions, although there was no post-processing [53]. In addition, this study not only uses the data from GRACE but also from GRACE-FO in HRB and covers a longer period. On the other hand, this study considers and compares the TWSA based on the latest mascon solutions (RL06) from CSR and JPL, and the SWSA based on three LSMs, and finally forms three combinations, the GRACE/FO with the CLSM, NOAH, and VIC, respectively.

The Spatial Trends of GWS Depletion in HRB
According to hydrogeological characteristics, HRB consists of the following two subzones: mountain area and plain area. The plain area consists of a piedmont pluvial plain, a central alluvial plain, and a coastal plain [50]. GWS decreases more quickly in the plain area (−2.36 cm·yr −1 ) than in the mountain area (−1.63 cm·yr −1 ) (Figure 8). The spatial patterns of GWS changes show significant groundwater depletion in the piedmont plain and central plain of HRB. The spatial patterns of GWS depletion are opposite to the precipitation recharge patterns. The maximum recharge rates are 250, 200, and 150 mm·yr −1 for the piedmont plain, the central plain, and the coastal plain, respectively [16].
The spatial pattern of the GWS trend is further estimated based on the administrative division. The northern part, including Neimeng and Liaoning Province, shows the lowest GWS depletion rate, less than 1 mm·yr −1 . Beijing City also shows a small GWS depletion rate of around 1 mm·yr −1 . The southern part, including Henan Province, Shanxi Province, and Shandong Province, shows serious GWS depletion. The mountain area in Shanxi Province is the largest coal-producing region in China. This large-scale coal mining would destroy the aquifers and result in the observed groundwater depletion [4].

Impacts of Climate Changes on GWS Depletion in HRB
Two driving factors account for the GWS depletion: climate changes (mainly the changes of temperature, precipitation, and evaporation) and human activities (water exploitation, hydraulic projects, land use/covers changes, etc.). Based on the method of Andersen et al. [62], the GWS changes (GWSCs) for a given area can be expressed as follows: We consider the impacts of climate changes on GWSC in HRB in the following four aspects: (1) Due to the weakening monsoon winds over the past 60 years, many studies have found decreasing trends of annual precipitation in HRB by hydro-meteorological stations [63][64][65]. However, the annual precipitation from LSMs indicates a non-significant increasing trend during 2000-2020 (6.44 mm·yr −1 , R 2 = 0.28). The increased annual precipitation is almost offset by evaporation, which also increases (5 mm·yr −1 , R 2 = 0.62). Net recharge is the subtraction between precipitation and evapotranspiration (P-ET). In HRB, the net recharge remained constant during 2000-2020 (1.44 mm·yr −1 , R 2 = 0.03). Therefore, the persistent downward trend of GWS in HRB year by year indicates that climate changes are not the major reasons for GWS depletion.
(2) Figure 9 shows that GWSC has a low correlation, even different phases, with precipitation, evaporation, and net recharge in HRB. Previous studies have shown that the maximum correlation coefficient between precipitation and GWS is low in HRB and has a time lag of nine months [66]. In this study, the maximum cross correlation coefficient between net recharge and GWSC is 0.16 with a lag of two months. The lag time corresponds with the seasonal pattern of GWSA in HRB. (3) As GWSC might have a strong relationship with dry/wet conditions, the groundwater storage deficit index (GWSDI) is compared with the self-calibrating Palmer drought severity index (scPDSI). The scPDSI is a comprehensive drought index derived from a variety of meteorological, agricultural, and hydrological variables to give a comprehensive picture of drought [67]. The scPDSI data are provided by the Climatic Research Unit at the University of East Anglia with a spatial resolution of 0.5 • × 0.5 • [68]. The low correlation between the scPDSI and GWSDI series also indicates that climate changes are not the major factors for GWS depletion in HRB ( Figure 9). (4) The spatial distribution of GWS changes is different from the precipitation spatial distribution [32,38,49]. Based on the above four aspects, climate changes are not the major reasons for GWS depletion in HRB.

Impacts of Human Activities on GWS Depletion in HRB
In this part, we consider the impacts of human activities on GWS depletion in HRB. The groundwater depletion is caused by the excess in groundwater exploitation compared to groundwater recharge [69]. The current global overexploitation of groundwater resources exceeds their recharge by 3.5 times [70]. Groundwater extraction increased by about 2.5 billion m 3 per year to meet the need in China [71]. According to the Haihe River Water Conservancy Commission [38], the mean water supply was 37.7 billion m 3 during 2000-2020, which is 1.3 times the mean total water resources (28.3 billion m 3 ), and 1.7 times the mean groundwater resources ( Figure 10). During 2003-2020, the accumulative water overexploitation was 198 billion m 3 . Therefore, overexploitation is one of the major reasons for groundwater depletion in HRB. Irrigation is an important factor affecting the change of terrestrial water storage [57]. Numerous studies have demonstrated that GRACE can detect GWS changes in arable areas, such as in northern India and Pakistan [3,72], the High Plains aquifer [73], the Central Valley [74], and the alluvial aquifer along the Mississippi River [27]. According to the latest version of the "Global Map of Irrigated Areas" (version 5) from FAO's Global Information System on Water and Agriculture, irrigation in HRB mainly depends on groundwater rather than surface water. Cao et al. estimated that irrigation has severely depleted groundwater resources at a rate of around 4 km 3 ·yr −1 since the 1960s in the North China Plain [75]. In addition, the annual two crops system is the most popular planting pattern in HRB, with winter wheat and summer maize. The annual average water consumption of winter wheat is 450 mm, which is two or three times the average precipitation amount during its growth period [76]. The mean irrigation was 24.6 billion m 3 , consuming 65% of the water supply, and exceeded the groundwater resources in 2000-2020 [38]. During this period, irrigation accumulatively used 45 billion m 3 of water more than the groundwater resources in HRB. The excess water is compensated by groundwater, leading to continuous groundwater depletion. Therefore, agricultural activities and the consequent irrigation are important factors for GWS depletion in HRB.
Land use/cover changes also show influences on GWS changes. Firstly, based on remote sensing images, the area of arable land and grass land decreased by 7.85% (160,354 to 147,763 km 2 ) and 4.54% (61,242 to 58,462 km 2 ) from 2000 to 2020 in HRB, while the more water-intensive forest and urban land increased by 2.75% (60,444 to 62,104 km 2 ) and 54.40% (25,120 to 38,784 km 2 ), respectively. Additionally, the growth of deep-rooted vegetation caused by the implementation of ecological restoration will consume more water and reduce the groundwater level. The switch toward more profitable crops also causes increased groundwater consumption.
Moreover, previous studies show that urban expansion is a non-ignorable factor in GWS depletion in the Jing-Jin-Ji region [20]. During the past few decades, HRB has been undergoing significant changes in urbanization due to social-economic development. By analyzing the data from the Resource and Environment Science and Data Center, we found that the urban area increased from 25,120 km 2 in 20,00 to 27,104 km 2 in 2005, 28092 km 2 in 2010, 29,376 km 2 in 2015, and 38,784 km 2 in 2020. On one hand, the obvious urban expansion in HRB leads to a decrease in infiltration with the low permeability of impervious surfaces. On the other hand, continuous urban expansion increases the water consumption, resulting in the overexploitation of groundwater and the aggravation of groundwater depletion. Figure 11 shows that there is an obviously negative relationship between the population and GDP (Gross Domestic Product) with GWS changes in Beijing City, Tianjin City, and Hebei Province. There are the following two major inter-basin water transfer projects in HRB: the Middle Route of the South to North Water Diversion Project (SNWDP) and the Yellow River transfer scheme. In total, 9.5 billion m 3 of water is delivered to the North China Plain annually, which will significantly change the water supply structure [77]. In this study, we found that only in 2019, when the watershed water regulation was more than 10 billion m 3 , was the trend of groundwater depletion alleviated. In addition, the GWS changes had no relation to the inter-basin water transfer (Figure 12), which proves from the side that the GWS depletion in HRB is serious and more inter-basin water transfers are needed. All in all, human activities influence GWS changes in the aspects of groundwater exploitation, agricultural activities and irrigation, land use/cover changes, and city expansion and inter-basin water transfer projects in HRB.

Limitations in this Study
Although this study uses the latest data and considers different data combinations, we have to admit that the GWSA estimated in this study is affected by a series of limitations. These limitations are from both data processing schemes and the propagation of model errors. For example, the uncertainties in the simulations of SM, SWE, CWS changes with different LSMs are manifested, but to what levels remains unknown. In addition, the typical thickness of the unsaturated zone is as deep as 10-15 m for the North China Plain [75], while the soil depths in LSMs range from 200 to 340 cm. Thus, SMS do not account for water storage variations in deep soil. The unmatched spatial resolution of the dataset and the processing methods for each dataset could result in substantial uncertainties. Recently, data assimilation techniques have been proposed to improve GWS by assimilating the GRACE/FO observation into LSMs. We will use these new techniques in our future studies to enhance the accuracy of the GWS estimation.

Conclusions
GRACE and GRACE Follow-On (GRACE/FO) mascon data were used to estimate groundwater storage variations after removing other components from land surface models in the HRB, China. The TWSA derived from the two GRACE/FO data products (CRS and JPL) show that there was good agreement and consistency in HRB during 2003-2020, and the mean value was used as the ensemble solution of TWSA in HRB. The monthly TWSA exhibited a decreasing trend at a rate of 19.2 mm·yr −1 from 2003 to 2020, with an uncertainty of 1.74 cm. The surface water storage based on land surface models does not contribute significantly to the decrease in the TWSA. In addition, three land surface models (CLSM, NOAH, and VIC) show different SWSA. Therefore, GWS decreases at the rates of 18.7, 18.6, and 19.2 mm·yr −1 for the combination of the GRACE/FO with CLSM, NOAH, and VIC, respectively. The GWS decreases more quickly in the plain area (−23.6 mm·yr −1 ) than in the mountain area (−16.3 mm·yr −1 ). During the period of 2003-2020, the groundwater depletion was equivalent to a net loss of 108 billion m 3 of water.
The GWS depletion in HRB has more of a human than natural cause. In the scenario of climate changes, the precipitation and evaporation in the HRB were stable without a significant increase or decrease during 2000-2020, and the changes in precipitation and evaporation almost offset each other. Compared with climate changes, human activities play a dominant role in the GWS depletion in the HRB. The continuous groundwater overexploitation, agricultural activities and irrigation, and urban expansion account for the groundwater resource shortage. Consequently, we should strengthen the management system of water resources in many aspects, and vigorously develop water-saving irrigation strategies, carry out groundwater recharge, and strictly control groundwater exploitation.
Author Contributions: Y.G., F.W., and R.J. processed the data; Y.G., F.G., B.Y., and J.B. conducted the analysis and drafted the paper. N.X. and Q.L. have worked collaboratively on structuring and reviewing the text, as well as discussing the results. All authors have read and agreed to the published version of the manuscript.