Quantifying Multi-Source Uncertainties in GRACE-Based Estimates of Groundwater Storage Changes in Mainland China

Abstract

The Gravity Recovery and Climate Experiment (GRACE) satellites have been widely used to estimate groundwater storage (GWS) changes, yet their uncertainties related to the multi-source datasets used are rarely investigated. This study focuses on quantifying the uncertainties of GRACE GWS estimates in mainland China during 2003–2015, by generating a total of 3456 solutions from the combinations of multiple GRACE products and auxiliary datasets. The Bayesian model averaging (BMA) approach is used to derive the optimal estimates of GWS changes under an uncertainty framework. Ten river basins are further identified to analyze the estimated annual GWS trends and uncertainty magnitudes. On average, our results show that the BMA-estimated annual GWS trend in mainland China is −1.93 mm/yr, whereas its uncertainty reaches 4.50 mm/yr. Albeit the estimated annual GWS trends and uncertainties vary across river basins, we found that the high uncertainties of annual GWS trends are tied to the large differences between multiple GRACE data and soil moisture products used in the GWS solutions. These findings highlight the importance of paying more attention to the existence of multi-source uncertainties when using GRACE data to estimate GWS changes.

1. Introduction

Groundwater plays an indispensable role in human society and in natural systems [1]. Unreasonable use of groundwater resources has led to serious groundwater storage (GWS) depletion in many regions of the world, thereby threatening the safety and security of human water supplies [2,3,4,5]. It is thus essential to gain insight into the changing regularity of GWS so as to promote the sustainable use of groundwater resources. Traditional monitoring of GWS changes often relies on well observations. Even with high accuracy and reliability, it is difficult to obtain information on large-scale GWS changes due to the limited number of monitoring wells and high maintenance costs [6]. The Gravity Recovery and Climate Experiment (GRACE) satellites, launched in 2002, are capable of detecting monthly terrestrial water storage (TWS) changes by measuring the time-variable gravity fields of the Earth. Further deducting the changes of non-GWS components (e.g., surface water storage (SWS), soil moisture storage (SMS), snow water equivalent (SWE), canopy water storage (CWS), glacier water equivalent (GWE), etc.) from GRACE TWS, regional GWS changes can be calculated accordingly [7,8]. Up to now, the approach of GRACE-based TWS decomposition has been widely used to investigate large-scale GWS changes in various regions, such as the northern India [9], the Central Valley of California [10], and the North China Plain [11].

The accuracy and reliability of GRACE-based GWS estimation, however, are directly related to a range of errors [12]. Specifically, the raw GRACE Level-2 data products are released in a format of spherical harmonic (SH) solution. Due to the limited measurement accuracy of satellite instruments and the degree/order truncation (usually 60 or 90) of SH coefficients, GRACE data contain inherent measurement errors [13]. The north–south striping errors are also tied in GRACE TWS inversion because unwanted correlation exists in odd and even order SH coefficients, respectively [14]. When a low-pass filtering is used to suppress the noise of SH coefficients, it could induce the mutual signal leakages between grid cells, that is, the leakage error [15]. In further analysis of multiple GRACE and ancillary datasets with inconsistent spatiotemporal resolutions, rescaling the datasets to reconcile their resolutions, if required, could lead to the resampling errors [16,17]. Moreover, existing studies that focused on quantifying the annual trends or periodic changes of TWS and GWS, may vary in the estimated results owning to the different fitting approaches used, thereby introducing the fitting errors [18,19,20]. Currently, some correction approaches (e.g., the decorrelation filtering, scale factor and the forward modeling) have been developed to effectively suppress the errors of north–south striping and signal leakage in GRACE TWS inversion. The later released mass concentration (Mascon) algorithm can also well derive TWS changes from GRACE data by correcting the relevant errors [21]. There are three main approaches to performing Mascon solutions. The first is to use a priori information from geophysical models and remote sensing observations to constrain the solution within a Bayesian framework, represented by the JPL Mascon solution [22]. The second is based on finite order SH coefficients truncated by the spacing velocity or spacing acceleration, represented by the GSFC Mascon solution [23]. The third is the calculation on an equal-area geodesic grid by defining each mass tile as a finite truncated SH solutions without using external geophysical models or data, represented by the CSR Mascon solution [24].

Another source of errors in further estimating GWS with GRACE data is the uncertainty contribution of the non-GWS component data used in TWS decomposition. As the large-scale changes of SWS and SMS in most regions are difficult to obtain from site observations, previous studies related to GRACE GWS estimation tend to use the simulated data of non-GWS components from model simulations. One thing to note is that the outputs of different models for the same hydrological variable may differ due to inconsistent driving data, biased parameter values, and imperfect physical structure used in the modelling. This further leads to varying GWS estimates when different model outputs are used for GRACE TWS decomposition [25]. Many previous studies have confirmed the existence of multi-source errors regarding non-GWS components and their larger impacts on GRACE GWS estimation [26,27]. For example, by using the Global Land Data Assimilation System (GLDAS) and the WaterGAP Global Hydrology Model (WGHM) to estimate non-GWS component changes, respectively, Jin and Feng [28] confirmed that the maximum difference of annual GWS trends estimated by GRACE can reach 1.66 mm/yr in Asia, Europe, and North America, with the opposite estimates in South America. Yin et al. [29] likewise documented that the choice of different model data had a remarkable influence on the GWS estimation in Australia, that is, the annual GWS trend estimated by GRACE-GLDAS combination is 1.8 times larger than that estimated by GRACE-WGHM combination. Given the importance of multi-source errors regarding non-GWS components and their non-negligibility in GRACE TWS decomposition, this study intends to perform the GRACE-based GWS estimation on the mainland of China. The commonly used model data in existing GRACE GWS studies are collected to quantify the multi-source uncertainties of non-GWS components, and eventually reassess the GWS changes of mainland China under an uncertainty framework [25,30,31].

Previous studies of GRACE GWS estimation in mainland China mainly focused on analyzing the spatiotemporal patterns and potential drivers of GWS changes [32,33,34] by using models to estimate non-GWS components, they often choose to calculate the average of multiple model outputs or validate with site observations for optimal selection, with the purpose of enhancing the credibility of GRACE GWS estimation. For example, Feng et al. [35] chose to decompose GRACE TWS by deducting the averaged non-GWS component changes simulated from four GLDAS models and estimated the GWS changes in three main aquifers of China. Zhang et al. [36] compared the SMS changes simulated by four GLDAS models with site observations and finally selected the best-fitting one to derive GWS changes from GRACE data in North China Plain. Huang et al. [37] adopted well-observed groundwater level data to validate the GRACE-estimated GWS changes and then analyzed the groundwater budget in karst area of southwest China. Zhang et al. [38] used the simulated GWS changes from well-established hydrological model as a validation of GRACE GWS estimates in Haihe River basin of China and confirmed the feasibility of satellite gravity technology in estimating GWS changes. Yi et al. [39] estimated the GWS changes in 10 major basins of China from 2003 to 2014 and likewise chose to average the outputs of non-GWS components from GLDAS models in GRACE TWS decomposition. Although these studies have attempted to estimate regional GWS changes using GRACE data throughout China, few or no studies have provided a holistic analysis on the multi-source uncertainties caused by different modelling data of non-GWS components used in GRACE GWS estimates. GRACE does have the capacity to provide a “big picture” view of large-scale GWS changes, but the existence of multi-source errors regarding non-GWS components may influence the robustness of GRACE-estimated GWS changes [40].

In this study, six GRACE datasets (three Mascon solutions and three SH solutions), three land surface model products (CLSM, VIC, and NOAH models from GLDAS-2.1), two global hydrological model products (WGHM v2.2d and PCR-GLOBWB 2), one reanalysis data product (ERA5-land), and two glacier datasets (degree-day glacier dataset and global glacier mass loss dataset), are collected to obtain a total of 3456 possible estimates of GWS changes in mainland China. By using the statistical metrics to quantify the differences of TWS and non-GWS component changes estimated by multi-source datasets, their uncertainty contributions to GRACE-estimated GWS changes are analyzed subsequently. Moreover, the Bayesian model averaging (BMA) ensemble approach is used to investigate the annual trends and posteriori information of GWS changes at the basin scale. Compared to previous studies, this study highlights the impacts of different model data choices on GRACE GWS estimates and reassesses the GWS changes in mainland China, which is expected to provide a scientific basis for further optimization of national groundwater resources management.

2. Study Area and Datasets

2.1. Study Area

China has a vast territory with diverse topography and climate, resulting in a complex hydrometeorological regimes over land. In terms of the topography and identified watersheds, the mainland of China can be divided into 10 major river basins (Figure 1): Songhua River basin (SRB), Liao River basin (LRB), Haihe River basin (HRB), Yellow River basin (YRB), Huaihe River basin (HHB), Yangtze River basin (YZRB), Pearl River basin (PRB), Southeast Rivers basin (SEB), Southwest Rivers basin (SWB), and Continental asin (CB) [41]. The eastern and southern coastal river basins receive abundant precipitation recharge with annual amounts above 1200 mm, whereas the northwest and interior basins are relatively dry with annual amounts below 200 mm. A few western river basins (e.g., SWB and CB) can receive water recharge from melting mountain glaciers. Similar to the pattern of precipitation, actual evapotranspiration in 10 river basins decreases gradually from southeast to northwest [42]. The uneven spatial distribution of hydrometeorological features has resulted in abundant water resources (accounting for ~80% of the national total) in southern China but scarce (accounting for ~20%) in northern China [43].

