Monitoring Terrestrial Water Storage Changes with the Tongji-Grace2018 Model in the Nine Major River Basins of the Chinese Mainland

: Data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission can be used to monitor changes in terrestrial water storage (TWS). In this study, we exploit the TWS observations from a new temporal gravity ﬁeld model, Tongji-Grace2018, which was developed using an optimized short-arc approach at Tongji University. We analyzed the changes in the TWS and groundwater storage (GWS) in each of the nine major river basins of the Chinese mainland from April 2002 to August 2016, using Tongji-Grace2018, the Global Land Data Assimilation System (GLDAS) hydrological model, in situ observations, and additional auxiliary data (such as precipitation and temperature). Our results indicate that the TWS of the Songliao, Yangtze, Pearl, and Southeastern River Basins are all increasing, with the most drastic TWS growth occurring in the Southeastern River Basin. The TWS of the Yellow, Haihe, Huaihe, and Southwestern River Basins are all decreasing, with the most drastic TWS loss occurring in the Haihe River Basin. The Continental River Basin TWS has remained largely unchanged over time. With the exception of the Songliao and Pearl River Basins, the GWS results produced by the Tongji-Grace2018 model are consistent with the in situ observations of these basins. The correlation coefﬁcients for the Tongji-Grace2018 model results and the in situ observations for the Yellow, Huaihe, Yangtze, Southwestern, and Continental River Basins are higher than 0.710. Overall, the GWS results for the Songliao, Yellow, Haihe, Huaihe, Southwestern, and Continental River Basins all exhibit a downward trend, with the most severe groundwater loss occurring in the Haihe and Huaihe River Basins. However, the Yangtze and Southeastern River Basins both have upward-trending modeled and measured GWS values. This study demonstrates the effectiveness of the Tongji-Grace2018 model for the reliable estimation of TWS and GWS changes on the Chinese mainland, and may contribute to the management of available water resources.


Introduction
The GRACE satellite mission, which was jointly developed by the National Aeronautics and Space Administration (NASA) and the German Space Flight Center (GSFC), was successfully launched in March 2002 [1]. As of October of 2017, GRACE satellites have collected~15 years of observational data, which provides important data support for monitoring large-scale mass changes in various surface features. Previous studies have used GRACE observational data as constraints in model inversions of temporal gravity fields.
Research groups working at locations such as the Center for Space Research (CSR) of the University of Texas at Austin, Geo Forschungs Zentrum Potsdam (GFZ), NASA Jet Propulsion Laboratory (JPL), and Graz University of Technology (UT Graz) have all employed a dynamic approach to the gravity field inversion [2,3]. Using a modified short-arc approach, Chen et al. produced the Tongji-GRACE01, Tongji-GRACE02, and Tongji-Grace2018 models [4][5][6]. For the HUST-Grace2016 model, the group at the Huazhong University of Science and Technology relied on a modified dynamic approach [7]. These temporal gravity field models are available on the ICGEM website [8]. Additionally, researchers at CSR and JPL developed GRACE mascon solutions that are based on the mass concentration method [9,10]. Due to the availability of the GRACE satellites to monitor changes in surficial features at a spatial resolution of~300 km and at an equivalent water height (EWH) of 1 cm, analyzing data from GRACE satellites has become one of the most effective methods of monitoring macro-scale TWS changes. Some of the most common methods of analyzing water storage changes include in situ observations, infrared or microwave remote sensing satellite data, and numerical simulations [11][12][13]. However, because there are very few in situ observations, it is challenging to continuously monitor the water storage changes in a large area. Remote sensing satellite images only capture information about storage changes in surficial water bodies, and do not provide any insight into storage changes in subterranean water bodies. Furthermore, numerical simulations are limited by the accuracy of the model parameters, which can be difficult to constrain. Using GRACE satellite data to monitor TWS changes is advantageous because this method addresses some of these shortcomings and limitations.
Because GRACE has observed TWS, including groundwater, soil moisture, and snow water, it is possible to estimate groundwater variations using other climate models and/or in situ observations for soil moisture and snow water. GRACE products combined with other hydrological data have been widely used in evaluating GWS changes in many regions, e.g., northern India [14,15], California's Central Valley and Mid-Atlantic Region of the United States [16][17][18], the North China Plain (NCP) and West Liaohe River Basin of China [19][20][21], and Australia [22]. However, there are currently few studies on TWS and GWS changes in the nine major river basins of the Chinese mainland. In this study, we used Tongji-Grace2018, the GLDAS hydrological model, in situ observations, and meteorological data (such as precipitation and temperature) to comprehensively analyze the TWS and GWS changes in each of the nine major river basins of the Chinese mainland from April 2002 to August 2016.

