Estimation of Terrestrial Water Storage Changes at Small Basin Scales Based on Multi-Source Data

: Terrestrial water storage changes (TWSCs) retrieved from the Gravity Recovery and Climate Experiment (GRACE) satellite mission have been extensively evaluated in previous studies over large basin scales. However, monitoring the TWSC at small basin scales is still poorly understood. This study presented a new method for calculating TWSCs at the small basin scales based on the water balance equation, using hydrometeorological and multi-source data. First, the basin was divided into several sub-basins through the slope runoff simulation algorithm. Secondly, we simulated the evapotranspiration (ET) and outbound runoff of each sub-basin using the PML_V2 and SWAT. Lastly, through the water balance equation, the TWSC of each sub-basin was obtained. Based on the estimated results, we analyzed the temporal and spatial variations in precipitation, ET, outbound runoff, and TWSC in the Ganjiang River Basin (GRB) from 2002 to 2018. The results showed that by comparing with GRACE products, in situ groundwater levels data, and soil moisture storage, the TWSC calculated by this study is in good agreement with these three data. During the study period, the spatial and temporal variations in precipitation and runoff in the GRB were similar, with a minimum in 2011 and maximum in 2016. The annual ET changed gently, while the TWSC ﬂuctuated greatly. The ﬁndings of this study could provide some new information for improving the estimate of the TWSC at small basin scales.


Introduction
Terrestrial water storage (TWS) is the sum of all forms of water storage over land surfaces and plays a vital role in the global and regional hydrological cycle [1][2][3]. It reflects all types of water stored on continents, including surface water, soil water, and groundwater [4]. TWS change (TWSC) has significant impacts on terrestrial ecosystems, human beings, and even the sea level [5]. Monitoring the spatial and temporal variabilities in TWSC is of great importance for effective water resources management and sustainable utilization [6].
The launch of the Gravity Recovery and Climate Experiment (GRACE) satellites has provided a new method for terrestrial hydrology research, which can efficiently improve the monitoring result of the changes in the water cycle at a large scale [7]. Numerous studies have shown that GRACE satellites play a key role in the large-scale monitoring of TWSC or TWS anomalies (TWSA) [8][9][10]. However, the applications of GRACE satellites are limited by the low spatial resolution that cannot be directly used in small-scale research [11][12][13]. To use GRACE data at small basins, some studies have derived TWS at small basin scales by downscaling GRACE data [14][15][16]. There are two typical types of methods to downscale GRACE data: statistical downscaling and dynamic downscaling. Scholars have used climate area and is one of the typical rainstorm areas in China. As the joint effects of the complex terrain and the East Asian monsoon climate, flood disasters frequently occur in the GRB [27]. The mean annual air temperature is around 18°C. It receives an annual precipitation between 1200 and 1700 mm, and the precipitation distribution is uneven throughout the season. More than 60% of the total precipitation is concentrated from April to September [28].