Groundwater is a vital source of human water use in mainland China. More than 100 km³ of groundwater resources are extracted annually throughout the nation, primarily for agricultural irrigation and residential water consumption. Long-term overexploitation of groundwater resources has led to a significant decline of groundwater level in some river basins [44]. For instance, the GWS in plain area of HRB was found to show a significantly decreasing trend in a rate of −32.0~−17.0 mm/yr during 2003–2018 [11]. Some localized areas in SRB, LRB, and YRB also experienced notable GWS depletion, which was confirmed to be related to the groundwater-fed agricultural irrigation and anthropogenic revegetation [35,45]. The excessive GWS depletion has caused a series of environmental problems (e.g., river blanking, land subsidence, and soil salinization), posing potential threats to economic development and food security. Hence, monitoring and understanding the changes of GWS in each basin is of great importance to identify areas where groundwater is over-exploited and aid in the development of sustainable water management strategies in China.

2.2. Datasets

2.2.1. GRACE Data

GRACE satellites, a joint mission of the National Aeronautics and Space Agency and the German Aerospace Center, were launched on 17 March 2002. Their major mission is to monitor the mass redistribution of the Earth’s surface by remotely sensing the time-variable gravity field of the Earth [46]. Changes in TWS induce mass redistribution and, therefore, can be estimated from GRACE data after removing atmospheric, ocean, and other time-variable gravity information. Initially, three official organizations, that is, Center for Space Research (CSR), GeoForschungs Zentrum (GFZ), and the Jet Propulsion Laboratory (JPL), are responsible for processing GRACE data into 0, 1A, 1B, and 2 levels to meet the requirements of various end-users [8]. To facilitate the investigation of TWS changes, the release of GRACE Tellus Level-3 data becomes a more convenient choice to provide the global monthly TWS changes with a spatial resolution of 1° × 1°. Another type of GRACE Level-3 product is the Mascon solution that claims to require no further error correction. Compared to SH processing, GRACE Mascon data show lower uncertainty and higher field observation correlation [24,47]. Currently, three types of Mascon products are generally used in existing studies, that is, CSR RL06 Mascon, JPL RL06 Mascon, and GSFC RL06 Mascon [48]. The spatial resolution of CSR RL06 Mascon data is 0.25° × 0.25°, whereas the resolutions of JPL and GSFC RL06 Mascon data are both 0.5° × 0.5°. In this study, multiple estimates of 2003–2015 monthly TWS changes in mainland China are achieved from three GRACE Tellus Level-3 products released by CSR, JPL, and GFZ, as well as three GRACE Mascon products released by CSR, JPL, and GSFC. To unify the spatial resolution of TWS data for subsequent analysis, we resampled all GRACE TWS data to 1° × 1° resolution using the nearest neighbor interpolation algorithm. Doing so may alter the spectral properties of raw TWS data to some extent, but is an essential process for performing subsequent TWS decomposition and uncertainty analysis.

2.2.2. Global Models and Reanalysis Data

The outputs of two global hydrological models (i.e., PCR-GLOBWB 2 [49] and WGHM [50]), three land surface models (i.e., CLSM-F2.5 [51], NOAH-3.6 [52], and VIC-4.1.2 [53] from GLDAS-2.1), and the ERA5-land [54] reanalysis product, are used to estimate monthly non-GWS component changes in mainland China during 2003–2015. A summary of the spatial resolution of each model and reanalysis product, as well as their estimated non-GWS components is listed in Table 1. Given the inconsistency of spatial resolution between datasets, all non-GWS data estimated by global models and reanalysis products are uniformly rescaled to 1° × 1° grid cells by using the nearest neighbor interpolation algorithm, similar to the processing of GRACE TWS data. In general, the global hydrological models quantify the human utilization of groundwater and surface water, as well as the water flow and storage, to assess the change of terrestrial water cycle across the globe [55,56]. The land surface models aim to generate the best land surface simulations by assimilating satellite and ground observation data through advanced land surface modeling and data assimilation techniques [57]. They are currently applied to investigate climate variability and hydrological cycles worldwide [58,59]. ERA5-Land is a reanalysis dataset to describe the evolution of water and energy cycles over land from 1950 to the present. Compared to previous ERA5, ERA5-Land combines model data with more ground observations, eventually forming a more complete and consistent dataset at the global scale [60,61].

2.2.3. Glacier Datasets

Asia’s mountain glaciers have been proven to be rapidly melting and shrinking under climate warming, while contributing to regional TWS changes [62]. Here, two glacier datasets, one from the degree-day model glacier dataset [63] and the other from the global glacier mass loss dataset at the beginning of the 21st century [64], are used to estimate the monthly GWE changes in mainland China during 2003–2015. The degree-day model is a function of precipitation and temperature that describes the changes of glacier mass balance with altitude by simulating the hydrological processes of precipitation, snow accumulation, and snow and ice melting and refreezing [65]. This model has been widely used for estimating GWE changes globally [66,67]. More details regarding the forcing data and parameter assignment of the degree-day model we used can be found in Table S1 of the Supplementary Materials. Another dataset of global glacier mass loss uses high-precision digital elevation models to estimate global glacier mass balance. Its reliability and uncertainty estimates can be accurate to the scale of individual glaciers on a global scale. The raw datasets of glacier mass changes are represented as mass unit at 0.1° × 0.1° spatial resolutions. We further aggregate the raw datasets in a summation mode to the spatial resolution of 1° × 1° and then translate into GWE changes in a unit of equivalent water height by combining with the areas of grid cells.

3. Methods

3.1. GRACE-Based Estimation of GWS Changes

Theoretically, TWS comprises SWS, SMS, SWE, CWS, GWE, and GWS [68]. Obtaining the GWS is as simple as subtracting other water storage components from the TWS via the following formula:

GWSA = TWSA − SWSA − SMSA − SWEA − CWSA − GWEA(1)

where letter A represents the anomaly values that are calculated by subtracting the long-term mean of the corresponding TWS or individual storage component data from 2003 to 2015. Given that the data of six TWSA, two SWSA, six SMSA, six SWEA, four CWSA, and two GWEA, are used in this study, a total of 3456 (=6 × 2 × 6 × 6 × 4 × 2) GWSA solutions are obtained by their combined operation.

3.2. Performance Metrics for Quantifying the Uncertainties of GWS Estimation

To further quantify the discrepancies of TWSA and its component estimates from multi-source datasets, as well as their uncertainty contributions to GWS estimates, we employed three performance metrics, namely, standard deviation (SD), coefficient of variation (CV), and component contribution ratio (CCR):

(2)$\mathrm{S}\mathrm{D}=\sqrt{\frac{1}{N-1}\underset{n=1}{\sum ^N}{(x_n-\bar{x})}^2}$

(3)$\mathrm{C}\mathrm{V}=\frac{|\bar{x}|}{\mathrm{S}{\mathrm{D}}_{\mathrm{i}}}\times 100\%$