Study Area
Due to factors such as location, water vapor source, and topographical conditions, the spatial distribution of water resources in China is extremely uneven; water is much more plentiful along the southeast coast than it is in the northwest inland region. The eastern and southern regions are in subtropical and temperate monsoon climate areas. The wet season encompasses summer and autumn, whereas the dry season consists of spring and winter. The western and northern regions of China are located in temperate continental and highland climate zones; these climate zones are characterized by relatively low annual precipitation values. According to the spatial distribution of water resources, the Resource and Environmental Science and Data Center (RESDC) at the Chinese Academy of Sciences [23] divides China into nine major river basins: the Songliao River Basin (SLRB), the Yellow River Basin (YERB), the Haihe River Basin (HARB), the Huaihe River Basin (HURB), the Yangtze River Basin (YARB), the Pearl River Basin (PERB), the Southeastern River Basin (SERB), the Southwestern River Basin (SWRB), and the Continental River Basin (CORB) (Figure 1). The acreage, vegetation coverage, and climatic conditions of each basin are shown in Table 1.

Terrestrial Water Storage Estimates Using the GRACE Model
In this study, we used the Tongji-Grace2018 model [6] from April 2002 to August 2016, and used spherical harmonic coefficients that are truncated to the degree and order of 60. Based on the work of Wahr [24], the EWH of the global surface mass changes can be calculated using the GRACE model: where is the radius of the Earth at the equator (~6378 km), is the average density of the Earth (~5517 kg/m 3 ), is the average density of water (1000 kg/m 3 ), is the Gaussian filter coefficient, and Δ Δ are the variations in the Stokes coefficients. These variations are calculated by subtracting the average value of the model between April 2002 and August 2016 from the model for each month in that time period. Further parameter descriptions can be found in Wahr [24].

Terrestrial Water Storage Estimates Using the GRACE Model
In this study, we used the Tongji-Grace2018 model [6] from April 2002 to August 2016, and used spherical harmonic coefficients that are truncated to the degree and order of 60. Based on the work of Wahr [24], the EWH of the global surface mass changes can be calculated using the GRACE model: where a is the radius of the Earth at the equator (~6378 km), ρ ave is the average density of the Earth (~5517 kg/m 3 ), ρ w is the average density of water (1000 kg/m 3 ), w n is the  [24].

Surface Water Storage Estimates with the GLDAS Model
The GLDAS was jointly developed by the Goddard Space Flight Center (GSFC) and the National Centers for Environmental Prediction (NCEP) [25]. This model uses remote sensing satellite data and surface observational data as the inputs for four land surface process hydrological models: NOAH, VIC, CLM, and MOSAIC. Based on model simulation and data assimilation algorithms, the model outputs are the global surface state variables and the flux data. In this study, we employed the soil moisture (SM) and snow water equivalent (SWE) outputs of the GLDAS-NOAH model. These outputs, which have a spatial resolution of 1 • and a temporal resolution of one month, cover the same time period as the GRACE model.

Groundwater Storage Estimates with the GRACE and GLDAS Models
According to the terrestrial water storage balance equation, the change in the groundwater storage can be calculated by subtracting the changes in the surface water storage (SWS) (from the GLDAS hydrological model) from the change in the TWS (from the GRACE model) [26,27]: where ∆GWS is the change in the groundwater storage, ∆TWS is the change in the terrestrial water storage (via the GRACE model), ∆SMS is the change in the soil moisture storage (via the GLDAS hydrological model), and ∆SWES is the change in the snow water equivalent (via the GLDAS hydrological model).