Data Description 2.2.1. Hydrometeorological Data
The meteorological and hydrological data used in this study mainly include the China meteorological forcing dataset (CMFD), CMA Land Data Assimilation System Version 2.0 (CLDAS-V2.0), and in situ data.
The CMFD is a near-surface meteorological and environmental reanalysis dataset provided by the Institute of Tibetan Plateau Research at the Chinese Academy of Science [29]. The dataset is based on the Princeton reanalysis data, Global Energy and Water Exchanges-Surface Radiation Budget (GEWEX-SRB), the goal of the Global Land Data Assimilation System (GLDAS), and Tropical Rainfall Measuring Mission (TRMM) data, combined with meteorological data from the China Meteorological Administration (CMA). Meteorological station observation data from 740 stations in the CMA were used to correct systematic departures in the background data. The spatial resolution of CMFD is 0.1 • × 0.1 • with temporal coverage from 1979 to 2018. The CMFD dataset has been widely used in the hydrometeorological forecasts and hydrological research [30][31][32].
CLDAS-V2.0 has a spatial resolution of 0.0625 • × 0.0625 • and a temporal resolution of 1 h, covering the whole East Asia region (0 • -65 • N, 60 • -160 • E). The CLDAS dataset includes a variety of land-surface forcing data, including precipitation, temperature, shortwave radiation, specific humidity, wind speed, surface pressure, and soil status variables. In addition, the CLDAS-V2.0 dataset was developed by CMA. The soil moisture storage (SMS) starts in April 2017, and the SMS data at 2 m were used as the validation data.
In situ groundwater level data of Ganzhou Hotel station (25.

GRACE Data
This study used the GRACE Mascon product from the Center for Space Research (CSR, the University of Texas at Austin) for validation of the mean TWSC estimate in the GRB. CSR Mascon (CSRM) is available at a spatial grid of 0.25 • × 0.25 • and monthly temporal resolution (www2.csr.utexas.edu/grace/ accessed on 1 July 2021). More details about this product can be found in Save et al. [33]. This product has been widely used for TWSC estimation at large scales [34][35][36].

Digital Elevation Model (DEM) Data
The global 1 arc-second SRTM V3.0 dataset (SRTMGL1) is a joint product of the National Aeronautics and Space Administration (NASA) and the National Geospatialintelligence Agency (NGA). The SRTMGL1 product can be downloaded at NASA's Earth Observing System Data and Information System (EOSDIS) website (http://reverb.echo. nasa.gov/ accessed on 1 July 2021). The World Geodetic System 84 (WGS84) ellipsoid and Earth Gravitational Model 1996 (EGM96) geoid were used in the SRTMGL data products [37]. Several studies have validated the STRMGL1 data and shown that the product has high accuracy in China [38,39]. The SRTMGL1 data in this study were used to divide small-scale basins.

Division of Small-Scale Basin
The slope runoff simulation algorithm proposed by Ocallaghan et al. [40] is widely used in extracting basin features based on DEM [41,42]. Based on the slope runoff simulation algorithm, the river network information of the GRB is extracted from the SRTMGL1. After the SRTMGL1 data are input into the model, the vector river network data of the basin are first generated by the slope runoff simulation algorithm. The basic principle is to simulate the movement path of real surface water flow to describe the catchment capacity of each pixel in the region. Then, the final river network of the region is the basin part whose cumulative flow is greater than the threshold. Using the regional hydrological station information as a constraint, the closest catchment location to these stations is searched for as the basin outlet in the calculated river network. Through the flow direction data calculated in the slope runoff simulation algorithm, all raster image elements controlled by each outlet are retraced one by one. According to the characteristics of pixel distribution, the small-scale basin data can be obtained by connecting the center points of the outermost grid in turn. In this study, we defined a basin an area of~600 km 2 as a small basin.
It is crucial to select appropriate accuracy evaluation indicators for evaluating the accuracy of small basin divisions. Among the river network accuracy evaluation methods, the match error of river system extraction is a common quantitative index, which is to superimpose the extracted river network with the real river network, calculate the fine polygon area produced by its superimposition, and then calculate its proportion in the total basin area. The calculation equation is shown as follows: D is the match error of basin system extraction, A i is the area of the finely divided polygons produced by the superposition of the two river networks, and S is the total area of the basin. The research shows that if D is small, it indicates that the river network is better consistent with the real river network. In general, when D is less than 0.02, the precision is high. When D is less than 0.03, the precision is reasonable. When D is more than 0.03, the precision is poor.

Simulation of High-Resolution Regional ET
ET is the critical link of the hydrological cycle and an essential part of water balance and energy balance [43]. Therefore, the accurate estimation of ET is very important for estimating the TWSC based on the water balance equation in this study.
In order to obtain high-precision ET in the GRB, we used the Penman-Monteith-Leuning Version 2 (PML_V2) model to simulate ET. The PML_V2 has the following features: (i) a relatively simple structure to estimate the coupling of ET and GPP while retaining reasonable biophysical meaning; (ii) includes the effects of CO 2 concentration changes on carbon assimilation and canopy conductance; and (iii) can be applied on regional and large scales easily because it only requires daily meteorological and remote sensing forcing data [44]. Furthermore, by comparing with the ground observation data of 95 global flux stations, this product has higher accuracy and is better than most available ET products [45][46][47].
Although PML_V2 performs well on a global scale, there is still room for improvement at regional scales. The MODIS and GLDAS data were used to produce PML_V2 ET data, but GLDAS data used less ground observation data in China [48,49]. Thus, in order to further improve the accuracy of the modeled ET in China, the input data were changed from GLDAS to CMFD in this study. CMFD has advantages of high precision and wide coverage, and it has been widely used in hydrological research [30,31]. Therefore, we used the PML_V2 model and CMFD dataset to estimate the regional ET for the GRB.

Simulation of Runoff Based on SWAT Model
Hydrological models are widely used to simulate the runoff process, and many scholars have employed some common models such as the Variable Infiltration Capacity (VIC) model [25,26,50], Soil and Water Assessment Tool (SWAT) model [51][52][53], and Identification of Hydrographs and Components from Rainfall, Evaporation, and Stream (IHACRES) model [50,51] to simulate the runoff in various river basins. However, the condition of each basin is different, especially the difference in climate, which makes the selection of hydrological models critical to the accuracy of runoff simulation [52]. In terms of this issue, Huang et al. [54] evaluated the performance of 9 hydrological models for 12 large-scale river basins around the world. In their conclusions, the SWAT model performs well in simulating runoff in the Yangtze River Basin. Some other studies in the literature have also shown that the SWAT model shows good results in the simulation of runoff in China compared with other models, e.g., the VIC model and Xinanjiang model [55][56][57]. As a physically based model, the SWAT model can simulate the short-term and long-term hydrological process, chemical process, erosion process, soil hydrological cycle, and crop production process of the basin rapidly [23,58]. Compared with other models, the SWAT model can simulate the response of runoff to LUCC in most basins in China more obviously and accurately [55,59]. The SWAT model considers the longitudinal movement of water flow in each small unit and the lateral exchange of water flow between each unit and each other, so it is suitable for the simulation of lateral runoff at small scales. Therefore, the SWAT model was used to estimate the outbound runoff at small-scale basins in this study.
The SWAT model of the GRB was established by using the DEM, soil type, land cover, precipitation, and meteorological data. CN2.mgt and 11 other parameters related to runoff were selected by referring to other studies [59]. SWAT-CUP 2012 was used to obtain the ranges of each parameter and sensitivity order through multiple iterations. The sequential uncertainty fitting (SUFI-2) algorithm was used to update parameter values. Daily precipitation and monthly runoff data in 2002-2011 and 2012-2018 were used for the calibration and validation, respectively. The applicability of the model was evaluated by the Nash-Sutcliffe Efficiency (NSE), coefficient of determination (R 2 ), and Percent Bias (PBIAS) indexes.

Estimation of TWSC Based on the Water Balance Equation
According to the law of conservation of mass, in the hydrological cycle, for any region and period, the difference between the inflow and the outflow of water must be equal to the change in the water storage, which is known as the basin water balance equation. Therefore, we calculated the terrestrial water flux (the slope of the change in terrestrial water storage over a given period) for the GRB using the basin water balance equation (Equation (2)) based on small-basin-scale precipitation, ET, and runoff. Detailed water balance components are shown in Table 1. The TWSC is obtained by integrating the land water flux.
dS/dt represents the TWSC for a specific period (t), P represents precipitation, ET represents evapotranspiration, and Q represents runoff.

Small-Scale Basin Division in the GRB
Based on the slope runoff simulation algorithm, the small-scale basin division of the GRB was completed. The small-scale basin division results of the GRB can be divided into two levels. The first level is the control layer of hydrological stations, each hydrological station controls a basin, and the only outlet in the basin is the corresponding hydrological station (as shown in Figure 2a). The second level is the area control layer, subdivided according to the topographic features and river direction based on the first level results. Finally, the GRB is divided into 172 sub-basins, with an average basin area of 591 km 2 (as shown in Figure 2b).
Because the real river network information of the GRB is difficult to obtain, our study used the river network data provided by Open Street Map (OSM) to calculate the match error of the extracted river network in the GRB. The results showed that the match error of the extracted river network in our study is 0.013, which meets the accuracy requirement.

Variation in Precipitation in the GRB
The spatial distribution of average annual precipitation in the GRB from 2002 to 2018 is shown in Figure 3a. The mean annual precipitation of the GRB during the study period ranges from 1500 mm to 1800 mm, with an average value of 1644 mm. The areas with large precipitation are mainly concentrated in the central and northern regions. The precipitation in the central and northern parts of the GRB is more than 1700 mm, while the average annual precipitation in other regions varies from 1500 mm to 1600 mm. Generally speaking, the average annual precipitation is quite different in spatial distribution. The precipitation in the sub-basins has a large inter-annual variation, and the maximum annual precipitation in the same sub-basin is 2 to 2.5 times that of the minimum annual precipitation.
For the monthly mean precipitation variation of 172 sub-basins in the GRB (Figure 4b), the higher precipitation in most sub-basins occurs in May and June, and the lower precipitation is in October. In May and June with high precipitation, the difference in precipitation in different sub-basins is more than 100 mm. In August, the difference in precipitation between different sub-basins increases further and the difference can reach more than 150 mm. The monthly mean precipitation of different sub-basins varies greatly in the same month. However, the changes in monthly precipitation in the 172 sub-basins of the GRB show a relatively consistent trend during the study period.

Variation in ET in the GRB
The spatial distribution of the average ET in the GRB from 2002 to 2018 is presented in Figure 3b. Figure 3b shows that the average annual ET of the GRB during the study period varies from 600 mm to 1000 mm, with an average of 780 mm. The regions with maximum ET are concentrated mostly in the middle and south, while the ET in the northwest is lower than those in other regions, and the spatial difference of annual average ET is obvious. Figure 4c,d show the variation in annual and monthly mean ET in the 172 sub-basins from 2002 to 2018. From Figure 4c, we notice that the ET in most sub-basins presents a relatively consistent fluctuation during the study period, but the annual ET values of different sub-basins in the same year are quite distinct, and the difference in annual ET among the sub-basins can reach 300 mm. Figure 4d shows that the monthly average ET of different sub-basins presents highly consistent fluctuations. The monthly ET increases from January, reaches a peak in July, and then begins to decrease. The average monthly ET of different sub-basins varies from 10 mm to 150 mm. Noting that precipitation and ET show different seasonal variations, the precipitation peaks from May to June, but the peak of ET occurs from July to August, and ET may be affected more by temperature in the GRB. Figure 5 shows the scatter plots of the simulated and measured monthly runoff data at Waizhou station during the parameter calibration period, with an overall good fitting accuracy and an R 2 of 0.88 between the simulated and measured values.     (Figure 4a,b) may cause the runoff changes during the same period. Figure 4f shows the change in monthly average runoff in 172 sub-basins from 2002 to 2018. Overall, it shows an upward trend in the first half of the year, reaches a peak in June, and shows a downward trend in the second half of the year. However, the monthly average runoff in the different sub-basins is quite different, and the difference can reach 100 mm. It should be noted that the seasonal runoff is similar to seasonal precipitation.

Spatial and Temporal Analysis of TWSC Variation in the GRB
Based on CMFD precipitation, simulated ET, and runoff, we obtained the TWSC estimates of the 172 sub-basins from 2002 to 2018 using the water balance equation. Figure 8 shows the spatial distribution of the average annual TWSC in the GRB from 2002 to 2018. It can be seen from Figure 8 that the TWSC is negative in some areas of the upper reaches, while it is positive in the middle and lower reaches of most sub-basins. One sub-basin in the northeast of the basin has an obvious deficit. Both the largest losses and the largest increases in the TWSC appear in the downstream area. In order to further analyze the variation in TWSC in the GRB, we calculated the spatial distribution of TWSC from 2002 to 2008, as shown in Figure 9. It can be seen that there is a significant loss of TWSC in the GRB in 2003. In 2003, all the TWSC in the GRB shows a negative value, and the loss decreases from the southeast to the northwest. In 2004, the TWSC in the upstream area of the GRB also presents a negative value in most sub-basins. to 2018, most sub-basins in GRB show an increase in TWSC. It should be noted that the interannual variation in TWSC is in good agreement with the interannual variation in precipitation, and the TWSC in the GRB may be more influenced by precipitation.  Figure 4h shows that in most sub-basins of the GRB, the interannual variation in TWSC is consistent. The monthly mean TWSC in most sub-basins increases from January to June, while the TWSC is in a state of loss from July to October. The difference in TWSC among different sub-basins is large in August but is little in other months. As described in Sections 4.2 and 4.4, the variations in precipitation and runoff in different sub-basins are also large in August, which may be the reason for the large difference in TWSC.

Comparison with Other TWSC Data
To verify the accuracy of the method in this study, we first compared the TWSC calculated in this study with two TWSC datasets in the whole GRB, considering the coarse resolution of GRACE. The first dataset was the CSRM product. The temporal gaps in CSRM were filled with datasets published by Zhong et al. [60]. With the second data, we used the method proposed by Humphrey et al. [61] to reconstruct the daily TWSA, and further calculated the TWSC of the natural month in the GRB.
Due to the low resolution of GRACE satellites, we evaluated the overall TWSC in the GRB (as shown in Figure 10). The RMSE for the TWSC calculated in this study compared with the CSRM from 2003 to 2018 is 37.1 mm/month. The RMSE between the TWSC calculated in this study and the reconstructed TWSC datasets further reduce to 21.9 mm/month, and the RMSE between CSRM and the reconstructed TWSC is 29.8 mm/month.

Comparison with In Situ Groundwater Level Data
We also used the in situ groundwater level data of the Ganzhou Hotel well site from January 2007 to December 2010, the Zhangjiawei well site from January 2013 to May 2015, and the Yichun well site from January 2009 to December 2012 to verify the TWSC estimates in this study. We compared the Ground Water Level Change (GWLC) with the TWSC estimate in the corresponding sub-basin indirectly, as shown in Figure 11. As a result, the in situ groundwater level data of the three stations and the TWSC of the corresponding sub-basin show a relatively similar variation. However, due to the lack of long time series in situ data, we only analyzed the groundwater level and the TWSC of the corresponding sub-basin of three stations with relatively short time series. In addition, the correlation coefficients (CC) of GWLC and TWSC for all three stations during the analysis period are greater than 0.4.

Comparison with Soil Moisture Storage Change (SMSC)
We adopted the CLDAS-V2.0 soil moisture storage (SMS) product to calculate the SMSC of the GRB. Then, SMSC was calculated as the difference between SMS in two consecutive months. We calculated the CC between TWSC and SMSC in 172 sub-basins over the common overlap period (April 2017-December 2018); the spatial distribution is shown in Figure 12. The CC between TWSC and SMSC for 172 sub-basins ranges from 0.26 to 0.88, and the median CC is 0.56. The CC of most sub-basins is greater than 0.5. The sub-basins with larger CC are mainly concentrated in the downstream area, while the sub-basins with smaller CC are mainly concentrated in the upstream area.
Two sub-basins with high CC (labeled (a) and (b)) and two sub-basins with minimum CC (labeled (c) and (d)) were selected for further analysis. It can be seen from Figure 12c,d that although the CC between TWSC and SMSC in these two sub-basins is small, their changes show a relatively good agreement. The most extreme values occurred in the same month. Overall, time series are similar in both data, and the main maxima and minima are found in a similar month. These three comparisons suggest that the TWSC calculated in our study is reliable.

Uncertainty
Uncertainties in ET and runoff can come from multiple sources: model misrepresentation or parameterization, forcing data, and other variables [62,63]. Due to the complexity of the model, it is difficult to estimate the uncertainty of ET and runoff caused by the uncertainty of the input precipitation forcing dataset. To quantify the uncertainty of the hydrological model output data, Koukoula et al. [64] evaluated water cycle components from 18 state-of-the-art water resources reanalysis (WRR) datasets derived from different hydrological models, meteorological forcing, and precipitation datasets. They found that the uncertainty in the estimation of the water cycle components from different products is mainly attributable to the differences in the schemes used by the different models, while different precipitation forcing datasets have less impacts on the precision of the WRR output data. Therefore, in this study, we assumed that each variable is independent and estimated the uncertainties of TWSC.
Uncertainties in the monthly estimates of TWSC were assessed by combining the individual uncertainties for precipitation, ET, and runoff with the propagation of errors [65] (Equation (3)).
where σ 2 P , σ 2 ET , and σ 2 Q represent the error of precipitation, ET, and runoff, respectively. Due to data limitations, we simply calculated the uncertainties of the TWSC for the whole sub-basins in the GRB. We used monthly precipitation measurements from three precipitation stations in the GRB from the CMA to evaluate the accuracy of the CMFD precipitation estimations. The uncertainty of ET was derived from the RMSE between PML_V2 ET and the ET estimations from the water balance method in the GRB [66]. We used monthly runoff measurements from Waizhou station to evaluate the accuracy of the runoff estimations. The σ P , σ ET , σ Q , and σ TWSC are ±17.3 mm/month, ±21.2 mm/month, ±28.7 mm/month, and ±39.6 mm/month, respectively. The uncertainty of TWSC estimated with this method is close to the RMSE between the mean TWSC of all sub-basins and TWSC from CSRM in the GRB.

Impact of ET and Runoff
In the process of ET simulation in small basins, we replaced the forcing data used in the PML_V2 model with the CMFD dataset, but the ET calculated by the PML_V2 model was globally calibrated. There was no region-specific calibration parameter and its accuracy may need further improvement in specific regions at small scales. Zhang et al. [44] showed that when using the PML_V2 model to simulate ET, due to the forcing data of Collection 6 MCD12 has spectral confusion between evergreen broadleaf forests (EBF) and evergreen needleleaf forests (ENF). It is difficult to distinguish several classes using MODIS coarse resolution data, which may lead to some potential uncertainty in ET estimation. Srivastava et al. [24] also showed that due to the MODIS algorithm not considering the effects of cloud cover and leaf shadow effects, MODIS-ET is more accurate in the dry season than in the flood season. The impact of ET by using satellite products also needs to be considered in future research.
We found that the CC between TWSC and SMSC in the sub-basins of the downstream region is generally larger, and the regions with smaller CC are mainly concentrated in the upstream region ( Figure 12). It can be seen from Figure 3a,c that the downstream region is also the area with large precipitation and runoff. We speculate that this phenomenon may be related to the accuracy of the runoff simulation. Muleta et al. [67] found that the sensitivity of SWAT model parameters and the accuracy of runoff simulation are significantly different between the periods of high and low precipitation. To address this issue, Levesque et al. [68] applied runoff data from periods of low precipitation to calibrate the SWAT model in a Canadian basin dominated by ice melt recharge and improved the simulation accuracy of the model. In the GRB, there are large spatial differences in precipitation between upstream and downstream regions. However, when we simulated runoff using the SWAT model, the model calibration was carried out for the whole basin due to limited in situ runoff data, which may lead to a higher accuracy of the runoff simulation in the downstream region compared to the upstream region. With more runoff data, we can determine the model parameters from upstream and downstream regions, and further improve the accuracy of the upstream runoff simulation.
Besides, when the SWAT model was applied to simulate the runoff process at a small basin, due to the lack of data about the groundwater process and aquifer system, especially the natural and artificial water storage systems closely related to hydrological conditions, such as dams and reservoirs, wetlands and ponds, the precision is constrained. Kumari et al. [69] also proved that when using the hydrological model to simulate runoff, the shape of a flow duration curve is mainly affected by reservoirs, land-use types, and upstream abstractions of water. The GRB is widely distributed, and several artificial water conservancy facilities such as reservoirs and dams play an important role in the hydrological cycle. If these data are coupled into the model, the accuracy of runoff simulation will be further improved.

Conclusions
In this study, we took the GRB as the study area and presented a new method for calculating TWSC at small-basin scales. First, we divided the GRB into 172 sub-basins based on the slope runoff simulation algorithm. Then, we simulated the ET and runoff of each sub-basin with the PML_V2 and SWAT. Finally, through the water balance equation, the TWSC of each sub-basin was obtained. We evaluated its effectiveness by comparing the estimated TWSC against the CSRM, in situ groundwater level data, and SMSC. We also analyzed the temporal and spatial variations in precipitation, ET, runoff, and TWSC in the GRB from 2002 to 2018. From the estimated results and spatio-temporal analyses, we can draw the following conclusions.
The results of this research demonstrate that compared with the in situ data, the simulated runoff is reliable, and the proposed method is capable of obtaining the TWSC at small scales. The TWSC calculated by our method agrees well with the CSRM, in situ groundwater level data, and SMSC. Compared with the CSRM and reconstructed natural monthly TWSC, the RMSE is 37.1 mm/month and 21.9 mm/month, respectively. By comparing with in situ groundwater level data, and SMSC, it shows higher CC. These results indicate that the proposed method makes sense and is beneficial.
The spatial and temporal variations in precipitation and runoff in the GRB are similar, with the minimum value in 2011 and the maximum value in 2016, and the high-value area is mainly concentrated in the downstream area. The annual ET changes gently, with the minimum value in 2007 and the maximum value in 2004. The TWSC fluctuates greatly in the study period, with the minimum value in 2003 and the maximum value in 2015, and the spatial distribution is uniform.
From this study, the accuracy of the evaluation method still needs to be further improved because of the lack of accurate and comparable TWSC observations. We used the GRACE product, the groundwater monitoring level observations, and soil moisture storage to evaluate the method's accuracy. However, this evaluation method is only qualitative and not quantitative. Nevertheless, we could conclude that in small basins, our method provides reliable data that are suitable for the study of small-basin TWSC. In future work, the quality of the TWSC could be improved by using more accurate runoff and ET simulations.  Data Availability Statement: The in situ data used in this study were mainly provided by the institution introduced in Section 2.2.1. The China meteorological forcing dataset (CMFD) is available at (http://data.tpdc.ac.cn/en/data/8028b944-daaa-4511-8769-965612652c49/; accessed on 13 August 2021). CMA Land Data Assimilation System Version 2.0 (CLDAS-V2.0) is available at (http://data.cma.cn/data/detail/dataCode/NAFP_CLDAS2.0_NRT/; accessed on 12 August 2021). The GRACE Mascon product from the Center for Space Research is available at (www2.csr.utexas. edu/grace/; accessed on 1 July 2021). The global 1 arc-second SRTM V3.0 dataset (SRTMGL1) is available at (https://lpdaac.usgs.gov/products/srtmgl1v003/; accessed on 1 July 2021). The PML_V2 evapotranspiration is available at (https://code.earthengine.google.com/e0453cf3e7a6e62513da409 89f29a029; accessed on 13 August 2021 from GEE). The soil map is available at (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/; accessed on 13 August 2021). The dataset used to fill the temporal gaps in the GRACE CSR Mascon product is available at (http://data.tpdc.ac.cn/en/data/71cf70ec-0858-499d-b7f2-63319e1087fc/; accessed on 13 August 2021).