(4)$\mathrm{C}\mathrm{C}\mathrm{R}=\frac{\mathrm{S}{\mathrm{D}}_{\mathrm{i}}}{{\mathrm{S}\mathrm{D}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\times 100\%$

Note that, as the purpose of this study is to analyze the uncertainty contribution of different data choice on the estimation of annual GWSA trend, $x_n$ and $\bar{x}$ hereby are one of the annual trends and their arithmetic mean for the total number of $N$ estimates of TWS or individual storage component, respectively. ${\mathrm{S}\mathrm{D}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ is the total SD value equaling to the sum of $\mathrm{S}{\mathrm{D}}_{\mathrm{i}}$ values for the $i^{th}$ set including TWSA and all components. Simply put, SD is an indicator to measure how large the uncertainties in estimated annual trends of TWSA or individual storage components. CV measures the variability of estimated annual trends relative to their arithmetic mean value. As GWSA is linearly related to TWSA and non-GWSA components in Equation (1), the SD value of GWSA is not equal to the sum of the SD values of TWSA and non-GWSA, but rather a square root relationship. Hence CCR here reflects the normalized contribution ratios of multi-source uncertainties quantified by SD values to the estimated annual GWSA trends.

3.3. Merging GWS Changes Based on BMA

After obtaining GWSA ensemble solutions, the BMA approach is used to adjust the weights of each GWSA member with the aim of obtaining a statistically optimal estimation. The reference model is used as an estimate of the prior probability in BMA modelling by empirically selecting an estimated GWSA with relatively small uncertainty, that is, the GWSA estimation by subtracting the WGHM-simulated non-GWSA components (including SWSA, SMSA, SWEA, and CWSA) and the averaged GWEA of two glacier datasets from the TWSA estimation of GRACE CSR mascon data. We set the 2003–2011 timespan as the training period for BMA model calibration and the 2012–2015 timespan as the validation period for BMA merging. The weight ratios of GWSA ensemble are adjusted at each grid cell in the study area. Based on the derived weights, the optimal GWSA solution is obtained by multi-source fusing all GWSA solutions from 2003 to 2015. We calculate the annul trends of all GWSA solutions and BMA estimation using least squares harmonic analysis and further investigate their magnitudes and uncertainties in combination with the probability density function (PDF) statistics across 10 river basins. The specific steps of BMA-merging are described as follows [69]:

Step 1 assumes that $y$ is the ensemble of 3456 GWSA solutions. ${f=f}_1,\cdots ,f_k$ represents the estimated data of $k$ members (in this study, $k$ = 3456). $y^T$ is the GWSA reference model. The optimal solution of GWSA ensemble is unknown. Hence, transforming the probability distribution $p\left(\left.f_k\right|y^T\right)$ to a normal distribution, with the mean as a simple linear function $a_i+b_i{f}_k$ of the original model data (variance is expressed as ${\sigma }_k^2$), as shown in Equation (5):

(5)$\left.y\right|f_1,\cdots ,f_k,y^T~N\left(a_i+b_i{f}_k,{\sigma }_k^2\right)$

where $a_i$ and $b_i$ are the bias correction terms obtained through linear regression.

Step 2 is to obtain the weights of candidate models by a probability-weighted average, as shown in the following Equation (6):

(6)$p\left[\left.y\right|\left(f_1,\cdots ,f_k,y^T\right)\right]=\sum _{k=1}^{k}p\left(\left.f_k\right|y^T\right)p_k\left(\left.y\right|\left(f_k,y^T\right)\right)$

The posterior probability $p\left(\left.f_k\right|y^T\right)$ of estimated values for each GWSA solution in Equation (6) is the weight ratio ${\omega }_k$ of $k$-th GWSA member. $p_k\left(\left.y\right|\left(f_k,y^T\right)\right)$ is the conditional probability density function associated with the GWSA solution $f_k$. $p\left(f_k\right)$ is the prior probability of optimal estimation that can be taken as a uniform distribution in the absence of specific information, that is $p\left(f_k\right)=1/k$.

Step 3 is to calculate the expected value of multi-source BMA merging, as shown in Equation (7):

(7)$E\left(\left.y\right|\left(f_1,\cdots ,f_k,y^T\right)\right)=\sum _{k=1}^{k}p\left(\left.f_k\right|y^T\right)·E\left[g\left(\left.y\right|f_k,{\sigma }_k^2\right)\right]=\sum _{k=1}^k{\omega }_k{f}_k$

where $g\left(\left.y\right|f_k,{\sigma }_k^2\right)$ represents a normal distribution regarding $f_k and {\sigma }_k^2$. In this study, ${\omega }_k$ and ${\sigma }_k^2$ are calculated by the Markov Chain Monte Carlo algorithm based on its advantage in applying to samples with arbitrary probability distributions [70]. Since the initial probability distributions of GWSA ensemble are different, the state transfer matrix of the Markov chain model can be made independent of the initial state probability distribution after convergence, so that the probability distributions of the final state can converge to the same stable probability distribution.

4. Results and Discussion

4.1. Uncertainty Estimation of TWS and Non-GWS Component Changes

Figure 2 uses SD values to show the differences of the annual trends of TWSA and its non-GWSA components (including SWSA, SMSA, CWSA, SWEA, GWEA) estimated by multiple datasets in mainland China during 2003–2015. The spatially averaged SD values, arithmetic mean and CV values of the annual trends for individual water storages at the basin scales are summarized in Table 1. The spatial distributions of annual trends for TWSA and non-GWSA components estimated by individual datasets are shown in Supplementary Figures S1–S6 in the Supplementary Materials.

As shown in Figure 2a, the notable discrepancies of annual TWSA trend estimates are found among six GRACE datasets, especially in regions with large annual trend magnitudes, such as the central–eastern, northwestern, and southwestern regions of China.

This is mainly because the Mascon algorithm can better constrain the leakage errors in regions with large TWS changes compared to the SH solutions [24,47,71], such that the magnitude of annual TWSA trend estimated by Mascon data is relatively larger in these regions than the estimated results of SH products (Supplementary Figure S1) [72]. At the basin scale (Table 1), the uncertainty of estimated annual TWSA trends is highest in SWB with an SD value of 5.18 mm/yr and lowest in PRB with an SD of only 1.55 mm/yr. The estimated annual TWSA trends in other river basins exhibit moderate uncertainties, with SD values varying in the range of 1.72–4.46 mm/yr. However, it is noteworthy that if comparing the percentage of SD value to the mean annual trend (i.e., the CV value showing the level of dispersion around the mean), the river basin where uncertainty has the greatest impact on the estimation of annual TWSA trend is SRB (CV = 188%). The mean annual trend of TWSA in SRB is only 0.91 mm/yr, whereas the SD value is as high as 1.72 mm/yr (Table 1). This indicates that the difference of estimated annual TWSA trends using different GRACE data is objectively existing in SRB. A similar case can be found in CB and YZRB where CV values reach 126% and 85%, respectively, likewise suggesting the need for caution when using GRACE data to estimate TWS changes in these river basins [73,74,75,76].

Conversely, the CV value in HRB is only 33%, which is mainly because the significant decreasing trend of TWS (−13.36 mm/yr) and the relatively consistent estimates of GRACE data with the SD value equal to 4.46 mm/yr (Figure S1). The LRB, YRB, and PRB all have the similar pattern that the SD values account for less than half of the mean annual TWSA trends (Table 2).

Differences are observed in the annual trends of non-GWS components estimated by multi-source datasets (Figure 2b–f). The SD values estimated for annual SWSA trends are larger in some scattered grid cells located in central and eastern China, with values ranging from 1.05 to 2.86 mm/yr. The underlying cause for the large SD values is that the number of reservoirs and lakes considered by different models in simulating SWS changes is not exactly the same. For example, the PCR-GLOBWB simulates SWS changes of Miyun Reservoir (116.88°E, 40.47°N) and Xidayang Reservoir (114.83°E, 38.84°N), but they are not simulated by the WGHM, leading to significant differences in estimating SWSA for the correspinding grid cells. Due to the largest number of reservoirs/lakes simulated by different models in YRB and the largest variability of SWS in this region, the largest SD value (2.86 mm/yr) of estimated annual SWSA trends occurs in YRB, corresponding to the mean annual trend of 2.66 mm/yr and the CV value of 107% (Table 1). Furthermore, even if all models simulate SWS changes for the same reservoir, the simulated results may differ due to the differences of model structure and parameterization schemes. Take the Three Gorges Reservoir (111.08°E, 30.84°N) as an example, the SWS simulated by WGHM increased in an rate of 4.9 mm/yr, while that simulated by PCR-GLOBWB showed a larger increasing trend (21.1 mm/yr) during 2003–2015 (Figure S2). Such a discrepancy somewhat reflects the shortcoming and inconsistency of current global hydrological models in simulating SWS changes [49,50,77].

GWEA and SWEA show a similar pattern to SWSA, with larger SD values of annual trends appearing in some scattered grid cells (Figure 2c,d). For GWEA, the annual trends estimated by global glacier dataset reveal the loss of glacier mass in CB and SWB, whereas those estimated by degree-day model show the slight increasing trends of glacier mass in these areas (Figure S4). This may be related to the uncertainties of parameter values in degree-day model [64,78,79,80,81]. Among 10 basins, the SD values of multi-source datasets in estimating annual GWEA and SWEA trends are both highest in SWB, with SD values of 0.65 and 1.08 mm/yr, corresponding to CV values of 103% and 275%, respectively. However, the estimated annual SWEA and GWEA trends in CB have higher CV values, 1161% and 210%, respectively, suggesting the requirement of carefully selecting models to assess SWE and GWE changes in CB. Further, as the magnitude of annual CWSA trends estimated by all models is so small (−1.1 × 10⁻³ to 7.8 × 10⁻⁴ mm/yr nationally), the magnitude of further calculated SD values is also small (6.4 × 10⁻⁴ to 1.4 × 10⁻³ mm/yr nationally) (Figure 2e). Nevertheless, the comparisons of annual CWSA estimates among individual models still reveal relatively similar patterns simulated by CLSM and VIC, but with significant differences compared to the simulations of NOAH and WGHM (Figure S5). NOAH clearly reflects decreasing and increasing trends of CWSA in southern and northeastern China, respectively, whereas the estimation of WGHM is not clear in spatial variability and even shows some randomness among grid cells (Figure S5).

The SD values regarding annual SMSA trends are significantly larger in central, western, and northwestern China (Figure 2f), mainly attributed to the overestimation of annual SMSA trends by ERA5-land (Figure S6). Due to the imperfect parameterization scheme of soil freeze–thaw processes, ERA5-land has difficulty in capturing the increased soil moisture caused by snowmelt in winter and early spring in the regions of northern HHB and HRB, thereby resulting in an overestimation of the decreasing trends of SMSA [82,83]. Compared to other models, ERA5-land simulates permafrost changes other than liquid soil water in SWB [84]. This may partly explain the relatively high trend of increasing SMS estimated by ERA5-land in this region compared to other models (Figure S6). On the basis of the statistics in Table 2, the largest dispersion of annual SMSA trends estimated by multi-source datasets appears in HHB, with mean annual trend and SD value of −4.41 and 4.70 mm/yr, respectively. The CV value in this river basin is equal to 106%. Other river basins all have lower uncertainty in multi-source estimated annual SMSA trends than HHB, with mean annual SMSA trends ranging from −2.23 to 2.63 mm/yr and SD values from 0.79 to 3.29 mm/yr (Table 2). Except for SEB, the CV values of annual SMSA trends are higher than the mean annual SMSA trends in the remaining river basins, indicating that the selection of which model to estimate SMS changes in these regions should be carefully considered. Such a case is particularly evident in LRB where multiple models estimate positive or negative annual SMSA trends but with arithmetic means so close to zero (0.06 mm/yr), thereby resulting in large SD and CV values of 2.37 mm/yr and 4037%, respectively (Figure 2f and Table 2).

4.2. Estimation of Multi-Source Uncertainties in GRACE-Based GWS Changes

Figure 3 depicts the SD values of estimated annual GWSA trends from 3456 combinations in mainland China during 2003–2015, and the basin-scale statistics of multi-source uncertainty contributions of TWSA and non-GWSA components to annual GWSA trend estimates. As seen, the uncertainties of annual GWSA trends expressed by SD values do not correspond well to the arithmetic mean of annual GWSA trends in terms of the spatial pattern (Figure S7), due to the complicated influences of uncertainty contributions from TWSA and non-GWSA components. The regions with significantly large uncertainties are mainly located in the central and western regions of mainland China, indicating the existence of high uncertainties in estimating GWS changes from GRACE data (Figure 3). The uncertainties of estimated annual GWSA trends are less than 5 mm/yr in the vast majority of grid cells. However, there are some scattered grid cells with extremely high uncertainties in annual GWSA trend estimation, which is directly related to the large uncertainties of annual trend estimates for SWSA, GWEA, and SWEA (Figure 2b–d). To sum up, the spatially mean uncertainty of annual GWSA trend estimated by multi-source datasets across mainland China is 4.50 mm/yr and the uncertainties of 10 river basins ranges from 2.41 mm/yr for PRB to 8.41 mm/yr for HHB (Figure 3).

The further quantification of uncertainty contribution percentage shows that the annual trend uncertainties of TWSA and SMSA contributes most to the uncertainty of GWSA in 8 out of 10 river basins (Figure 3). Their averaged CCR values for uncertainty contribution reach 42% and 36%, respectively, and the total contribution of both is close to 80% in 8 river basins. Moreover, SWSA has largest contribution to the uncertainty of annual GWSA trend estimation in YRB with a CCR value of 36%, followed by TWSA and SMSA with both CCR values of 31%. SWSA ranks 2nd in SEB following the uncertainty contribution of TWSA to GWSA (CCR values of 44% for TWSA and 38% for SWSA). The uncertainty contribution of SWSA to GWSA is equally significant in HHB, ranking 3rd with a CCR value of 24%. Except in SWB and CB, the CCR values of SWEA and GWEA in the remaining river basins are relatively low (average values of 4% and 1%, respectively), and can be almost ignored. In SWB, apart from the primary uncertainty contribution from TWSA (CCR = 49%) and SMSA (CCR = 26%), the total contribution of SWEA and GWEA to the uncertainty of annual GWSA trend estimates exceeds 20%. The level of uncertainty contribution from SWEA (CCR = 6%) and GWEA (CCR = 2%) in CB is lower than in SWB, but can also account for a 10th of the total CCR in this region. These indicate that the impacts of SWEA and GWEA uncertainties on GWSA estimation is objective presence and cannot be ignored in SWB and CB [35,85]. On average in the entire mainland China, the uncertainty contribution of TWSA to GWSA estimation is always the largest, followed by SMSA and SWSA, which correspond to the CCR values of 42%, 36%, and 17%, respectively (Figure 3).

4.3. BMA-Based Estimation of GWS Changes under Uncertainty Framework

Figure 4 shows the spatial pattern of 2003–2015 annual GWSA trends in mainland China estimated by merging multi-source annual trend solutions using BMA algorithm, The PDF fitting of estimated annual GWSA trends and the optimal estimation of BMA- merging in each of the 10 river basins are plotted together. As seen, the annual GWSA trends in north, northwest, and southwest mainland China decrease significantly, whereas those in central and west mainland China increase obviously. The annual trends of GWSA in other regions are moderate relatively (Figure 4). Such a spatial pattern of GWS changes is quite consistent with the findings of Yi et al. [39] and Yin et al. [33]. However, if comparing the PDF-fittings of ensembles for mean annual GWSA trends in each basin, two curve types are found to confirm the existence of uncertainties in multisource-based estimation of GWS changes in mainland China. For the first type, the fitted PDF curves of annual GWSA trends for CB, SRB, YZRB, PRB, and SEB show only one peak feature with steep kurtosis, suggesting that the multisource-estimated annual GWSA trends in these river basins are in relatively good agreement with low uncertainty. There is a 95% probability that the annual GWSA trends for CB and SRB are in the −4 to 1 mm/yr interval, a 95% probability of being in the 0 to 5 mm/yr interval for YZRB and PRB, and an 95% probability of being in the −4 to 1 mm/yr interval for SEB due to a slightly lower concentration of ensemble estimation (Figure 4). The PDF curves of annual GWSA trends in SWB and LRB also exhibit one-peak feature but with slightly flat kurtosis, which means the larger uncertainty of estimating GWS changes in these river basins than aforementioned basins. The estimated annual GWSA trend of SWB has a 95% probability of being in a −14 to 0 mm/yr interval, whereas that in YRB has a 95% probability of being in the value interval of −12 to 0 mm/yr.

Regarding the second type of PDF-fitting, the fitted curves of annual GWSA trends show a multi-peak feature, that is, double-peaks in YRB, and triple-peaks in HRB and SEB (Figure 4). It is noted that more peak features mean flatter kurtoses and larger value intervals with the same conditional probability for equal sample sizes in this study. YRB has slightly flat peaks with a relatively long interval between two peaks, suggesting the large uncertainty in estimating annual GWSA trend in this area. The distribution of estimated annual GWSA trends in YRB is quite dispersed with a 95% probability of occurrence in the range of −14 to −1 mm/yr. The triple-peaks of PDF curves appearing in HRB and HHB imply the larger uncertainty of GWSA estimation. For HRB, the highest peak value corresponds to an annual GWSA trend of about −12 mm/yr, but there is a 95% probability that it will vary from −20 to −4 mm/yr. For HHB, however, the highest peak corresponds to an annual GWSA trend value of about −3 mm/yr, but a 95% probability of varying from −8 to 12 mm/yr.

By comparing the BMA-estimated annual GWSA trends of this study with previous ones, our study shows a decreasing trend of GWSA (−11.54 mm/yr) in HRB during 2003–2015, smaller than previous estimate of −22.00 mm/yr (2003–2010) by Feng et al. [35], but similar with Cao et al. [86] −14.67 mm/yr (1960–2008) and Yin et al. [33] −10.30 mm/yr (2002–2016). The annual GWSA trend of YRB estimated by BMA-merging is −7.38 mm/yr, comparable to the estimation of Xie et al. [87] (−6.50 mm/yr) during 2005–2014, but differing greatly to the estimation of Zhang et al. [88] (−3 mm/yr) during 2005–2013. Regarding the annual GWSA trends of LRB and HHB, the estimated annual GWSA trend of LRB in this study is −6.46 mm/yr that is significantly larger than the finding of Zhong et al. [89] (−3.16 mm/yr during 2005–2011). However, the estimated annual GWSA trend in HHB is −0.52 mm/yr, smaller than the −1.14 mm/yr estimated by Su et al. [90] from 2005 to 2012. For the annual GWSA trend of SWB, the BMA output of this study is −5.95 mm/yr, which is contrary to the result (5.59 Gt/yr) estimated by Zou et al. [91] in the Tibetan Plateau during 2013–2016. The discrepancies are likely related to the differences of datasets, methodology, and the size of study area. Overall, the BMA-estimated annual GWSA trends shown in Figure 4 are approximately the same as previous studies in terms of spatial pattern, but differs in the magnitude of annual trends in some localized areas.

5. Conclusions

By collecting the multi-source datasets of TWSA and non-GWSA components from GRACE data, global models, and reanalysis products, this study generated an ensemble estimation of GWSA totaling 3456 members in mainland China from 2003 to 2015. The BMA algorithm was used to derive the optimal estimates of GWS changes for the entire mainland China, and the mean annual trends of GWSA in 10 major river basins were further explored under an uncertainty framework. The main objective of this study is to investigate how large the uncertainties in the estimation of GWS changes caused by different choices of datasets in GRACE TWS decomposition, in particular for the calculated annual GWSA trends.

The results show that TWSA and SMSA contributes most to the uncertainties of GWSA estimates in 8 out of the 10 river basins in mainland China, whereas TWSA and SWSA have the largest contribution to the uncertainties of estimated GWSA in the remaining two river basins (i.e., YRB and SEB). The BMA-based estimation of GWSA for the entire mainland China shows a slightly decreasing trend in a rate of −1.93 mm/yr, but the corresponding uncertainty is as high as 4.50 mm/yr, arising from the different choices of multi-source datasets used in GRACE TWSA decomposition. The findings should be valued because the uncertainty magnitude is already larger than the estimated annual GWSA trend value, potentially leading to the low robustness of estimated GWS changes from GRACE data. Albeit the quantitative relationship between estimated annual GWSA trend and its uncertainty may vary in various river basins and grid cells, the differences of estimated GWS changes cannot be ignored when decomposing GRACE TWS using different datasets.

To sum up, our findings highlight the importance of accounting for multi-source uncertainties caused by dataset selection when estimating GWS changes using GRACE data. There are still considerable discrepancies in current datasets of TWS and storage components published by various official institutions. It is expected that these discrepancies can be reduced in the coming future with the continuously improved models and algorithms, thereby allowing for better estimates of GWS changes by using GRACE data.

Author Contributions

Conceptualization, Y.P.; methodology, validation, formal analysis, data curation, Q.L.; visualization, writing—original draft preparation, Q.L. and C.Z.; writing—review and editing, Q.L., Y.P., C.Z. and H.G.; supervision, C.Z., Y.P. and H.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data for GRACE Tellus Level-3 are available at https://podaac.jpl.nasa.gov/dataset, CSR RL06 mascon data can be found at https://www2.csr.utexas.edu/grace/RL06_mascons.html, JPL RL06 mascon data can be found at https://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons, GSFC RL06 mascon data can be found at https://earth.gsfc.nasa.gov/geo/data/grace-mascons. Other data can be found in the references.

Acknowledgments

The authors express their deep gratitude to the funding agency for supporting this research.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Spatial distribution of 10 major river basins in mainland China. The text in the lower left corner describes the full names, abbreviations and corresponding areas (unit: 104 km2) of 10 river basins.

Figure 2. Spatial distributions of differences in the annual trends of TWSA and non-GWSA components estimated by multiple datasets in mainland China from 2003 to 2015, expressed as standard deviation (SD) values. (a) TWSA; (b) SWSA; (c) GWEA; (d) SWEA; (e) CWSA; (f) SMSA.

Figure 3. Spatial distribution of SD values for annual GWSA trend estimates and the multi-source uncertainty contribution of annual GWSA trends from TWSA and non-GWSA components in mainland China. The percentage contribution of TWSA and non-GWSA components to the uncertainty of GWSA in each of the 10 river basins and mainland China is separately counted and illustrated in a ring chart. CWSA is not shown in the ring chart as its contribution to the uncertainty of annual GWSA trend estimate is too small (≤1%).

Figure 4. Spatial distribution of the 2003–2015 annual GWSA trends (unit: mm/yr) estimated by the combination of GRACE multi-source ensemble and BMA-merging algorithm in mainland China. The PDF of estimated annual GWSA trends in each of the 10 river basins and the optimal estimation of trend value by BMA-merging are shown in histograms and inserted values, respectively.

Table 1

List of non-groundwater components provided by global hydrological models, land surface models, and reanalysis data product.

Data	Spatial Resolution	SWS	SMS	SWE	CWS
WGHM	0.5° × 0.5°	√	√	√	√
PCR-GLOBWB 2	0.5° × 0.5°	√	√	√
CLSM	1° × 1°		√	√	√
NOAH	1° × 1°		√	√	√
VIC	1° × 1°		√	√	√
ERA5-land	0.1° × 0.1°		√	√

Table 2

Spatial mean statistics of multi-source estimated annual water storage trends as well as their SD and CV values in 10 river basins of mainland China.

Basin ID	TWSA			SWSA			GWEA			SWEA		CWSA a				SMSA
	Trend(mm/yr)	SD(mm/yr)	CV(%)	Trend(mm/yr)	SD(mm/yr)	CV(%)	Trend(mm/yr)	SD(mm/yr)	CV(%)	Trend(mm/yr)	SD(mm/yr)	CV(%)	Trend(10−3 mm/yr)	SD(10−3 mm/yr)	CV(%)	Trend(mm/yr)	SD(mm/yr)	CV(%)
SRB	0.91	1.72	188	0.51	0.49	96	—	—	—	0.45	0.24	53	0.77	0.89	115	1.33	1.92	145
LRB	−6.09	2.21	36	0.13	0.56	442	—	—	—	0.20	0.15	76	0.08	0.86	1100	0.06	2.37	4037
HRB	−13.36	4.46	33	0.38	0.77	201	—	—	—	0.01	0.06	996	−0.30	0.64	210	−2.19	3.29	150
YRB	−6.01	2.48	41	2.66	2.86	107	−0.01	0.01	106	0.02	0.10	560	−0.15	1.00	676	−1.24	2.50	201
HHB	−7.86	3.76	48	−3.37	2.65	79	—	—	—	−0.01	0.02	242	−1.00	1.16	113	−4.41	4.70	106
YZRB	2.78	2.36	85	0.71	1.05	149	−0.06	0.07	108	−0.01	0.11	786	0.86	1.43	165	−0.69	1.85	268
PRB	3.63	1.55	43	0.38	0.37	97	—	—	—	0.01	0.01	124	−0.30	1.00	331	1.44	1.56	108
SEB	4.31	1.95	45	2.08	1.67	80	—	—	—	−0.01	0.02	119	0.36	1.00	277	2.63	0.79	30
SWB	−9.60	5.18	54	0.03	0.24	818	−0.63	0.65	103	−0.65	1.80	275	−0.44	1.23	280	−2.23	2.74	123
CB	−2.31	2.91	126	0.55	0.72	130	−0.08	0.16	210	0.03	0.37	1161	0.09	0.78	830	−0.48	2.27	475

^a Statistics for annual CWSA trends and SD values are presented as scientific counts due to their small magnitude.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15112744/s1. Figure S1: Spatial distribution of annual trends for six GRACE products from 2003–2015 (in mm/yr): (a) CSR mascon; (b) GSFC mascon; (c) JPL mascon; (d) CSR SH; (e) GFZ SH; (f) JPL SH; Figure S2: Spatial distribution of annual trends for 2 SWSA products from 2003–2015 (in mm/yr): (a) WGHM; (b) PCR-GLOBWB; Figure S3: Spatial distribution of annual trends for two GWEA products from 2003–2015 (in mm/yr): (a) Degree-Day model (b) Global Glacier Dataset; Figure S4: Spatial distribution of annual trends for six SWEA products from 2003–2015 (in mm/yr): (a) NOAH (b) VIC; (c) CLSM; (d) ERA5-land; (e) WGHM; (f) PCR-GLOBWB; Figure S5: Spatial distribution of annual trends for six CWSA products from 2003–2015 (in mm/yr): (a) NOAH (b) VIC; (c) CLSM; (d) WGHM; Figure S6: Spatial distribution of annual trends for six SMSA products from 2003–2015 (in mm/yr): (a) NOAH (b) VIC; (c) CLSM; (d) ERA5-land; (e) WGHM; (f) PCR-GLOBWB; Figure S7: Spatial distribution of annual trends means of 3456 GWSA results from 2003–2015 (in mm/yr). Table S1: Inventory of degree-day model-driven data and global glacier datasets.

References

1. Aeschbach-Hertig, W.; Gleeson, T. Regional Strategies for the Accelerating Global Problem of Groundwater Depletion. Nat. Geosci.; 2012; 5, pp. 853-861. [DOI: https://dx.doi.org/10.1038/ngeo1617]

2. Dangar, S.; Asoka, A.; Mishra, V. Causes and Implications of Groundwater Depletion in India: A Review. J. Hydrol.; 2021; 596, 126103. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126103]

3. Diffenbaugh, N.S.; Swain, D.L.; Touma, D. Anthropogenic Warming Has Increased Drought Risk in California. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 3931-3936. [DOI: https://dx.doi.org/10.1073/pnas.1422385112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25733875]

4. Le Brocque, A.F.; Kath, J.; Reardon-Smith, K. Chronic Groundwater Decline: A Multi-Decadal Analysis of Groundwater Trends under Extreme Climate Cycles. J. Hydrol.; 2018; 561, pp. 976-986. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2018.04.059]

5. Scanlon, B.R.; Duncan, I.; Reedy, R.C. Drought and the Water–Energy Nexus in Texas. Environ. Res. Lett.; 2013; 8, 045033. [DOI: https://dx.doi.org/10.1088/1748-9326/8/4/045033]

6. Rodell, M.; Velicogna, I.; Famiglietti, J.S. Satellite-Based Estimates of Groundwater Depletion in India. Nature; 2009; 460, pp. 999-1002. [DOI: https://dx.doi.org/10.1038/nature08238] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19675570]

7. Rodell, M.; Famiglietti, J.S. The Potential for Satellite-Based Monitoring of Groundwater Storage Changes Using GRACE: The High Plains Aquifer, Central US. J. Hydrol.; 2002; 263, pp. 245-256. [DOI: https://dx.doi.org/10.1016/S0022-1694(02)00060-4]

8. Wahr, J.; Molenaar, M.; Bryan, F. Time Variability of the Earth’s Gravity Field: Hydrological and Oceanic Effects and Their Possible Detection Using GRACE. J. Geophys. Res.; 1998; 103, pp. 30205-30229. [DOI: https://dx.doi.org/10.1029/98JB02844]

9. Vissa, N.K.; Anandh, P.C.; Behera, M.M.; Mishra, S. ENSO-Induced Groundwater Changes in India Derived from GRACE and GLDAS. J. Earth Syst. Sci.; 2019; 128, 115. [DOI: https://dx.doi.org/10.1007/s12040-019-1148-z]

10. Scanlon, B.R.; Longuevergne, L.; Long, D. Ground Referencing GRACE Satellite Estimates of Groundwater Storage Changes in the California Central Valley, USA. Water Resour. Res.; 2012; 48, W04520. [DOI: https://dx.doi.org/10.1029/2011WR011312]

11. Gong, H.; Pan, Y.; Zheng, L.; Li, X.; Zhu, L.; Zhang, C.; Huang, Z.; Li, Z.; Wang, H.; Zhou, C. Long-Term Groundwater Storage Changes and Land Subsidence Development in the North China Plain (1971–2015). Hydrogeol. J.; 2018; 26, pp. 1417-1427. [DOI: https://dx.doi.org/10.1007/s10040-018-1768-4]

12. Kiss, A.; Földváry, L. Uncertainty of GRACE-Borne Long Periodic and Secular Ice Mass Variations in Antarctica. Acta Geod. Geophys.; 2017; 52, pp. 497-510. [DOI: https://dx.doi.org/10.1007/s40328-016-0185-1]

13. Chen, J.L.; Wilson, C.R.; Tapley, B.D.; Ries, J.C. Low Degree Gravitational Changes from GRACE: Validation and Interpretation. Geophys. Res. Lett.; 2004; 31, L22607. [DOI: https://dx.doi.org/10.1029/2004GL021670]

14. Swenson, S.; Wahr, J. Post-Processing Removal of Correlated Errors in GRACE Data. Geophys. Res. Lett.; 2006; 33, L08402. [DOI: https://dx.doi.org/10.1029/2005GL025285]

15. Swenson, S.; Wahr, J. Methods for Inferring Regional Surface-Mass Anomalies from Gravity Recovery and Climate Experiment (GRACE) Measurements of Time-Variable Gravity. J. Geophys. Res.; 2002; 107, pp. ETG 3-1-ETG 3-13. [DOI: https://dx.doi.org/10.1029/2001JB000576]

16. Colombo, O.L. Numerical Methods for Harmonic Analysis on the Sphere; Ohio State Univerisy of Columbus Deparment of Geodetic Science and Surveying: Columbus, OH, USA, 1981.

17. Sideris, M.G.; Tziavos, I.N. FFT-Evaluation and Applications of Gravity-Field Convolution Integrals with Mean and Point Data. Bull. Geod.; 1988; 62, pp. 521-540. [DOI: https://dx.doi.org/10.1007/BF02520242]

18. Eicker, A.; Mayer-Guerr, T.; Kurtenbach, E. Challenges in Deriving Trends from GRACE. Geodesy for Planet Earth; Kenyon, S.; Pacino, M.C.; Marti, U. International Association of Geodesy Symposia Springer: Berlin/Heidelberg, Germany, 2012; Volume 136, pp. 153-160. ISBN 978-3-642-20337-4

19. Baur, O. On the Computation of Mass-Change Trends from GRACE Gravity Field Time-Series. J. Geodyn.; 2012; 61, pp. 120-128. [DOI: https://dx.doi.org/10.1016/j.jog.2012.03.007]

20. Földváry, L. Desmoothing of Averaged Periodical Signals for Geodetic Applications. Geophys. J. Int.; 2015; 201, pp. 1235-1250. [DOI: https://dx.doi.org/10.1093/gji/ggv092]

21. Scanlon, B.R.; Zhang, Z.; Save, H.; Wiese, D.N.; Landerer, F.W.; Long, D.; Longuevergne, L.; Chen, J. Global Evaluation of New GRACE Mascon Products for Hydrologic Applications. Water Resour. Res.; 2016; 52, pp. 9412-9429. [DOI: https://dx.doi.org/10.1002/2016WR019494]

22. Tangdamrongsub, N.; Hwang, C.; Shum, C.K.; Wang, L. Regional Surface Mass Anomalies from GRACE KBR Measurements: Application of L-Curve Regularization and a Priori Hydrological Knowledge. J. Geophys. Res.; 2012; 117, B11406. [DOI: https://dx.doi.org/10.1029/2012JB009310]

23. Lemoine, F.G.; Luthcke, S.B.; Rowlands, D.D.; Chinn, D.S.; Klosko, S.M.; Cox, C.M. The Use of Mascons to Resolve Time-Variable Gravity from GRACE. Dynamic Planet; Tregoning, P.; Rizos, C. International Association of Geodesy Symposia Springer: Berlin/Heidelberg, Germany, 2007; Volume 130, pp. 231-236. ISBN 978-3-540-49349-5

24. Save, H.; Bettadpur, S.; Tapley, B.D. High-Resolution CSR GRACE RL05 Mascons: High-Resolution CSR GRACE RL05 Mascons. J. Geophys. Res. Solid Earth; 2016; 121, pp. 7547-7569. [DOI: https://dx.doi.org/10.1002/2016JB013007]

25. Chen, J.; Famigliett, J.S.; Scanlon, B.R.; Rodell, M. Groundwater Storage Changes: Present Status from GRACE Observations. Surv. Geophys.; 2016; 37, pp. 397-417. [DOI: https://dx.doi.org/10.1007/s10712-015-9332-4]

26. Joodaki, G.; Wahr, J.; Swenson, S. Estimating the Human Contribution to Groundwater Depletion in the Middle East, from GRACE Data, Land Surface Models, and Well Observations. Water Resour. Res.; 2014; 50, pp. 2679-2692. [DOI: https://dx.doi.org/10.1002/2013WR014633]

27. Bonsor, H.; Shamsudduha, M.; Marchant, B.; MacDonald, A.; Taylor, R. Seasonal and Decadal Groundwater Changes in African Sedimentary Aquifers Estimated Using GRACE Products and LSMs. Remote Sens.; 2018; 10, 904. [DOI: https://dx.doi.org/10.3390/rs10060904]

28. Jin, S.; Feng, G. Large-Scale Variations of Global Groundwater from Satellite Gravimetry and Hydrological Models, 2002–2012. Glob. Planet. Change; 2013; 106, pp. 20-30. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2013.02.008]

29. Yin, W.; Li, T.; Zheng, W.; Hu, L.; Han, S.-C.; Tangdamrongsub, N.; Šprlák, M.; Huang, Z. Improving Regional Groundwater Storage Estimates from GRACE and Global Hydrological Models over Tasmania, Australia. Hydrogeol. J.; 2020; 28, pp. 1809-1825. [DOI: https://dx.doi.org/10.1007/s10040-020-02157-3]

30. Elsaka, B.; Abdelmohsen, K.; Alshehri, F.; Zaki, A.; El-Ashquer, M. Mass Variations in Terrestrial Water Storage over the Nile River Basin and Mega Aquifer System as Deduced from GRACE-FO Level-2 Products and Precipitation Patterns from GPCP Data. Water; 2022; 14, 3920. [DOI: https://dx.doi.org/10.3390/w14233920]

31. Shao, C.; Liu, Y. Analysis of Groundwater Storage Changes and Influencing Factors in China Based on GRACE Data. Atmosphere; 2023; 14, 250. [DOI: https://dx.doi.org/10.3390/atmos14020250]

32. Liu, F.; Kang, P.; Zhu, H.; Han, J.; Huang, Y. Analysis of Spatiotemporal Groundwater-Storage Variations in China from GRACE. Water; 2021; 13, 2378. [DOI: https://dx.doi.org/10.3390/w13172378]

33. Yin, Z.; Xu, Y.; Zhu, X.; Zhao, J.; Yang, Y.; Li, J. Variations of Groundwater Storage in Different Basins of China over Recent Decades. J. Hydrol.; 2021; 598, 126282. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126282]

34. Zhang, J.; Liu, K.; Wang, M. Seasonal and Interannual Variations in China’s Groundwater Based on GRACE Data and Multisource Hydrological Models. Remote Sens.; 2020; 12, 845. [DOI: https://dx.doi.org/10.3390/rs12050845]

35. Feng, W.; Shum, C.; Zhong, M.; Pan, Y. Groundwater Storage Changes in China from Satellite Gravity: An Overview. Remote Sens.; 2018; 10, 674. [DOI: https://dx.doi.org/10.3390/rs10050674]

36. Zhang, C.; Duan, Q.; Yeh, P.J.-F.; Pan, Y.; Gong, H.; Moradkhani, H.; Gong, W.; Lei, X.; Liao, W.; Xu, L. et al. Sub-Regional Groundwater Storage Recovery in North China Plain after the South-to-North Water Diversion Project. J. Hydrol.; 2021; 597, 126156. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126156]

37. Huang, Z.; Yeh, P.J.-F.; Jiao, J.J.; Luo, X.; Pan, Y.; Long, Y.; Zhang, C.; Zheng, L. A New Approach for Assessing Groundwater Recharge by Combining GRACE and Baseflow With Case Studies in Karst Areas of Southwest China. Water Resour. Res.; 2023; 59, e2022WR032091. [DOI: https://dx.doi.org/10.1029/2022WR032091]

38. Zhang, X.; Ren, L.; Feng, W. Comparison of the Shallow Groundwater Storage Change Estimated by a Distributed Hydrological Model and GRACE Satellite Gravimetry in a Well-Irrigated Plain of the Haihe River Basin, China. J. Hydrol.; 2022; 610, 127799. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.127799]

39. Yi, S.; Wang, Q.; Sun, W. Basin Mass Dynamic Changes in China from GRACE Based on a Multibasin Inversion Method: Basin Mass in China. J. Geophys. Res. Solid Earth; 2016; 121, pp. 3782-3803. [DOI: https://dx.doi.org/10.1002/2015JB012608]

40. Alley, W.M.; Konikow, L.F. Bringing GRACE Down to Earth. Groundwater; 2015; 53, pp. 826-829. [DOI: https://dx.doi.org/10.1111/gwat.12379]

41. Herman, T.; Tregear, T.R. A Geography of China. Geogr. J.; 1966; 132, 101. [DOI: https://dx.doi.org/10.2307/1793089]

42. Gao, G.; Chen, D.; Xu, C.; Simelton, E. Trend of Estimated Actual Evapotranspiration over China during 1960–2002. J. Geophys. Res.; 2007; 112, D11120. [DOI: https://dx.doi.org/10.1029/2006JD008010]

43. Liu, B.; Speed, R. Water Resources Management in the People’s Republic of China. Int. J. Water Resour. Dev.; 2009; 25, pp. 193-208. [DOI: https://dx.doi.org/10.1080/07900620902868596]

44. Jia, X.; O’Connor, D.; Hou, D.; Jin, Y.; Li, G.; Zheng, C.; Ok, Y.S.; Tsang, D.C.W.; Luo, J. Groundwater Depletion and Contamination: Spatial Distribution of Groundwater Resources Sustainability in China. Sci. Total Environ.; 2019; 672, pp. 551-562. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.03.457] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30965267]

45. Zheng, L.; Pan, Y.; Gong, H.; Huang, Z.; Zhang, C. Comparing Groundwater Storage Changes in Two Main Grain Producing Areas in China: Implications for Sustainable Agricultural Water Resources Management. Remote Sens.; 2020; 12, 2151. [DOI: https://dx.doi.org/10.3390/rs12132151]

46. Tapley, B.D.; Bettadpur, S.; Watkins, M.; Reigber, C. The Gravity Recovery and Climate Experiment: Mission Overview and Early Results. Geophys. Res. Lett.; 2004; 31, L09607. [DOI: https://dx.doi.org/10.1029/2004GL019920]

47. Watkins, M.M.; Wiese, D.N.; Yuan, D.-N.; Boening, C.; Landerer, F.W. Improved Methods for Observing Earth’s Time Variable Mass Distribution with GRACE Using Spherical Cap Mascons: Improved Gravity Observations from GRACE. J. Geophys. Res. Solid Earth; 2015; 120, pp. 2648-2671. [DOI: https://dx.doi.org/10.1002/2014JB011547]

48. Loomis, B.D.; Luthcke, S.B.; Sabaka, T.J. Correction to: Regularization and Error Characterization of GRACE Mascons. J Geod; 2019; 93, 1399. [DOI: https://dx.doi.org/10.1007/s00190-019-01286-2]

49. Sutanudjaja, E.H.; van Beek, R.; Wanders, N.; Wada, Y.; Bosmans, J.H.C.; Drost, N.; van der Ent, R.J.; de Graaf, I.E.M.; Hoch, J.M.; de Jong, K. et al. PCR-GLOBWB 2: A 5 Arcmin Global Hydrological and Water Resources Model. Geosci. Model Dev.; 2018; 11, pp. 2429-2453. [DOI: https://dx.doi.org/10.5194/gmd-11-2429-2018]

50. Müller Schmied, H.; Cáceres, D.; Eisner, S.; Flörke, M.; Herbert, C.; Niemann, C.; Peiris, T.A.; Popat, E.; Portmann, F.T.; Reinecke, R. et al. The Global Water Resources and Use Model WaterGAP v2.2d: Model Description and Evaluation. Geosci. Model Dev.; 2021; 14, pp. 1037-1079. [DOI: https://dx.doi.org/10.5194/gmd-14-1037-2021]

51. Li, B.; Beaudoing, H.; Rodell, M. NASA/GSFC/HSL, GLDAS Catchment Land Surface Model L4 monthly 1.0 × 1.0 Degree V2.1; Goddard Earth Sciences Data and Information Services Center (GES DISC): Greenbelt, MD, USA, 2020; [DOI: https://dx.doi.org/10.5067/FOUXNLXFAZNY]

52. Beaudoing, H.; Rodell, M. NASA/GSFC/HSL, GLDAS Noah Land Surface Model L4 Monthly 1.0 × 1.0 Degree V2.0; Goddard Earth Sciences Data and Information Services Center (GES DISC): Greenbelt, MD, USA, 2019; [DOI: https://dx.doi.org/10.5067/QN80TO7ZHFJZ]

53. Beaudoing, H.; Rodell, M. NASA/GSFC/HSL, GLDAS VIC Land Surface Model L4 Monthly 1.0 × 1.0 Degree V2.1; Goddard Earth Sciences Data and Information Services Center (GES DISC): Greenbelt, MD, USA, 2020; [DOI: https://dx.doi.org/10.5067/VWTH7S6218SG]

54. Muñoz Sabater, J. ERA5-Land Hourly Data from 1950 to Present; Copernicus Climate Change Service (C3S) Climate Data Store (CDS): Brussels, Belgium, 2019; [DOI: https://dx.doi.org/10.24381/cds.e2161bac]

55. Döll, P.; Kaspar, F.; Lehner, B. A Global Hydrological Model for Deriving Water Availability Indicators: Model Tuning and Validation. J. Hydrol.; 2003; 270, pp. 105-134. [DOI: https://dx.doi.org/10.1016/S0022-1694(02)00283-4]

56. Müller Schmied, H.; Eisner, S.; Franz, D.; Wattenbach, M.; Portmann, F.T.; Flörke, M.; Döll, P. Sensitivity of Simulated Global-Scale Freshwater Fluxes and Storages to Input Data, Hydrological Model Structure, Human Water Use and Calibration. Hydrol. Earth Syst. Sci.; 2014; 18, pp. 3511-3538. [DOI: https://dx.doi.org/10.5194/hess-18-3511-2014]

57. Rodell, M.; Houser, P.R.; Jambor, U.; Gottschalck, J.; Mitchell, K.; Meng, C.-J.; Arsenault, K.; Cosgrove, B.; Radakovich, J.; Bosilovich, M. et al. The Global Land Data Assimilation System. Bull. Amer. Meteor. Soc.; 2004; 85, pp. 381-394. [DOI: https://dx.doi.org/10.1175/BAMS-85-3-381]

58. Kumar, S.V.; Jasinski, M.; Mocko, D.M.; Rodell, M.; Borak, J.; Li, B.; Beaudoing, H.K.; Peters-Lidard, C.D. NCA-LDAS Land Analysis: Development and Performance of a Multisensor, Multivariate Land Data Assimilation System for the National Climate Assessment. J. Hydrometeorol.; 2019; 20, pp. 1571-1593. [DOI: https://dx.doi.org/10.1175/JHM-D-17-0125.1]

59. Kumar, S.V.; Zaitchik, B.F.; Peters-Lidard, C.D.; Rodell, M.; Reichle, R.; Li, B.; Jasinski, M.; Mocko, D.; Getirana, A.; De Lannoy, G. et al. Assimilation of Gridded GRACE Terrestrial Water Storage Estimates in the North American Land Data Assimilation System. J. Hydrometeorol.; 2016; 17, pp. 1951-1972. [DOI: https://dx.doi.org/10.1175/JHM-D-15-0157.1]

60. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D. et al. The ERA5 Global Reanalysis. Q. J. R. Meteorol. Soc.; 2020; 146, pp. 1999-2049. [DOI: https://dx.doi.org/10.1002/qj.3803]

61. Muñoz-Sabater, J.; Dutra, E.; Agustí-Panareda, A.; Albergel, C.; Arduini, G.; Balsamo, G.; Boussetta, S.; Choulga, M.; Harrigan, S.; Hersbach, H. et al. ERA5-Land: A State-of-the-Art Global Reanalysis Dataset for Land Applications. Earth Syst. Sci. Data; 2021; 13, pp. 4349-4383. [DOI: https://dx.doi.org/10.5194/essd-13-4349-2021]

62. Yao, T.; Thompson, L.G.; Mosbrugger, V.; Zhang, F.; Ma, Y.; Luo, T.; Xu, B.; Yang, X.; Joswiak, D.R.; Wang, W. et al. Third Pole Environment (TPE). Environ. Dev.; 2012; 3, pp. 52-64. [DOI: https://dx.doi.org/10.1016/j.envdev.2012.04.002]

63. Wu, Y. Estimating Glacial Material Balance in the Asian High Sierra 1999–2019 Using Degree-Day Models and Gravity Satellites. Master’s Thesis; Capital Normal University: Beijing, China, 2021.

64. Hugonnet, R.; McNabb, R.; Berthier, E.; Menounos, B.; Nuth, C.; Girod, L.; Farinotti, D.; Huss, M.; Dussaillant, I.; Brun, F. et al. Accelerated Global Glacier Mass Loss in the Early Twenty-First Century. Nature; 2021; 592, pp. 726-731. [DOI: https://dx.doi.org/10.1038/s41586-021-03436-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33911269]

65. Jóhannesson, T.; Sigurdsson, O.; Laumann, T.; Kennett, M. Degree-Day Glacier Mass-Balance Modelling with Applications to Glaciers in Iceland, Norway and Greenland. J. Glaciol.; 1995; 41, pp. 345-358. [DOI: https://dx.doi.org/10.3189/S0022143000016221]

66. Anderson, B.; Lawson, W.; Owens, I.; Goodsell, B. Past and Future Mass Balance of ‘Ka Roimata o Hine Hukatere’ Franz Josef Glacier, New Zealand. J. Glaciol.; 2006; 52, pp. 597-607. [DOI: https://dx.doi.org/10.3189/172756506781828449]

67. Möller, M.; Schneider, C. Calibration of Glacier Volume–Area Relations from Surface Extent Fluctuations and Application to Future Glacier Change. J. Glaciol.; 2010; 56, pp. 33-40. [DOI: https://dx.doi.org/10.3189/002214310791190866]

68. Yeh, P.J.-F.; Swenson, S.C.; Famiglietti, J.S.; Rodell, M. Remote Sensing of Groundwater Storage Changes in Illinois Using the Gravity Recovery and Climate Experiment (GRACE). Water Resour. Res.; 2006; 42, W12203. [DOI: https://dx.doi.org/10.1029/2006WR005374]

69. Vrugt, J.A.; Robinson, B.A. Treatment of Uncertainty Using Ensemble Methods: Comparison of Sequential Data Assimilation and Bayesian Model Averaging. Water Resour. Res.; 2007; 43, W01411. [DOI: https://dx.doi.org/10.1029/2005WR004838]

70. Liu, J.S. Monte Carlo Strategies in Scientific Computing; Springer Series in Statistics; Springer: New York, NY, USA, 2004; ISBN 978-0-387-76369-9

71. Schmidt, R.; Flechtner, F.; Meyer, U.; Neumayer, K.-H.; Dahle, C.; König, R.; Kusche, J. Hydrological Signals Observed by the GRACE Satellites. Surv. Geophys.; 2008; 29, pp. 319-334. [DOI: https://dx.doi.org/10.1007/s10712-008-9033-3]

72. Chen, Q.; Shen, Y.; Kusche, J.; Chen, W.; Chen, T.; Zhang, X. High-Resolution GRACE Monthly Spherical Harmonic Solutions. J. Geophys. Res. Solid Earth; 2021; 126, e2019JB018892. [DOI: https://dx.doi.org/10.1029/2019JB018892]

73. Hu, L.; Jiao, J.J. Calibration of a Large-Scale Groundwater Flow Model Using GRACE Data: A Case Study in the Qaidam Basin, China. Hydrogeol. J.; 2015; 23, pp. 1305-1317. [DOI: https://dx.doi.org/10.1007/s10040-015-1278-6]

74. Hu, Z.; Zhou, Q.; Chen, X.; Chen, D.; Li, J.; Guo, M.; Yin, G.; Duan, Z. Groundwater Depletion Estimated from GRACE: A Challenge of Sustainable Development in an Arid Region of Central Asia. Remote Sens.; 2019; 11, 1908. [DOI: https://dx.doi.org/10.3390/rs11161908]

75. Long, D.; Yang, Y.; Wada, Y.; Hong, Y.; Liang, W.; Chen, Y.; Yong, B.; Hou, A.; Wei, J.; Chen, L. Deriving Scaling Factors Using a Global Hydrological Model to Restore GRACE Total Water Storage Changes for China’s Yangtze River Basin. Remote Sens. Environ.; 2015; 168, pp. 177-193. [DOI: https://dx.doi.org/10.1016/j.rse.2015.07.003]

76. Wang, X.; De Linage, C.; Famiglietti, J.; Zender, C.S. Gravity Recovery and Climate Experiment (GRACE) Detection of Water Storage Changes in the Three Gorges Reservoir of China and Comparison with in Situ Measurements. Water Resour. Res.; 2011; 47, W12502. [DOI: https://dx.doi.org/10.1029/2011WR010534]

77. Scanlon, B.R.; Zhang, Z.; Save, H.; Sun, A.Y.; Müller Schmied, H.; Van Beek, L.P.H.; Wiese, D.N.; Wada, Y.; Long, D.; Reedy, R.C. et al. Global Models Underestimate Large Decadal Declining and Rising Water Storage Trends Relative to GRACE Satellite Data. Proc. Natl. Acad. Sci. USA; 2018; 115, pp. E1080-E1089. [DOI: https://dx.doi.org/10.1073/pnas.1704665115]

78. Deng, C.; Zhang, W. Spatial Distribution Pattern of Degree–Day Factors of Glaciers on the Qinghai–Tibetan Plateau. Env. Monit. Assess.; 2018; 190, 475. [DOI: https://dx.doi.org/10.1007/s10661-018-6860-7]

79. Hughes, P.D.; Braithwaite, R.J. Application of a Degree-Day Model to Reconstruct Pleistocene Glacial Climates. Quat. Res.; 2008; 69, pp. 110-116. [DOI: https://dx.doi.org/10.1016/j.yqres.2007.10.008]

80. Wu, L.; Li, H.; Wang, L. Application of a Degree-Day Model for Determination of Mass Balance of Urumqi Glacier No. 1, Eastern Tianshan, China. J. Earth Sci.; 2011; 22, pp. 470-481. [DOI: https://dx.doi.org/10.1007/s12583-011-0201-x]

81. Zhang, S.; Ye, B.; Liu, S.; Zhang, X.; Hagemann, S. A Modified Monthly Degree-Day Model for Evaluating Glacier Runoff Changes in China. Part I: Model Development. Hydrol. Process.; 2012; 26, pp. 1686-1696. [DOI: https://dx.doi.org/10.1002/hyp.8286]

82. Wu, Z.; Feng, H.; He, H.; Zhou, J.; Zhang, Y. Evaluation of Soil Moisture Climatology and Anomaly Components Derived From ERA5-Land and GLDAS-2.1 in China. Water Resour Manag.; 2021; 35, pp. 629-643. [DOI: https://dx.doi.org/10.1007/s11269-020-02743-w]

83. Wang, S.; Cui, G.; Li, X.; Liu, Y.; Li, X.; Tong, S.; Zhang, M. GRACE Satellite-Based Analysis of Spatiotemporal Evolution and Driving Factors of Groundwater Storage in the Black Soil Region of Northeast China. Remote Sens.; 2023; 15, 704. [DOI: https://dx.doi.org/10.3390/rs15030704]

84. Cheng, M.; Zhong, L.; Ma, Y.; Zou, M.; Ge, N.; Wang, X.; Hu, Y. A Study on the Assessment of Multi-Source Satellite Soil Moisture Products and Reanalysis Data for the Tibetan Plateau. Remote Sens.; 2019; 11, 1196. [DOI: https://dx.doi.org/10.3390/rs11101196]

85. Jiao, J.J.; Zhang, X.; Wang, X. Satellite-Based Estimates of Groundwater Depletion in the Badain Jaran Desert, China. Sci Rep; 2015; 5, 8960. [DOI: https://dx.doi.org/10.1038/srep08960]

86. Cao, G.; Han, D.; Song, X. Evaluating Actual Evapotranspiration and Impacts of Groundwater Storage Change in the North China Plain. Hydrol. Process.; 2014; 28, pp. 1797-1808. [DOI: https://dx.doi.org/10.1002/hyp.9732]

87. Xie, X.; Xu, C.; Wen, Y.; Li, W. Monitoring Groundwater Storage Changes in the Loess Plateau Using GRACE Satellite Gravity Data, Hydrological Models and Coal Mining Data. Remote Sens.; 2018; 10, 605. [DOI: https://dx.doi.org/10.3390/rs10040605]

88. Zhang, K.; Xie, X.; Zhu, B.; Meng, S.; Yao, Y. Unexpected Groundwater Recovery with Decreasing Agricultural Irrigation in the Yellow River Basin. Agric. Water Manag.; 2019; 213, pp. 858-867. [DOI: https://dx.doi.org/10.1016/j.agwat.2018.12.009]

89. Zhong, Y.; Zhong, M.; Feng, W.; Zhang, Z.; Shen, Y.; Wu, D. Groundwater Depletion in the West Liaohe River Basin, China and Its Implications Revealed by GRACE and In Situ Measurements. Remote Sens.; 2018; 10, 493. [DOI: https://dx.doi.org/10.3390/rs10040493]

90. Su, Y.; Guo, B.; Zhou, Z.; Zhong, Y.; Min, L. Spatio-Temporal Variations in Groundwater Revealed by GRACE and Its Driving Factors in the Huang-Huai-Hai Plain, China. Sensors; 2020; 20, 922. [DOI: https://dx.doi.org/10.3390/s20030922] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32050517]

91. Zou, Y.; Kuang, X.; Feng, Y.; Jiao, J.J.; Liu, J.; Wang, C.; Fan, L.; Wang, Q.; Chen, J.; Ji, F. et al. Solid Water Melt Dominates the Increase of Total Groundwater Storage in the Tibetan Plateau. Geophys. Res. Lett.; 2022; 49, e2022GL100092. [DOI: https://dx.doi.org/10.1029/2022GL100092]