Groundwater Storage Estimates from In Situ Observations
The in situ observations were sourced from the "China Geological Environment Monitoring Groundwater Level Yearbook", which was published by the China Institute of Geological Environment Monitoring (CIGEM) [28]. We analyzed data that was collected from 349 groundwater wells between 2005 and 2016. As shown in Figure 1, the number of groundwater wells is higher in eastern China than it is in western China. When using these in situ observations to analyze the changes in the regional groundwater, it is necessary to multiply the water level change by the reference specific yield value to convert that water level change into a regional groundwater storage change [18,29]: where ∆h is the groundwater level change value observed in a given well (m), N is the number of subareas with groundwater wells in the study region, C is the area of the grid where the groundwater wells are located, and S is the reference specific yield value for a given basin. For each basin, Table 2 shows the number of groundwater wells and the reference specific yield value [19,[30][31][32][33].

Precipitation and Temperature
Precipitation and temperature are the main factors that cause temporal and spatial changes in the regional TWS. The monthly surface precipitation and temperature data used in this study were downloaded from the China Meteorological Data Network [34]. These data are based on the precipitation and temperature values recorded at more than 2400 meteorological stations. To translate these individual data points into monthly grid cell surface temperatures or precipitation values (resolution of 0.5 • × 0.5 • ), spatial interpolation was performed using the Thin Plate Spline (TPS) algorithm in the ANUSPLIN software [30,35,36]. Figure 2 shows the spatial variation in the monthly average precipitation values and the monthly average temperature values from April 2002 to August 2016 on the Chinese mainland. cell surface temperatures or precipitation values (resolution of 0.5° × 0.5°), spatial interpolation was performed using the Thin Plate Spline (TPS) algorithm in the ANUSPLIN software [30,35,36]. Figure 2 shows the spatial variation in the monthly average precipitation values and the monthly average temperature values from April 2002 to August 2016 on the Chinese mainland.

The Scaling Factor
By processing the GRACE model via Gaussian filtering [24,37] and a de-striping error algorithm [38][39][40], we effectively suppressed most of the signal noise. However, this also reduced the overall strength of the signal. A scale factor can be applied to correct this leakage error in the model. We calculated a scaling factor k that minimizes the sum of squared residuals of and [41]: where is the time series of regional SWS changes calculated by the GLDAS hydrological model before filtering, and is the time series of regional SWS changes calculated by the GLDAS model after processing with the same filtering method as the GRACE model. By multiplying the time series of regional terrestrial water changes calculated by the GRACE model by the scaling factor k, the signal that has been weakened by the filtering process can be partially restored.

Workflow
The data processing workflow for this study is shown in Figure 3.
(1) Due to issues such as poor sensor performance and insufficient energy supply [42], there are some data gaps in the GRACE satellite data; however, using data from adjacent months, we used interpolation to retrieve these missing data points [43][44][45]. Because GRACE data cannot obtain the degree-1 coefficients, the degree-1 coefficients calculated by Swenson et al. can be used [46]. Due to their poor quality, the C20 coefficients obtained by the GRACE satellite can be replaced by the C20 coefficients determined by satellite laser ranging (SLR) observations [47,48]. To further process the GRACE model, we made

The Scaling Factor
By processing the GRACE model via Gaussian filtering [24,37] and a de-striping error algorithm [38][39][40], we effectively suppressed most of the signal noise. However, this also reduced the overall strength of the signal. A scale factor can be applied to correct this leakage error in the model. We calculated a scaling factor k that minimizes the sum of squared residuals of σ original and kσ f iltered [41]: where σ original is the time series of regional SWS changes calculated by the GLDAS hydrological model before filtering, and σ f iltered is the time series of regional SWS changes calculated by the GLDAS model after processing with the same filtering method as the GRACE model. By multiplying the time series of regional terrestrial water changes calculated by the GRACE model by the scaling factor k, the signal that has been weakened by the filtering process can be partially restored.

Workflow
The data processing workflow for this study is shown in Figure 3. lated the change in the regional groundwater level, and multiplied that value by the reference specific yield of that region to obtain the change in the local groundwater storage (△GWS). Ultimately, this value was compared to the regional GWS changes estimated with the GRACE and GLDAS models.
(4) Finally, we analyzed the precipitation, temperature, and water resources bulletin data to try to determine the factors driving the regional TWS and GWS variations.

Data Processing
Without applying any filters, the Tongji-Grace2018 model calculated the trend of the spatial changes in China's regional TWS from April 2002 to August 2016 ( Figure 4a). We observed distinct hydrological signals in the Tianshan Mountains, the North China Plain, and the southwestern part of China. However, there are still some stripe signal errors in the north-south direction, indicating that noise reduction algorithms are necessary, even in the model with low-noise. We then processed the Tongji-Grace2018 model by combining the Duan de-striping method with a Gaussian filter with radius of 150, 200, 250, and 300 km. After applying the 150 km Gaussian filter, the high-frequency noise in the model was noticeably diminished (Figure 4b). After using the Duan algorithm to remove the model stripe signal errors, our study area showed a relatively smooth hydrological signal ( Figure 4c). As shown in Figure 4c-f, the strength of the hydrological signal in the Tongji-Grace2018 model was inversely proportional to the Gaussian filter radius. Furthermore, when the Gaussian filter radius was greater than or equal to 200 km, there was no appreciable difference in the spatial variations in the hydrological signal. (1) Due to issues such as poor sensor performance and insufficient energy supply [42], there are some data gaps in the GRACE satellite data; however, using data from adjacent months, we used interpolation to retrieve these missing data points [43][44][45]. Because GRACE data cannot obtain the degree-1 coefficients, the degree-1 coefficients calculated by Swenson et al. can be used [46]. Due to their poor quality, the C20 coefficients obtained by the GRACE satellite can be replaced by the C20 coefficients determined by satellite laser ranging (SLR) observations [47,48]. To further process the GRACE model, we made glacial isostatic adjustment (GIA) corrections [49]. Due to the satellite orbit and the data post-processing methods, the temporal gravity field model includes high-frequency noise and "stripe" errors in the north-south direction. To reduce this noise, we applied Gaussian filtering [24] and Duan de-striping error algorithms [38] to the data. Additionally, we used the scale factor to correct the data leakage error. Because the aforementioned data processing removes the tidal and non-tidal atmospheric and oceanic influences from the gravity field model, it is possible to determine the EWH values for the global TWS changes. Finally, we used the latitude cosine weighted average method to calculate the regional changes in the terrestrial water storage ( TWS) [50].
(2) After expanding the GLDAS gridded data so that its degree and order match those of the GRACE model, we calculated the regional SWS changes ( SMS + SWES).
(3) For the groundwater well observational data, we pre-processed the data, calculated the change in the regional groundwater level, and multiplied that value by the reference specific yield of that region to obtain the change in the local groundwater storage ( GWS). Ultimately, this value was compared to the regional GWS changes estimated with the GRACE and GLDAS models.
(4) Finally, we analyzed the precipitation, temperature, and water resources bulletin data to try to determine the factors driving the regional TWS and GWS variations.

Data Processing
Without applying any filters, the Tongji-Grace2018 model calculated the trend of the spatial changes in China's regional TWS from April 2002 to August 2016 (Figure 4a). We observed distinct hydrological signals in the Tianshan Mountains, the North China Plain, and the southwestern part of China. However, there are still some stripe signal errors in the north-south direction, indicating that noise reduction algorithms are necessary, even in the model with low-noise. We then processed the Tongji-Grace2018 model by combining the Duan de-striping method with a Gaussian filter with radius of 150, 200, 250, and 300 km. After applying the 150 km Gaussian filter, the high-frequency noise in the model was noticeably diminished (Figure 4b). After using the Duan algorithm to remove the model stripe signal errors, our study area showed a relatively smooth hydrological signal ( Figure 4c). As shown in Figure 4c-f, the strength of the hydrological signal in the Tongji-Grace2018 model was inversely proportional to the Gaussian filter radius. Furthermore, when the Gaussian filter radius was greater than or equal to 200 km, there was no appreciable difference in the spatial variations in the hydrological signal.
After applying the Gaussian filters to the model, we determined the scale factors required to recover the leaked signals in each of the nine major river basins (Table 3). Multiplying the TWS time series from the GRACE model by the corresponding scale factor partially restored the strength of the signals that were weakened by post-processing. The magnitude of the scale factor for each basin increased with the Gaussian filtering radius ( Table 3); this observation corresponds with our earlier assertion that a larger filter radius corresponded to more signal leakage in the model. The optimal combination of filtering algorithm minimizes the signal leakage of the model. Ultimately, we selected the processing combination that includes the 150 km Gaussian filter and the Duan de-striping algorithm, applied these algorithms to the Tongji-Grace2018 model, and calculated the TWS changes in the nine major river basins of the Chinese mainland.

The Relationship between Precipitation and the GRACE and GLDAS Models
It is important to understand how the precipitation and the GLDAS model in each basin affect the calculated TWS changes in the Tongji-Grace2018 model. First, we calculated the annual phase and the annual amplitude of the TWS time series produced by the Tongji-Grace2018 model, the SWS time series produced by the GLDAS model, and the precipitation time series. Additionally, we calculated the correlation coefficient between the TWS time series and the SWS time series (i.e., the Tongji-Grace2018 and GLDAS models), the correlation coefficient between the TWS time series and the precipitation time series (i.e., the Tongji-Grace2018 model and the precipitation data), and the annual phase difference between the TWS time series and the precipitation time series (Table 4). With the exception of the Haihe and Continental River Basins, the correlation between the Tongji-Grace2018 model and the GLDAS model was greater than 0.55. In general, the correlation coefficient between the Tongji-Grace2018 model and the GLDAS model increased with the precipitation in each basin.   With respect to the precipitation changes, we found that the TWS changes in the Yellow, Haihe, and Huaihe River Basins have longer lag times of 2.6, 3.7, and 3.2 months, respectively; the corresponding correlation coefficients between the TWS time series and the precipitation time series for these basins are 0.07, −0.08, and 0.17. These low correlation coefficients could be caused by the high agricultural demand for water in the Huang-Huai-Hai area in the spring and summer, where the precipitation cannot offset the water loss due to crop irrigation [20,51]. In the Yangtze, Pearl, and Southwestern River Basins, the time lag between the precipitation and the TWS change is 1.5 months. The correlation coefficients between the precipitation time series and the TWS time series are 0.62, 0.49, and 0.54, respectively, indicating that rainfall is an important factor affecting the TWS changes in these three basins. In the Southeastern and Continental River Basins, changes in the TWS are consistent with the observed changes in the precipitation. The TWS changes in the Songliao River Basin actually precede the precipitation by~1.5 months. In this case, the cooler conditions found at higher latitudes may cause the water resources to be stored in the form of ice and snow during the winter [52].

Terrestrial Water Storage in the Nine Major River Basins
Distinct variations were observed in the TWS changes in the nine major river basins of China (Figure 4c). The TWS changes observed in the northern Songliao River Basin, the middle Yangtze River Basin, the Pearl River Basin, the Southeastern River Basin, the southeastern Southwestern River Basin, and the southeastern Continental River Basin exhibited an increasing trend, whereas the TWS changes in the southern Songliao Basin, the Haihe River Basin, the eastern Yellow River Basin, the northwestern Southwestern River Basin, and the northwestern Continental River Basin exhibited a long-term decreasing trend (Figure 4c). Figure 5 shows the changes in the TWS (Tongji-Grace2018 model), the SWS (GLDAS model), the precipitation, and the temperature for each basin. During the study period, the TWS of the Songliao, Yellow River, Haihe, Huaihe, and the Southwestern River Basins declined, whereas the TWS of the Yangtze, Pearl, and Southeastern River Basins increased. The TWS of the Continental River Basin was largely unchanged.    The TWS time series for each river basin. The red curve represents the regional TWS changes (Tongji-Grace2018 model), the purple curve represents the 13 month moving average of the TWS changes, the black straight line is the TWS trend line, the green curve is the SWS (SMS + SWES) (GLDAS hydrological model), and the blue column is the monthly rainfall. The grey column is the annual rainfall anomaly, which is obtained by subtracting the average yearly rainfall value from the annual rainfall. The yellow column is the annual temperature anomaly, which is obtained by subtracting the average yearly temperature from the annual temperature.
Although the overall TWS trend in the Songliao River Basin was consistent with the TWS increasing slowly at a rate of 1.23 ± 0.93 mm/year, the TWS in this area fluctuated wildly during certain time periods. The TWS was very low at the end of 2007 and 2012, and very high in 2013. The monthly precipitation in the Songliao River Basin in 2007 did not exceed 95 mm; the total annual precipitation that year was significantly lower than the average yearly precipitation. Additionally, the annual temperatures in 2007 and 2008 were abnormally high, resulting in increased ground evapotranspiration. Most parts of the Songliao River Basin received very little rainfall from October 2011 to February 2012. The "China Flood and Drought Disaster Bulletin" [53] indicates that the Songliao River Basin was plagued by severe droughts during these periods. After June 2012, the rainfall gradually increased, which alleviated the earlier drought conditions, caused the annual temperature to decline precipitously, and reduced the amount of ground evapotranspiration. By 2013, the yearly rainfall in the Songliao River Basin had increased significantly. In August 2013, the lower reaches of the Heilongjiang River, the upper reaches of the Nenjiang River, and the upper reaches of the Hunhe River all experienced severe flooding. During this period, the TWS of the Songliao River Basin reached its highest value since 2002. According to the data in the "Songliao River Basin Water Resources Bulletin", the total water supply of the Songliao River Basin declined slightly from 2000 to 2003, increased from 2003 to 2010, and has remained stable since 2010. The surface water supply was the single largest contributor to the total water supply. Because both the agricultural irrigation water consumption and the total water supply increased at similar growth rates, we infer that the total water supply was heavily influenced by the agricultural irrigation demand for water. From 2003 to 2012, the proportion of the total water consumption in Heilongjiang Province in the Songliao River Basin that was attributed to farm water consumption increased from 65.5% to 78.6% [52].
Due to the similar climate and environmental conditions, the TWS changes observed in the Yellow, Haihe, and Huaihe River Basins all followed the same trend of initially increasing and then decreasing. As shown in Figure 5 In July 2016, the average rainfall in the Haihe River Basin reached 140 mm. Due to the heavy rains in this time period, the Haihe River Basin quickly transitioned from a dry region to a region experiencing significant flooding. The terrestrial water supply was quickly replenished; the heavy rainfall allowed the TWS to rebound by mid-2016 ( Figure 5). Between 2004 and 2009, the rainfall surplus and deficit years in the Huaihe River Basin alternated, resulting in a relatively static terrestrial water storage trend. From 2010 to 2015, the rainfall anomalies in the Huaihe River Basin were all negative. In this area, the rainfall anomaly was not positive until 2016. During this period, the TWS of the Huaihe River Basin exhibited a downward trend. From 2004 to 2016, the TWS of the Yellow, Haihe, and Huaihe River Basins decreased at rates of −5.60 ± 0.68 mm/year, −11.61 ± 1.11 mm/year, and −12.10 ± 0.82 mm/year, respectively. Throughout the entire study period, the TWS of the Yellow River, Haihe, and Huaihe River Basins decreased at rates of −4.15 ± 0.74 mm/year, −10.96 ± 1.00 mm/year, and −9.40 ± 1.10 mm/year, respectively.
Because the Yangtze, Pearl, Southeastern, and Southwestern River Basins are all located in areas with subtropical monsoon climate conditions, the precipitation is highly seasonal; the precipitation is high in the summer and autumn, and low in the winter and spring. Unsurprisingly, the TWS of these basins also exhibit significant seasonality. In the seasons with abundant rainfall, the TWS increases; in the seasons with less rainfall, the TWS decreases. The correlation coefficients for the TWS and precipitation time series in the Yangtze, Pearl, Southeastern, and Southwestern River Basins are 0.62, 0.49, 0.62, and 0.54, respectively, indicating that rainfall is an important factor that heavily influences the TWS changes in these river basins.
According to the "China Flood and Drought Disaster Bulletin" [36], severe floods occurred in 2010, 2012, and 2016 in and around the Yangtze River Basin; the largest flood since 1998 occurred in 2016. The water level in both Dongting Lake and Poyang Lake exceeded the flood warning line for up to 29 days. Although the TWS of the Yangtze River Basin had been at a relatively high level for the previous three years (Figure 5), the low annual rainfall and unusually high temperatures resulted in very low TWS values in Yangtze River Basin in 2006, 2011, and 2013. In May 2011, the water level of Poyang Lake, Dongting Lake, and Honghu Lake were 85%, 24%, and 31% lower than the annual average water level, respectively. From 2003 to 2004, the rainfall in the Pearl and the Southeastern River Basins was relatively low, which resulted in lower TWS values in these two river basins. After 2004, the increasing precipitation allowed the TWS to gradually rebound. The Pearl River Basin suffered continuous heavy rainfall in June 2008, causing severe floods in the cities of Jiangxi, Hunan, Guangdong, Guangxi, and Guizhou. In 2009 and 2011, the lack of rainfall resulted in severe drought conditions in the Pearl River Basin. Because the drought conditions persisted for multiple years, the water storage of the main reservoirs in this area in 2011 was dangerously low; these low water storage values culminated in the saltiest tide in the Pearl River Estuary since 2005 [55]. A rare winter flood occurred in south China in November of 2015, with the Pearl River Basin showing extremely high levels of rainfall; more than 30 rivers (e.g., the Xijiang, Luoqing, Guijiang, and Xiangjiang Rivers) had water levels that surpassed their warning water level. This period of flooding corresponds to a jump in the TWS values from the Tongji-Grace2018model in this region. Terrestrial water storage in the Pearl River Basin decreased at a rate of −44.31 ± 9.68 mm/year between 2002 and 2004, and increased at a rate of 7.58 ± 1.22 mm/year between 2005 and 2016. For the entire study period, the TWS of the Yangtze, Pearl, and the Southeastern River Basins increased at rates of 4.50 ± 0.46 mm/year, 5.00 ± 1.11 mm/year, and 7.34 ± 1.02 mm/year, respectively.
Conversely, during that same time period, the TWS in the Southwestern River Basin exhibited a distinct downward trend. There has been a TWS deficit in this area since 2009. From 2002 to 2004, both the annual rainfall and the TWS of the Southwestern River Basin were relatively stable. In 2005 and 2006, both the precipitation and the TWS in the Southwestern River Basin were relatively low. Heavy rainfall in 2008 slowed the decline of the terrestrial water storage in the Southwestern River Basin. However, with very little precipitation, high annual temperatures, and intense ground evapotranspiration, droughts plagued this area from autumn of 2009 to spring of 2010 [56]. In March of 2010, the TWS of the Southwestern River Basin dropped to its lowest level since 2002.
Because the Continental River Basin is located in an area with arid and semi-arid climate conditions, this river basin is usually very dry. As shown in Figure 4, the TWS of the southern and northwestern Continental River Basin (i.e., the Tianshan region) has experienced a severe TWS loss. In contrast, the TWS of the southeastern Qinghai Province has slightly increased. The TWS variation in the Continental River Basin is within ±3 cm and is characterized by a certain interannual periodicity ( Figure 5). From 2002 to 2016, the monthly rainfall in the Continental River Basin did not exceed 60 mm, indicating that rainfall had a relatively small impact on the observed TWS changes in the Continental River Basin. In 2003, 2005, and from 2011 to 2012, the annual temperature anomaly in the Continental River Basin was negative, causing a corresponding increase in the TWS. From 2006 to 2007 and from 2013 to 2016, the TWS decreased when the yearly temperature anomaly was positive. As such, we conclude that the TWS of the Continental River Basin is greatly affected by temperature changes. Higher temperatures result in more intense ground evapotranspiration, which reduces the terrestrial water storage. From 2002 to 2016, the TWS of the Continental River Basin slowly decreased at a rate of −0.15 ± 0.33 mm/year.

GWS Results from the GRACE Models and In-Situ Well Observations
We quantified the modeled regional GWS changes by subtracting the SWS time series (GLDAS model) from the TWS time series (Tongji-Grace2018 model). We compared these GWS results to the in situ GWS observational data collected from groundwater wells in terms of RMSE, correlation coefficient, and long-term trends ( Figure 6). It can be seen from Figure 6 that the RMSE between the modeled and measured groundwater storages changes in the Haihe, Pearl, and Southeastern River Basins are relatively large, and are 6.045, 5.089, and 6.621 cm, respectively. The RMSE between the modeled and measured groundwater storages changes in other river basins are all less than 3.4 cm.
The drought conditions in 2012 caused the GWS of the Songliao River Basin to fall rapidly ( Figure 6). However, plentiful precipitation and low temperatures allowed the groundwater of the Songliao River Basin to be mostly replenished at the beginning of 2014. The modeled and measured GWS shrinkage rates in the Songliao River Basin from 2005 to 2016 are −6.52 ± 1.15 mm/year and −4.25 ± 0.76 mm/year, respectively. According to the "China Water Resources Bulletin", the ongoing demand for agricultural irrigation caused the groundwater in the basin to decrease gradually over time.
Similar to the TWS trends in these regions, the modeled GWS values for the Yellow, Haihe, and Huaihe River Basins all initially increased and then decreased. In 2003, a positive precipitation anomaly and a negative temperature anomaly led to increased surface water infiltration and groundwater replenishment in the Huang-Huai-Hai region. (GLDAS model) from the TWS time series (Tongji-Grace2018 model). We compared these GWS results to the in situ GWS observational data collected from groundwater wells in terms of RMSE, correlation coefficient, and long-term trends ( Figure 6). It can be seen from Figure 6 that the RMSE between the modeled and measured groundwater storages changes in the Haihe, Pearl, and Southeastern River Basins are relatively large, and are 6.045, 5.089, and 6.621 cm, respectively. The RMSE between the modeled and measured groundwater storages changes in other river basins are all less than 3.4cm.  The drought conditions in 2012 caused the GWS of the Songliao River Basin to fall rapidly ( Figure 6). However, plentiful precipitation and low temperatures allowed the groundwater of the Songliao River Basin to be mostly replenished at the beginning of 2014. The modeled and measured GWS shrinkage rates in the Songliao River Basin from 2005 to 2016 are −6.52 ± 1.15 mm/year and −4.25 ± 0.76 mm/year, respectively. According to the "China Water Resources Bulletin", the ongoing demand for agricultural irrigation caused the groundwater in the basin to decrease gradually over time.
Similar to the TWS trends in these regions, the modeled GWS values for the Yellow, Haihe, and Huaihe River Basins all initially increased and then decreased. In 2003, a positive precipitation anomaly and a negative temperature anomaly led to increased surface water infiltration and groundwater replenishment in the Huang-  Since 2005, the GWS of the Yellow River and Haihe River Basins has continued to decrease. A GWS deficit in the Yellow River Basin has existed since 2011. From 2013 to 2017, higher temperatures in the Yellow River Basin led to increased evapotranspiration, lower surface water levels, and a lack of groundwater recharge. From 2002 to 2014, the TWS and SWS time series of the Haihe River Basin were largely consistent with one another ( Figure 5). However, since 2014, the difference between the SWS and the TWS time series has gradually increased, indicating that the consumption of groundwater is increasing. Overall, the correlation coefficient between the TWS and the SWS time series in the Haihe River Basin is 0.350 (Table 4). It is likely that this low correlation coefficient is related to the severe groundwater deficit.
In 2014 and 2015, the Haihe River Basin experienced relatively high temperatures and very little rainfall. Specifically, from July to August 2015, most of northern China was plagued by drought. The evapotranspiration process consumes a large amount of water stored in the soil and causes the demand for agricultural irrigation and groundwater pumping to intensify. In the following three years, the GWS of the Haihe River Basin declined rapidly. According to the "China Water Resources Bulletin", from 2004 to 2016, the proportion of the water supply in the Haihe River Basin that came from groundwater sources increased from 53% to 67%. In the Haihe River Basin, the modeled and measured GWS values are largely consistent with one another. From 2005 to 2016, the correlation coefficients for the modeled and measured GWS in the Yellow River and Haihe River Basins are 0.796 and 0.608, respectively. The modeled and measured GWS shrinkage rates for the Yellow River Basin are −8.76 ± 0.72 mm/year and −5.88 ± 0.29 mm/year, respectively, and those of the Haihe River Basin are −17.65 ± 1.69 mm/year and −8.35 ± 1.76 mm/year, respectively.
Both the Huaihe and the Haihe River Basins are located in important grain-producing regions. As such, both of these regions require considerable water resources. According to the "China Water Resources Bulletin", the annual groundwater supply of the Huaihe River Basin during our study period was about 17 × 10 9 m 3 . However, because this groundwater supply only accounts for 20-30% of the total annual water supply, we conclude that the Huaihe River Basin is not heavily dependent on groundwater sources. Furthermore, the rainfall in the Huaihe River Basin is twice that of the Haihe River Basin; as such, the latter location must rely more heavily on groundwater sources for replenishment. Between 2005 and 2012, the GWS of the Huaihe River Basin decreased slowly; the modeled and measured groundwater shrinkage rates are −2.26 ± 2.59 mm/year and −1.89 ± 0.72 mm/year, respectively. With little rainfall and high temperatures, both the groundwater storage and the groundwater recharge in the Huaihe River Basin decreased precipitously; the modeled and measured GWS shrinkage rates for this period are −15.67 ± 7.12 mm/year and −5.52 ± 1.49 mm/year, respectively. From 2005 to 2016, the correlation coefficient for the modeled and measured GWS time series in the Huaihe River Basin is 0.710.
The GWS of the Yangtze, Pearl, Southeastern, and Southwestern River Basins all exhibit obvious seasonality, with lower GWS values in winter and spring, and higher GWS values in summer and autumn. From 2002 to 2004, the modeled GWS growth rates in the Yangtze and the Southwestern River Basins are 20.94 ± 7.23 mm/year and 13.71 ± 9.86 mm/year, respectively. In 2004 and 2005, the very low annual rainfall values in the Pearl and Southeastern River Basins caused the GWS to shrink at rates of −5.39 ± 11.71 mm/year and −16.11 ± 21.78 mm/year, respectively. From 2005 to 2016, due to changes in the TWS, the GWS in the Yangtze, Pearl, and Southeastern River Basins increased and the groundwater storage in the Southwestern River Basin decreased; the rates of change for the modeled GWS time series in the Yangtze, Pearl, Southeastern, and Southwestern River Basins are 4.67 ± 0.98 mm/year, 6.43 ± 1.40 mm/year, 11.35 ± 1.89 mm/year, and −9.04 ± 1.07 mm/year; the rates of change for the measured GWS time series in the Yangtze, Pearl, Southeastern, and Southwestern River Basins are 3.19 ± 0.47 mm/year, −1.37 ± 0.79 mm/year, 6.65 ± 2.37 mm/year, and −2.61 ± 0.42 mm/year. The correlation coefficients for the modeled and measured GWS time series in the Yangtze, Pearl, Southeastern, and Southwestern River Basins are 0.731, 0.147, 0.307, and 0.732, respectively. The low correlation coefficients of the Pearl and Southeastern River Basins are attributed to the lack of groundwater well stations and the poor data quality in these areas.
Our GWS results for the Continental River Basin indicate that the GWS increased at a rate of 9.98 ± 3.71 mm/year from 2002 to 2004. However, the implementation of ecological restoration projects in this region has increased the local water consumption in recent years [57]. Furthermore, the groundwater demand from both plants and humans has increased. According to statistical data from the "China Water Resources Bulletin", the groundwater supply of the Continental River Basin has been steadily increasing on a yearly basis, resulting in a decrease in groundwater in this basin. From 2005 to August 2016, the modeled and measured GWS values changed at rates of −5.31 ± 0.54 mm/year and −4.42 ± 0.32 mm/year; the correlation coefficient between these two GWS data sets is 0.741.

Conclusions
GRACE gravity satellites provide valuable data sets that can be used to monitor changes in terrestrial water resources. In this study, we utilized the Tongji-Grace2018 model, the GLDAS model, in situ well observations, and precipitation and temperature time series to comprehensively analyze the changes in the terrestrial water and groundwater storage of the nine major river basins on the Chinese mainland from April 2002 to August 2016.
The seasonal precipitation and temperature variations are the major drivers of TWS changes in the Yangtze, Pearl, Southeastern, and Southwestern River Basins. In these basins, the correlation coefficients for the TWS and precipitation time series are larger than 0.49. The TWS changes in the Songliao, Yellow, Haihe, and Huaihe River Basins are more moderate, with little to no evidence of seasonality. Other than the presence of certain interannual variations, the TWS in the Continental River Basin has been relatively stable. Overall, the TWS values of the Songliao, Yangtze, Pearl, and Southeastern River Basins are increasing, the TWS values of the Yellow, Haihe, Huaihe, and Southwestern River Basins are decreasing, and the TWS value of the Continental River Basin is balanced.
The modeled and measured groundwater storage time series of the Yangtze, Southeastern, and Southwestern River Basins all exhibit obvious seasonal changes. Except for the Pearl River Basin, the modeled and measured GWS time series are consistent with one another. From January 2005 to August 2016, the correlation between the modeled and measured groundwater storage changes in the Songliao, Pearl, and Southeastern River Basins was low, whereas the correlation coefficients for the other basins were all larger than 0.60. Due to factors such as precipitation, temperature, and human activity, the groundwater storage values of the Songliao, Yellow, Haihe, Huaihe, Southwestern, and Continental River Basins have all declined, whereas the groundwater storage values of the Yangtze and Southeastern River Basins have increased over time. Due to the large number of groundwater wells located in the Yellow, Huaihe, Yangtze, and Continental River Basins, the modeled and measured groundwater storage changes are in agreement with one another.
Through the above research results, we obtained the changes in TWS and GWS in the nine major river basins of the Chinese mainland. By comparing the groundwater changes calculated by the model with the in situ well observations, we demonstrated the effectiveness of using the Tongji-Grace2018 model combined with the hydrological model to estimate regional GWS. However, due to the limitation of the spatial resolution of the GRACE satellite gravity field model, we are currently analyzing the TWS and GWS on the Chinese mainland from a macro perspective. In future studies, it is hoped that techniques such as downscaling and data assimilation can be used to improve the resolution of the GRACE model, thereby improving the accuracy of estimating changes in terrestrial water and groundwater storage.
Funding: This research was funded by the National Natural Science Foundation of China (41674006).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.