Improving the Accuracy of Groundwater Storage Estimates Based on Groundwater Weighted Fusion Model

Abstract

It is an effective measure to estimate groundwater storage anomalies (GWSA) by combining Gravity Recovery and Climate Experiment (GRACE) data and hydrological models. However, GWSA results based on a single hydrological model and GRACE data may have greater uncertainties, and it is difficult to verify in some regions where in situ groundwater-level measurements are limited. First, to solve this problem, a groundwater weighted fusion model (GWFM) is presented, based on the extended triple collocation (ETC) method. Second, the Shiyang River Basin (SYRB) is taken as an example, and in situ groundwater-level measurements are used to evaluate the performance of the GWFM. The comparison indicates that the correlation coefficient (CC) and Nash-Sutcliffe efficiency coefficient (NSE) are increased by 9–40% and 23–657%, respectively, relative to the original results. Moreover, the root mean squared error (RMSE) is reduced by 9–28%, which verifies the superiority of the GWFM. Third, the spatiotemporal distribution and influencing factors of GWSA in the Hexi Corridor (HC) are comprehensively analyzed during the period between 2003 and 2016. The results show that GWSA decline, with a trend of −2.37 ± 0.38 mm/yr from 2003 to 2010, and the downward trend after 2011 (−0.46 ± 1.35 mm/yr) slow down significantly compared to 2003–2010. The spatial distribution obtained by the GWFM is more reliable compared to the arithmetic average results, and GWFM-based GWSA fully retain the advantages of different models, especially in the southeastern part of the SYRB. Additionally, a simple index is used to evaluate the contributions of climatic factors and human factors to groundwater storage (GWS) in the HC and its different subregions. The index indicates that climate factors occupy a dominant position in the SLRB and SYRB, while human factors have a significant impact on GWS in the Heihe River Basin (HRB). This study can provide suggestions for the management and assessments of groundwater resources in some arid regions.

1. Introduction

As an important component of terrestrial water storage (TWS), groundwater plays a key role in domestic, agriculture and industrial use, as well as ecosystems [1,2]. More than 38% of the world’s population lives in arid or semi-arid zones [3], where groundwater is usually the dominant freshwater resource, supplying domestic use and irrigation water [4]. Especially in northwest China, groundwater resources have been facing the risk of depletion, which may lead to the ecological environment of the region losing its ability to self-repair and endangering local ecological security [5]. Therefore, accurate estimation of groundwater storage anomalies (GWSA) is essential for the effective use of local groundwater resources. The traditional method of monitoring groundwater level mainly uses monitoring wells. However, monitoring wells are scarce, and the observation records are short and discontinuous, restricting research related to GWSA [2,6,7]. Therefore, it is important to seek an alternative method to obtain ground-based network data for monitoring of large-scale groundwater storage (GWS) variations.

Since March 2002, Gravity Recovery and Climate Experiment (GRACE) satellites have provided an opportunity to assess global TWS changes, with a resolution of ~300 km [8,9,10]. Currently, the GRACE gravity satellites are the only way to sense water storage at all levels, including soil moisture (SM), snow-water equivalent (SWE), canopy water storage (CWS), and GWS [10]. To isolate the GWS component from TWS, water storage changes of other components have to be estimated based on hydrological models [11,12]. At present, there are several frequently employed hydrological models and reanalysis datasets, such as the Global Land Data Assimilation System (GLDAS) [13], the WaterGAP Global Hydrology model (WGHM) [14], and the ERA5 reanalysis dataset [15,16]. Furthermore, previous studies have demonstrated the effectiveness of GRACE observations to estimate GWSA in many typically regions of the world, e.g., the Central Valley of California [17,18], northwest India [9,19,20], and the North China Plain [1,21,22].

Currently, most studies mainly rely on a single hydrological model to separate GWS components from GRACE-derived TWS [23,24,25]. However, the accuracy of these models is restricted by uncertainties in climate forcing (particularly precipitation), model parameters, and deficiencies in model structure [26,27,28,29,30,31]. Therefore, the effective combination of multiple models can improve the performance of hydrological simulations relative to a single model. For instance, Shamseldin et al. [32] used the method of multi-model ensemble to develop more skillful and reliable probabilistic hydrologic prediction. The results confirmed that better estimates of water storage can be obtained by combining the model outputs of different hydrological models. Long et al. [33] used the Bayesian model-averaging technique, which can merge multiple TWS products to analyze the spatiotemporal variability of TWS. Mehrnegar [27] presented the dynamic model-data-averaging method, which can be used to merge multiple TWS simulations. The result indicated that linear trends and seasonality within global hydrological models can be improved by using the dynamic model-data-averaging method. These multi-model techniques prove to provide accurate estimates by combining different models according to the different weighting strategies [32,34].

Triple collocation (TC) is a statistical method to estimate the random-error variance of three independent datasets [35]. Currently, the TC method has been used to estimate measurement errors of GRACE data [36]. Specifically, Khaki et al. [37] and Nigatu et al. [38] estimated the changes in key water-storage components by using the GRACE data and soil-moisture data based on the TC analysis method. Yin and Park [39] proposed a simple least-square merging approach using error characteristics quantified from the TC approach to estimate weight. Compared to the classic TC approach, the extended triple collocation (ETC), proposed by McColl [40], can obtain an additional evaluation index, that is, the correlation coefficient relative to the unknown true value. Up to now, there are few studies that have merged datasets from different sources based on the ETC method.

The Hexi Corridor (HC) is one of the most agriculturally rich areas of northwest China, which is characterized as an irrigation district of “no irrigation, no agriculture” [5,41]. Moreover, groundwater resources have been depleted on a large scale in the area due to poor management of groundwater exploitation [42]. The policy of building a water-saving society was introduced in Zhangye City of Gansu Province in 2001 [43]. The government initiated a policy called the Key Governance Planning Project of the Shiyang River Basin in 2007, which aimed to improve the ecological conditions of the area [44,45]. Accurate estimation of GWS is essential for understanding the complex hydrological process and formulating sustainable management policies for groundwater resources in the region.

The purpose of this study is to improve the accuracy of groundwater storage estimates in some regions where in situ groundwater-level measurements are limited and to quantify the impact of climate change and human activities. Specifically, a weighted fusion model is proposed, based on the squared correlation coefficient and error variance calculated by the ETC [40,46] method. The ratio of these two indicators is used to develop the groundwater weighted fusion model (GWFM), which is helpful in merging GWSA based on the GRACE and multiple hydrological models, and compare it with the original results. In addition, a simple and effective method is used to evaluate the contribution of climate factors and human factors to GWS.

2. Materials and Methods

2.1. Study Area

The HC is located in Northwest China (92°12′–104°20′ E, 37°17′–42°48′ N) (Figure 1a), including the five prefecture-level cities of Wuwei, Jinchang, Zhangye, Jiuquan, and Jiayuguan (Figure 1c, http://srtm.csi.cgiar.org (accessed on 10 August 2021)) [47]. The three major rivers of the HC, from west to east, are the Shule River, the Heihe River and the Shiyang River, originating from the Qilian Mountains. The HC is mainly covered by bare land and gobi, and agricultural land is concentrated in the Shiyang River Basin (SYRB) and Heihe River Basin (HRB) (Figure 1d, https://www.resdc.cn (accessed on 10 August 2021)) [48]. It belongs to the arid continental climate, with an average annual temperature ranging from 8 to 10 °C [42], and the annual evaporation is 1500–3200 mm [49]. The rainfall shows a decreasing trend from east to west (Figure 1b, http://data.cma.cn (accessed on 10 August 2021)), and ~80% of rainfall occurs during July–August [50]. The Qilian Mountains in the south nurture extensive glaciers and snow, which is an important freshwater resource for downstream cities critical to easing agricultural drought and other human activities in the HC [51].

2.2. Materials

2.2.1. GRACE Data

The GRACE RL05 Mascon solutions are utilized to derive TWS anomalies in this study, which are provided by the Center for Space Research (CSR) [52]. Monthly TWS anomalies are provided from April 2002 to June 2017, with a spatial resolution of 0.5 × 0.5°. The regularization constraint on mascon solutions is derived from original GRACE information with no empirical filtering post-processing [52,53]. Therefore, the product can capture all the signals observed by GRACE within the measurement noise level and be used without further processing [54]. Missing data in the CSR Mascon are filled by linear interpolation of the nearby monthly mean values [24].

2.2.2. Soil Moisture Datasets

GLDAS was jointly developed by the National Aeronautics and Space Administration and the National Oceanic and Atmospheric Administration, which can obtain land-surface state and flux with high time resolution (https://disc.gsfc.nasa.gov/ (accessed on 1 July 2021)) [13]. In this study, the monthly SM product provided by the GLDAS Noah model with a spatial resolution of 1.0 × 1.0° is used to estimate SM over the HC. For consistency of data resolution, the related datasets are interpolated into a spatial resolution of 0.5 × 0.5°. More details on various soil-moisture data used are summarized in

Table 1. Summary of soil-moisture products from GLDAS, WGHM, and ERA5-Land.

Datasets	Spatial Resolution	Temporal Resolution	Soil Layer	Depth (cm)
GLDAS-Noah	1.0 × 1.0°	monthly	4	0–10, 10–40, 40–100, 100–200
WGHM	0.5 × 0.5°	monthly	-	100–200
ERA5-Land	0.1 × 0.1°	monthly	4	0–7, 7–28, 28–100, 100–289

.

WGHM [14] was developed by the Institute of Physical Geography at the University of Frankfurt and provides information on spatiotemporal water-storage changes for most hydrological processes. This model accounts for four of the most important terrestrial water-storage components: surface water, snow, soil water, and groundwater storage [55]. The WGHM data were retrieved from https://doi.pangaea.de/10.1594/PANGAEA.918447 (accessed on 1 July 2021). The SM product provided by WGHM is used in this study, which is monthly data from January 2003 to December 2016 at a spatial resolution of 0.5 × 0.5°.

ERA5-Land [56] is a reanalysis dataset produced by replaying the land component of the ERA5 climate reanalysis (https://cds.climate.copernicus.eu/ (accessed on 1 July 2021)). It is one of the most modern and finest reanalysis datasets produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) within the Copernicus Climate Change Service. In this study, the SM product of ERA5-Land is employed, which is the monthly datasets, with a spatial resolution of 0.1 × 0.1° from 2003 to 2016. To maintain the same spatial resolution, the related datasets are interpolated into the 0.5 × 0.5° spatial resolution.

2.2.3. Groundwater Level from Wells

Groundwater monitoring data are collected from the groundwater yearbooks compiled by the China Institute of Geological Environment Monitoring (CIGEM), which is published by the Ministry of Land and Resources of the People’s Republic of China. Due to the sparse number of stations and a lack of continuous data at individual stations, the measured groundwater-level data of five wells from 2007 to 2014 are selected in the SYRB to verify the performance of GWFM-based GWSA in this study (shown in Figure 1b). The groundwater level can be converted to groundwater storage by multiplying by specific yield values. However, specific yield values are unknown, and the groundwater level is only used to verify the performance of GWFM. Therefore, there is no need to covert the levels to groundwater storage in this study to avoid possible errors associated within unknown specific yield values.

2.2.4. Auxiliary Data

The precipitation dataset is collected from the China Meteorological Data Service Center, based on the precipitation data of high-density ground stations in China (2472 national meteorological observatories). It uses the thin-plate splines method [57] of ANUSPLIN software for spatial interpolation to generate monthly grid data from 1961 to the present, with a spatial resolution of 0.5 × 0.5°. Additionally, evapotranspiration and temperature data from the GLDAS and ERA5-Land during 2003–2016 are collected to evaluate the impact of climate factors on GWS.

To evaluate the impact of human factors, the annual groundwater withdrawal data of the HC from 2003 to 2016 are also collected. They are collected from the Water Resources Bulletin of Gansu Province, which is published by the Gansu Provincial Department of Water Resources, China.

2.3. Methods

The flow chart of the study is shown in Figure 2. First, the GWSA results are estimated based on GRACE and three hydrological models (selected as GLDAS, ERA5, and WGHM). Second, the error variance and correlation coefficient of three GWSA derived by ETC are utilized for weight estimation, and then the GWSA results from different sources are merged based on the least squares framework. Furthermore, the merged GWSA are verified by the original results and in situ groundwater-level measurements. Finally, the impact of different factors on GWS is analyzed by combining climatic factors and water-consumption data.

2.3.1. Construction of GWFM

Hydrological simulation typically needs to be calibrated by the “true” value of the target variable. However, it is difficult to obtain measured data in some areas because of the uneven distribution of monitoring points. ETC is a statistical method to estimate the correlation coefficient and the random-error variance of three independent datasets. The prerequisites for the TC approach are: (i) linearity between the true hydrological signal and the observations, (ii) signal and error stationarity, (iii) independency between the errors and the hydrological signal (error orthogonality), and (iv) independency between the errors of each dataset (zero-error cross-correlation) [35]. This study uses the error variance and squared correlation coefficient calculated by the ETC method to develop the GWFM. The error model is given by [35,40]:

$$S_i={\alpha }_i+{\beta }_i\Theta +{\epsilon }_i$$

(1)

where $S_i (i=1,2,3)$ represents GWSA based on GRACE and three hydrological models; $\Theta$ denotes the unknown true hydrological signal; $S_i$ represents collocated measurement systems linearly related to the true unknown value, $\Theta$; ${\alpha }_i$ and ${\beta }_i$ represent the least-squares intercepts and slope, respectively; and ${\epsilon }_i$ represents additive zero-mean random noise.

Covariance estimation is used to solve random-error variance in this study. The covariances between the different datasets are given by [35,40]:

$$\begin{array}{l}cov(S_i,S_j)=E(S_i{S}_j)-E(S_i)E(S_j)={\beta }_i{\beta }_j{\sigma }_{\Theta }^2\\ +{\beta }_{i}cov(\Theta ,{\epsilon }_j)+{\beta }_{j}cov(\Theta ,{\epsilon }_i)+cov({\epsilon }_i,{\epsilon }_j)\end{array}$$

(2)

According to these four prerequisites of the TC approach, $cov({\epsilon }_i,{\epsilon }_j)=0 (i\ne j)$, $cov({\epsilon }_i,\Theta )=0$; the equation reduces to [35,40]:

$$C_{ij}=cov(S_i,S_j)=\left\{\begin{array}{c}{\beta }_i{\beta }_j{\sigma }_{\Theta }^2 (i\ne j)\\ {\beta }_i^2{\sigma }_{\Theta }^2+{\sigma }_{{\epsilon }_i}^2 (i=j)\end{array}\right.$$

(3)

where ${\sigma }_{{\epsilon }_i}^2$ represents the variance of random-error variance, ${\epsilon }_i$; ${\beta }_i^2{\sigma }_{\Theta }^2$ denotes the sensitivity of datasets, $S_i$, to changes in true signal. In other words, the higher ${\beta }_i$, the stronger the response of datasets,$S_i$, to hydrological signal. The sensitivity of each dataset can be calculated by combining their covariances [35,40]:

$${\beta }_i^2{\sigma }_{\Theta }^2=\left[\begin{array}{c}\frac{C_{12}{C}_{13}}{C_{23}}\\ \frac{C_{21}{C}_{23}}{C_{13}}\\ \frac{C_{31}{C}_{32}}{C_{12}}\end{array}\right]$$

(4)

The error variance can be obtained by subtracting the sensitivity ${\beta }_i^2{\sigma }_{\Theta }^2$ of each dataset from their total variance [35,40]:

$${\sigma }_{{\epsilon }_i}^2=\left[\begin{array}{c}{C}_{11}-\frac{C_{12}{C}_{13}}{C_{23}}\\ C_{22}-\frac{C_{21}{C}_{23}}{C_{13}}\\ C_{33}-\frac{C_{31}{C}_{32}}{C_{12}}\end{array}\right]$$

(5)

$${\rho }_{i,\Theta }^2=\frac{{\beta }_i^2{\sigma }_{\Theta }^2}{{\beta }_i^2{\sigma }_{\Theta }^2+{\sigma }_{{\epsilon }_i}^2}=\frac{SN{R}_i}{SN{R}_i+1}$$

(6)

$$SN{R}_i=\frac{var(S_i)}{var({\epsilon }_i)}=\frac{{\beta }_i^2{\sigma }_{\Theta }^2}{{\sigma }_{{\epsilon }_i}^2}$$

(7)

where ${\rho }_{i,\Theta }^2$ represents the squared correlation coefficient and $SNR$ represents the un-biased signal-to-noise ratio.

In this study, the error variance and squared correlation coefficient calculated by the above method are used to develop the GWFM. The detailed formula is as follows:

$$W_i=\left\{ \begin{array}{c}\frac{{\rho }_{i,\Theta }^2/{\sigma }_{{\epsilon }_i}^2}{{\displaystyle \sum _{i=1}^n{\rho }_{i,\Theta }^2/{\sigma }_{{\epsilon }_i}^2}} ({\rho }_{i,\Theta }^2>0,{\sigma }_{{\epsilon }_i}^2>0)\\ 1/3\hspace{1em}\hspace{1em}({\rho }_{i,\Theta }^2<0,{\sigma }_{{\epsilon }_i}^2<0)\end{array}\right.$$

(8)

$$Y_{model}=W_1{M}_1+W_2{M}_2+W_3{M}_3$$

(9)

where $M_i (i=1,2,3)$ denotes the time series of the same position for the three datasets and $W_i$ denotes the weight of the corresponding time series, $M_i$.

2.3.2. Estimation of GWSA Based on GRACE

Generally, GWSA can be estimated by subtracting CWS, SWE, and SM simulated by hydrologic models from the GRACE-derived TWS anomalies [4]. The detailed formula is as follows:

$$GWSA=TWSA-SMA-SWEA-CWSA$$

(10)

where $SMA$, $SWEA$, and $CWSA$ represent the storage anomalies of SM, SWE and CWS, respectively, relative to a reference period (the reference period is 2004–2009). Previous studies indicated that SM and GWS are the primary contributors to TWS changes and that variations in snow and ice, biomass, and surface water are relatively minor [9,19,58,59]. In addition, the selected study area is located in the arid region of northwest China and is mainly covered by bare land and gobi. Therefore, $SWE$ and $CWS$ can be ignored in this study.

2.3.3. Multiple Linear Regression of Time Series

To analyze the seasonal and secular trend of GWSA, multiple linear regression is used to analyze the temporal variability of GWSA. The regression model is given by [60]:

$$Y(t)={\beta }_1+{\beta }_{2}t+{\beta }_{3}sin(\pi t)+{\beta }_{4}cos(\pi t)+{\beta }_{5}sin(2\pi t)+{\beta }_{6}cos(2\pi t)+\epsilon$$

(11)

where $Y(t)$ denotes GWSA at time t; ${\beta }_1$ and ${\beta }_2$ denote the constant offset and secular trend, respectively; ${\beta }_3$ and ${\beta }_4$ represent the annual signal; ${\beta }_5$ and ${\beta }_6$ represent semi-annual signals; and $\epsilon$ represents the model error. Meanwhile, the annual and semi-annual amplitude are computed as [60]:

$$annual amplitude=\sqrt{{\beta }_3^2+{\beta }_4^2}$$

(12)

$$semi-annual amplitude=\sqrt{{\beta }_5^2+{\beta }_6^2}$$

(13)

2.3.4. Estimation of the Contribution to GWS

GRACE can monitor the temporal and spatial changes of TWS, including human factors and climate factors. In order to evaluate the contribution of different factors to GWS, a method is used to evaluate the contribution of these two factors, which can be computed as follows [61]:

$$GWS{C}_c=GWS{C}_{GRACE}-GWS{C}_H$$

(14)

$$GWSC=GWS{A}_t-GWS{A}_{t-1}$$

(15)

$${\eta }_H=\frac{GWS{C}_H}{\left|GWS{C}_H\right|+\left|GWS{C}_C\right|}$$

(16)

$${\eta }_C=\frac{GWS{C}_C}{\left|GWS{C}_H\right|+\left|GWS{C}_C\right|}$$

(17)

where $GWS{C}_C$ represents climate-driven GWS changes; $GWS{C}_{GRACE}$ represents the annual GWS changes estimated by GRACE data; and $GWS{C}_H$ denotes the part of GWS changes induced by human factors. ${\eta }_H$ and ${\eta }_C$ denote the contribution of human and climatic factors to GWS changes, respectively. If $\eta$ is positive, it provides a positive impact on GWS; otherwise, the opposite is true.

2.3.5. Evaluation Index

In this study, the correlation coefficient (CC), the root mean squared error (RMSE), and the Nash-Sutcliffe efficiency coefficient (NSE) are utilized to test the performance of this result [61,62,63].

$$CC=\frac{cov(X(t),Y(t))}{\sqrt{var[X(t)]var(Y(t))}}$$

(18)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(X(t)-Y(t))}^2}}$$

(19)

$$NSE=1-\frac{{\displaystyle \sum _{i=1}^n{(X(t)-Y(t))}^2}}{{\displaystyle \sum _{i=1}^n{(X(t)-X^{mean})}^2}}$$

(20)

where $n$ denotes the total number of observations; $X(t)$ and $Y(t)$ denote measurements and simulated values, respectively; and $X^{mean}$ represents the mean of measurements. Considering inconsistent scales between different results, in situ groundwater-level measurements and simulated results should be normalized to [−1, 1].

3. Results

3.1. Experimental Verifications of the GWFM

It is necessary to test the performance of the GWFM before it is applied to the study area and the SYRB is chosen as the study area. The GWSA results based on GRACE observations and three hydrological models (namely GLDAS, ERA5, and WGHM) are introduced as the input data of the GWFM, and GWFM-based results are verified against in situ groundwater-level measurements.

Figure 3 represents the time series of GWSA estimated from the GWFM and three GRACE-based GWSA (hereafter GRACE−GLDAS, GRACE−ERA5, and GRACE−WGHM). Moreover, in situ groundwater-level measurements are also shown, and the shaded areas represent the uncertainties of GWFM-based GWSA. From the long-term variation of GWSA point of view, GWFM-based GWSA agree well with that, based on GRACE in terms of periodicity and seasonality. Furthermore, it is obvious that the long-term trends of GRACE-based and GWFM-based GWSA and in situ groundwater-level measurements are generally similar, showing a decreasing trend. However, there is a clear difference in phase between them (Figure 3). In other words, there is a clear time lag between GRACE-based GWSA and in situ groundwater-level measurements. Many previous studies have reported the time lag; for example, Thomas et al. [64] indicated that when the lag time was two months, the correlation between GRACE-based GWSA and in situ groundwater-level measurements reached a maximum in the Central Valley of California. Abou et al. [65] reported that there was a clear time lag between in situ groundwater-level measurements and GWSA based on GRACE in the Bakhtegan catchment. In order to explore the best lag time, this study uses GRACE−GLDAS and in situ groundwater-level measurements (S1–S5) (

Table 2. Lagged CC between in situ groundwater-level measurements and GWSA in S1–S5 from 2007 to 2014.

Wells	Time Lag (Month)
	0	1	2	3	4	5	6
S1	0.34	0.41	0.51	0.60	0.65	0.63	0.55
S2	0.31	0.38	0.51	0.61	0.68	0.72	0.70
S3	0.59	0.60	0.64	0.66	0.69	0.67	0.62
S4	0.48	0.49	0.52	0.54	0.58	0.59	0.57
S5	0.51	0.53	0.58	0.62	0.64	0.61	0.56

). The result shows that the highest correlation (CC = 0.59–0.72) can be found when the lag time is 4–5 months.

In order to explore the reliability of the GWFM, the lag time is set to 4 months, and GRACE-based (GRACE−GLDAS, GRACE−ERA5, and GRACE−WGHM) and GWFM-based results are compared with in situ groundwater-level measurements (Figure 4). The comparison indicates that the seasonality of GRACE-based and GWFM-based GWSA is consistent with in situ groundwater-level measurements, and the annual amplitude of GRACE−ERA5 is greater than that of GRACE−GLDAS and GRACE−WGHM. Furthermore, the amplitudes of in situ groundwater-level measurements (S1–S5) also display larger differences. For example, S3 shows a small amplitude change after 2011; the amplitude of S1 varies from −3 m to 5 m, and the amplitude of S5 is between −8 m and 16 m from the perspective of the long-term average.

In order to quantify the agreement between GWFM-based GWSA and in situ groundwater-level measurements, three metrics (CC, RMSE, and NSE) are calculated over the SYRB, as shown in Figure 5. This result shows that the agreement between GWFM-based GWSA and in situ groundwater-level measurements is much better than that based on GRACE. Additionally, the GWFM effectively improves the CC and NSE and decreases the RMSE. Specifically, the CC between GWFM-based GWSA and in situ groundwater-level measurements, S1, increases from 0.54 to 0.74, and the RMSE decreases from 0.44 to 0.39, the NSE increases from 0.11 to 0.54 relative to the original results. Compared with the mean value of in situ groundwater-level measurements (expressed by S6), similar improvements can also be seen for CC (9–40%), NSE (23–657%), and RMSE (9–28%). The above verification results indicate that reasonable GWSA estimates can be obtained through the GWFM in the SYRB. Therefore, it can give us confidence in applying this developed method to the HC, so as to better understand the temporal and spatial characteristics.

3.2. Comparison of GWSA

Figure 6 shows the annual, monthly, and seasonal scales of GRACE-based and GWFM-based GWSA from January 2003 to December 2016. The long-term trend of GRACE-based and GWFM-based GWSA shows a reasonable agreement, and four results also have a similar annual cycle. For the intra-annual changes of GWSA, the GWSA time series have a reasonable agreement; the anomalies are positive from May to August (Figure 6c). However, the values are negative for other months, and the only exception is the GRACE−GLDAS, always remaining negative. The main reason is that the HC is dry and less rainy, and ~80% of rainfall occurs during the period from May to September [66], which effectively recharges the groundwater. The HC is a well-known irrigated agricultural area in northwest China, but surface-water resources are scarce, and irrigation water mainly comes from groundwater. Therefore, a large amount of groundwater is pumped in spring and summer due to irrigation needs, which leads to a decrease in groundwater storage.

From the perspective of the long-term trend, GWSA reveals a significant downward trend over the study period. Notably, 2011 is a turning point, and the downward trend before 2011 is significantly higher than the trend after 2011. In order to clarify the difference between the four results, multivariate statistical analysis is used for the GWSA time series in the two time periods of 2003–2010 and 2011–2016 in the HC.

Table 3. Comparison of annual amplitude and trend between GWSA from GRACE−GLDAS, GRACE−ERA5, GRACE−WGHM, and GWFM.

Datasets	2003–2010		2011–2016
	Annual Amplitude (mm)	Trend (mm/yr)	Annual Amplitude (mm)	Trend (mm/yr)
GRACE−GLDAS	6.59 ± 1.54	−4.17 ± 0.47	7.07 ± 3.34	−0.48 ± 1.38
GRACE−ERA5	9.44 ± 1.42	−1.08 ± 0.44	10.10 ± 3.37	−0.70 ± 1.39
GRACE−WGHM	7.12 ± 1.16	−1.67 ± 0.36	9.72 ± 3.20	0.07 ± 1.32
GWFM	7.40 ± 1.22	−2.37 ± 0.38	8.67 ± 3.26	−0.46 ± 1.35

summarizes the evaluation indexes for the GWSA time series, including the long-term trend and annual amplitude. From 2003 to 2010, four GWSA downtrends range from 1.08 (GRACE−ERA5) to 4.17 (GRACE−GLDAS) mm/yr, which clearly indicates groundwater depletion in the HC during the period from 2003 to 2011. From 2011 to 2016, the results also show a downward trend (except GRACE−WGHM), and decline rates range from 0.46 (GWFM) to 0.70 (GRACE−ERA5) mm/yr, while GRACE−WGHM increased at a rate of 0.07 mm/yr. The rate of decline from 2011 to 2016 is significantly slower than that from 2003 to 2010, while the annual amplitude, compared with the previous period, increases significantly. This result may be related to the water policy in the area, such as Gansu Province gradually implementing the most stringent water-resource management system and measures of “points to areas, Hexi first” in 2011 [67].

The annual and semi-annual changes of time series can be analyzed by phasor diagrams, which show their amplitude and phase based on a reference period (in this study, the reference period is 2004–2009). The length of each vector represents the magnitude of amplitude, while the vector direction represents a phase. The bigger the difference between two vector directions and length, the greater the phase and amplitude difference between two time series. In other words, the phase difference can affect the magnitude of the correlation, and a difference in amplitude can affect variance agreement. In the case of annual amplitude and phases, GRACE−ERA5 has a higher annual amplitude than other results (shown in Figure 7a). Although the best amplitude agreement exists between GWFM-based result and GRACE−WGHM, the phase correspondence is poor, while GWFM-based GWSA agree well with GRACE−ERA5 at phase. In terms of semi-annual amplitude and phases (shown in Figure 7b), the phase and amplitude in GWSA from the GWFM show favorable agreement with GRACE−GLDAS, while the semi-annual phase and amplitude of GRACE−ERA5 are different from other results.

3.3. Spatial Pattern of Variation Trends in GWSA

Figure 8 shows the spatial distribution of GWSA based on GRACE and GWFM over the HC from 2003 to 2016. Among these results, GRACE−GLDAS shows that the area of GWS depletion has a higher downward trend and coverage, but the characteristics of spatial distribution in the Shule River Basin (SLRB) and the HRB are not distinct. GRACE−ERA5 shows obvious spatial-change characteristics, such as the D1 of the SLRB, the D2 of the HRB, and the D3 of the SYRB as the main GWS depleted areas (shown in Figure 8b). GRACE−WGHM shows a downward trend, high in the north and low in the south, but the overall spatial distribution shows no significant characteristic changes relative to other results. The GWFM highlights remarkable GWS depletion in the d1, d2, d3, and d4 (shown in Figure 8d), with the rates of about −4.22 mm/yr, −2.67 mm/yr, −3.77 mm/yr, and −5.06 mm/yr, respectively. It should be noted that, GRACE−ERA5 and the GWFM show opposite trends in the southeast of the SYRB, with the rates of 1.91 and −1.76 mm/yr, respectively.

The GWSA based on the GWFM is shown in Figure 9a, and the simple average result (average) of both GRACE-based GWSA is shown in Figure 9b. The main depletion areas of the two results are basically similar, but there is a large difference between the GWFM and average in the southeast of the SYRB (Df1 and df1), where the two trends are −1.76 mm/yr and 1.11 mm/yr, respectively. The government of Gansu Province announced Gulang County and Wuwei City as a GWS over-exploitation area, that is, the Df4 in Figure 9a, which is more consistent with the results of the GWFM. Therefore, the GWFM can effectively integrate the advantages of multiple models, retain the characteristics of specific regional changes, and provide a more accurate GWSA result relative to simple average results.

3.4. Response of GWSA to Climate Change

GWS changes are closely related to climate and human factors. Therefore, it is necessary to evaluate the relationship between GWFM-based GWSA and influencing factors in order to better understand the causes and development of groundwater depletion in the HC.

In the context of climate change, precipitation and evapotranspiration are the dominant factors that have the greatest impact on GWS [68]. Temperature changes lead to changes in evapotranspiration, which, in turn, lead to changes in GWS. The net recharge of groundwater is the difference between recharge and discharge [2]. Groundwater in the HC piedmont plain is mainly from the infiltration of surface runoff, which accounts for ~80% of total recharge [69], followed by precipitation and seepage of irrigation water. Groundwater discharge is mainly groundwater pumping and evapotranspiration of shallow groundwater [70].

To further analyze the detailed relationship between climatic factors and GWSA, cross-wavelet analysis is used in this study. Cross-wavelet transforms between GWSA and climatic factors in the HC are displayed in Figure 10. Figure 10a indicates that the correlations between precipitation and GWSA are strong in the HC during the period of 2003–2016, and it shows a statistically positive correlation between precipitation and GWSA at the 95% confidence level. Evapotranspiration and temperature also exhibit a strong positive correlation. In addition, GWSA and climatic factors all have a main resonance period of about 1 month.

Under the background of climate change, precipitation is the input of water, and evapotranspiration is the output of water in a region. Therefore, precipitation minus evapotranspiration (P−ET) can represent the net recharge of surface water and groundwater [71]. During 2003–2016, the maximum P−ET occurred in June-September, and the minimum in January-April and October-December (shown in Figure 11). P−ET shows a significant downward trend relative to other time periods during 2007–2016, resulting in a significant decrease in net recharge. Such shortage of precipitation will directly hinder the growth of vegetation and human production, and excessive evapotranspiration will further accelerate the loss of available water resources and disrupt the balance of the water cycle [72]. Although the net recharge in summer is positive during this period, groundwater is still in a state of declining, indicating that groundwater is not effectively recharged.

In addition, snowmelt is also an important factor in groundwater replenishment. In the context of climate change, snowmelt will increase. This impact means less snow accumulation in the winter and an earlier peak runoff in the spring [73]. Meanwhile, snowmelt is an important water source in Northwest China, which is of great significance to maintenance of ecological balance and sustainable development [74]. Li et al. [75] showed that from 1960 to 2010, the average annual runoff in the arid area of northwest China was increasing. Among them, the increased rate of runoff in the northern mountainous area of the Qilian Mountains was 1.48 × 108 m³/10 a. Therefore, an increase in snowmelt will have a greater impact on runoff, which, in turn, affects recharge of groundwater.

3.5. Response of GWSA to Human Factors

In addition to climate factors, the impact of human factors on groundwater cannot be ignored. In the HC, water resources are scarce and unevenly distributed. During the crop-water demand season, a large amount of water resources is used for irrigation, when surface water for irrigation is limited, which will lead to a prominent contradiction between water supply and demand and inevitably lead to groundwater depletion. The HC has been undergoing tremendous changes over the past few decades. Niu et al. [76] showed that in the past 30 years, the increase in irrigation water consumption of farmland in the HRB had led to an average drop of about 1.86 m in groundwater. Zhou et al. [77] showed that as the area of farmland increased by 11.0%, the total irrigation water demand increased by 6.3% during the period from 2000 to 2010.

Figure 12 shows changes in groundwater withdrawal in the HC. During the entire survey period (2003–2016), the amount of groundwater withdrawal in the HRB is on a continuous upward trend, while other regions show a trend of a first decline and then an increase, reaching the lowest in 2009–2011, which is consistent with the change in precipitation. It is worth noting that even in the rainy season, GWS is in a state of decline. However, even in drier years, the amount of groundwater withdrawal is lower than in previous years, and the rate of decline in GWSA is higher than in the rainy season. This means that even when the amount of groundwater withdrawn is lower than in previous years, the amount of groundwater withdrawn is much higher than the amount of groundwater replenishment. If the use of water resources cannot be well improved, the area will continuously face the problem of groundwater depletion in the future.

4. Discussion

4.1. Spatial Distribution of Weight Index

The weighted fusion model is presented in this study, which can merge three GRACE-based GWSA. Specifically, the error variance and correlation coefficient of three GWSA derived by ETC are used for weight estimation. Then, three GWSA from different sources are merged by the least-squares framework. Therefore, it is necessary to discuss the weight of different original results.

Figure 13 shows the spatial distribution of the weights. This weight represents the relative contribution to the merged result. Among the three GRACE-based GWSA (including GRACE−GLDAS, GRACE−ERA5, and GRACE−WGHM), the largest average weight can be obtained by GRACE−ERA5 (0.38), followed by GRACE−WGHM (0.32) and GRACE−GLDAS (0.30). It is worth noting that there are apparent differences in spatial distribution, although the average weights of GRACE−GLDAS and GRACE−WGHM are relatively close. For example, GRACE−WGHM has larger weights than GRACE−GLDAS in C1, C2, and C3, accounting for about 0.50, 0.45, and 0.49, respectively. GRACE−GLDAS matches well with GRACE−WGHM in other regions. GRACE−ERA5 has higher weights relative to other results in B1 and B2. Furthermore, GRACE−ERA5 has a relatively high relative contribution to the merger result in most regions. These differences may be caused by different forced data and the parameters of the hydrological model [7,78]. In general, GRACE−ERA5 can accurately describe GWS changes in most areas of the HC from the perspective of single-model results.

4.2. Contributions of Different Factors to GWS

The contribution of climate factors and human factors to GWS has been evaluated using the method proposed in Section 2.3.4. Generally, the larger the value of η, the greater contribution it makes to GWS changes. The respective contributions are shown in

Table 4. Contributions (%) of climate factors and human factors to GWS changes in the HC and its three subregions.

	HC		SLRB		HRB		SYRB
	ηC	ηH	ηC	ηH	ηC	ηH	ηC	ηH
2003	53.73	−46.27	81.23	−18.77	37.28	−62.72	47.37	−52.63
2004	51.86	−48.14	27.78	−72.22	63.97	−36.03	71.66	−28.34
2005	−15.23	−84.77	−69.58	−30.42	−88.34	−11.66	−70.40	−29.60
2006	63.40	−36.60	75.86	−24.14	84.24	−15.76	73.67	−26.33
2007	−40.64	−59.36	−85.64	−14.36	−90.70	−9.30	−42.28	−57.72
2008	32.43	−67.57	−38.21	−61.79	12.85	−87.15	−58.48	−41.52
2009	65.73	−34.27	74.17	−25.83	81.21	−18.79	84.84	−15.16
2010	−61.61	−38.39	−87.49	−12.51	−87.27	−12.73	−90.15	−9.85
2011	70.55	−29.45	90.63	−9.37	84.20	−15.80	71.01	−28.99
2012	−11.98	−88.02	−88.78	−11.22	15.48	−84.52	−39.99	−60.01
2013	56.40	−43.60	62.68	−37.32	54.32	−45.68	80.03	−19.97
2014	55.33	−44.67	78.58	−21.42	46.74	−53.26	47.75	−52.25
2015	54.91	−45.09	78.35	−21.65	64.81	−35.19	−30.96	−69.04
2016	42.46	−57.54	−14.58	−85.42	8.81	−91.19	15.24	−84.76

.

As shown in Table 4, human factors are an important factor affecting GWS among climate and human factors in the HC and its subregions. For example, the impact of human activities shows a downward trend before 2010 and then begin to rise in the SYRB, which matches well with the change in trend of groundwater withdrawal in Figure 12. This is consistent with the conclusion drawn from other research conducted in this region, e.g., Liu et al. [44]. As for HRB, GWS changes affected by human factors show an upward trend with the increase in groundwater withdrawal. In the HC, the effects of climate change on GWS changes account for ~48%, while those of human activities contributed ~52%. This indicates that human activity has been the dominant factor driving the continuous reduction in groundwater. Wang et al. [5] reported that irrigation was continuously increased during the period of 2000–2016 in the HC. Moreover, Niu et al. [76] and Zhou et al. [77] also reported a similar situation in the subregions of the HC. This is consistent with our conclusion that human factors have become the dominant factor affecting GWS in the HC.

4.3. Limitation and Furture Work

In this study, there are some limitations, although reasonable results have been achieved. First, the GWFM proposed in this study is based on the same assumptions as the ETC method, such as error orthogonality and zero-error cross-correlation. Second, this study only uses GLDAS, ERA5, and WGHM, without considering other data sources. Third, there is only an eight-year overlap period with the GRACE data due to the short and discontinuous in situ groundwater-level measurement data. It is impossible to test the performance of the GWFM over a longer time frame.

In follow-up work, more data will be introduced to develop the model to compare and analyze GWSA in arid regions where in situ groundwater-level measurements are scarce. Furthermore, we will conduct a comparative analysis with the existing multi-model combination of technology to further improve the GWFM.

5. Conclusions

This study estimates GWSA by combining GRACE data with hydrological models. These data have their own unique characteristics and are developed for a global scale. Therefore, estimation of regional water-storage state generated by a single model may have greater uncertainty. In response to this problem, a GWFM is presented that can merge multi-source GWSA. The useful conclusions are as follows:

(1)

To obtain an accurate estimation of GWSA, this paper proposes a groundwater weighted fusion model. A comprehensive example is defined to verify the performance of the GWFM, and the superiority of the GWFM is verified by in situ groundwater-level measurements. The results show that the GWFM can effectively integrate the advantages of each data set sand produce a more reliable GWSA than the original results. Compared with GRACE-based GWSA, GWFM-based GWSA can obtain higher CC and NSE, CC increases by 9–40%, NSE increases by 23–657%, while RMSE decreases by 9–28%.

(2)

The GWSA result of the HC from 2003 to 2016 is calculated based on the GWFM. GWFM-based GWSA show an overall downward trend from 2003 to 2016, but 2011 is a turning point. From 2003 to 2010, there is a rapid downward trend, which is −2.37 ± 0.38 mm/yr, while the downward trend from 2011 to 2016 is significantly slowed, at −0.46 ± 1.35 mm/yr. This may be related to the local implementation of corresponding water-saving policies. In terms of spatial changes, in the central and southern part of the SLRB, the central part of the HRB and the northern part of the SYRB, which are the main GWS depleted areas, have a large downward trend. Furthermore, GWFM-based GWSA can better retain the characteristics of regional GWSA relative to the arithmetic average result, especially in the southeast of the SYRB.

(3)

A simple and effective method is used to evaluate the contribution of climate factors and human factors to GWS. The results show that the amount of groundwater withdrawal has a significant impact on GWS, especially in the HRB, where the amount of groundwater withdrawal is increasing every year. As for the HC, the effects of climate change on GWS changes account for ~48%, while those of human activities contributed ~52%. In general, human activities, especially agricultural irrigation, have become the main reason for GWS decline in the HC.