*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.

**Figure 2.** The flow chart of the study.

#### 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\_{\bar{i}} = \alpha\_{\bar{i}} + \beta\_{\bar{i}} \Theta + \varepsilon\_{\bar{i}} \tag{1}$$

where *Si* (*i* = 1, 2, 3) represents GWSA based on GRACE and three hydrological models; Θ denotes the unknown true hydrological signal; *Si* represents collocated measurement systems linearly related to the true unknown value, Θ; *α<sup>i</sup>* and *β<sup>i</sup>* represent the least-squares intercepts and slope, respectively; and *ε<sup>i</sup>* 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} \text{cov}(\mathcal{S}\_{i\boldsymbol{\nu}}\mathcal{S}\_{\boldsymbol{j}}) = E(\mathcal{S}\_{i}\mathcal{S}\_{\boldsymbol{j}}) - E(\mathcal{S}\_{i})E(\mathcal{S}\_{\boldsymbol{j}}) = \beta\_{i}\beta\_{\boldsymbol{j}}\sigma\_{\boldsymbol{\Theta}}^{2} \\ + \beta\_{i}\text{cov}(\boldsymbol{\Theta}, \boldsymbol{\varepsilon}\_{\boldsymbol{j}}) + \beta\_{\boldsymbol{j}}\text{cov}(\boldsymbol{\Theta}, \boldsymbol{\varepsilon}\_{\boldsymbol{i}}) + \text{cov}(\boldsymbol{\varepsilon}\_{i\boldsymbol{\nu}}\boldsymbol{\varepsilon}\_{\boldsymbol{j}}) \end{array} \tag{2}$$

According to these four prerequisites of the TC approach, cov(*εi*,*εj*) = 0 (*i* = *j* , cov(*εi*, Θ) = 0; the equation reduces to [35,40]:

$$\mathbf{C}\_{\mathbf{i}\mathbf{j}} = \text{cov}(\mathbf{S}\_{\mathbf{i}}, \mathbf{S}\_{\mathbf{j}}) = \begin{cases} \begin{array}{c} \beta\_{\hat{i}} \beta\_{\hat{j}} \sigma\_{\Theta}^{2} \left( \mathbf{i} \neq \mathbf{j} \right) \\ \beta\_{\hat{i}}^{2} \sigma\_{\Theta}^{2} + \sigma\_{\varepsilon\_{i}}^{2} \left( \mathbf{i} = \mathbf{j} \right) \end{array} \tag{3}$$

where *σ*<sup>2</sup> *<sup>ε</sup><sup>i</sup>* represents the variance of random-error variance, *<sup>ε</sup>i*; *<sup>β</sup>*<sup>2</sup> *i σ*2 <sup>Θ</sup> denotes the sensitivity of datasets, *Si*, to changes in true signal. In other words, the higher *βi*, the stronger the response of datasets,*Si*, to hydrological signal. The sensitivity of each dataset can be calculated by combining their covariances [35,40]:

$$
\beta\_i^2 \sigma\_{\Theta}^2 = \begin{bmatrix}
\frac{\underline{C\_{12}} \underline{C\_{13}}}{\underline{C\_{23}}} \\
\frac{\underline{C\_{21}} \underline{C\_{23}}}{\underline{C\_{12}}} \\
\frac{\underline{C\_{31}} \underline{C\_{32}}}{\underline{C\_{12}}}
\end{bmatrix} \tag{4}
$$

The error variance can be obtained by subtracting the sensitivity *β*<sup>2</sup> *i σ*2 <sup>Θ</sup> of each dataset from their total variance [35,40]:

$$
\sigma\_{\varepsilon\_i}^2 = \begin{bmatrix}
\mathbf{C}\_{11} - \frac{\mathbf{C}\_{12}\mathbf{C}\_{13}}{\mathbf{C}\_{23}} \\
\mathbf{C}\_{22} - \frac{\mathbf{C}\_{21}\mathbf{C}\_{23}}{\mathbf{C}\_{12}} \\
\mathbf{C}\_{33} - \frac{\mathbf{C}\_{31}\mathbf{C}\_{32}}{\mathbf{C}\_{12}}
\end{bmatrix} \tag{5}
$$

$$\rho\_{i,\Theta}^2 = \frac{\beta\_i^2 \sigma\_{\Theta}^2}{\beta\_i^2 \sigma\_{\Theta}^2 + \sigma\_{\varepsilon\_i}^2} = \frac{SNR\_i}{SNR\_i + 1} \tag{6}$$

$$SNR\_i = \frac{\text{var}(S\_i)}{\text{var}(\varepsilon\_i)} = \frac{\beta\_i^2 \sigma\_{\Theta}^2}{\sigma\_{\varepsilon\_i}^2} \tag{7}$$

where *ρ*<sup>2</sup> *<sup>i</sup>*,<sup>Θ</sup> 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:

$$\mathcal{W}\_{i} = \begin{cases} \frac{\rho\_{i,\Theta}^{2}/\sigma\_{\varepsilon\_{i}}^{2}}{\frac{\pi}{n}\rho\_{i,\Theta}^{2}/\sigma\_{\varepsilon\_{i}}^{2}} & (\rho\_{i,\Theta}^{2} > 0, \sigma\_{\varepsilon\_{i}}^{2} > 0) \\\ 1/3 & (\rho\_{i,\Theta}^{2} < 0, \sigma\_{\varepsilon\_{i}}^{2} < 0) \end{cases} \tag{8}$$

$$Y\_{\text{model}} = \mathcal{W}\_1 \mathcal{M}\_1 + \mathcal{W}\_2 \mathcal{M}\_2 + \mathcal{W}\_3 \mathcal{M}\_3 \tag{9}$$

where *Mi* (*i* = 1, 2, 3) denotes the time series of the same position for the three datasets and *Wi* denotes the weight of the corresponding time series, *Mi*.

#### 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 \tag{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) + \varepsilon \tag{11}$$

where *Y*(*t*) denotes GWSA at time *t*; *β*<sup>1</sup> and *β*<sup>2</sup> denote the constant offset and secular trend, respectively; *β*<sup>3</sup> and *β*<sup>4</sup> represent the annual signal; *β*<sup>5</sup> and *β*<sup>6</sup> represent semi-annual signals; and *ε* represents the model error. Meanwhile, the annual and semi-annual amplitude are computed as [60]:

$$
gammaal\ amplitude = \sqrt{\beta\_3^2 + \beta\_4^2} \tag{12}$$

$$\text{Isemi} - \text{annual amplitude} = \sqrt{\beta\_5^2 + \beta\_6^2} \tag{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]:

$$GWSC\_c = GWSC\_{GRACE} - GWSC\_H \tag{14}$$

$$GWSC = GWSA\_t - GWSA\_{t-1} \tag{15}$$

$$\eta\_H = \frac{GWSC\_H}{|GWSC\_H| + |GWSC\_C|} \tag{16}$$

$$\eta\_{\mathbb{C}} = \frac{GWSC\_{\mathbb{C}}}{|GWSC\_{H}| + |GWSC\_{\mathbb{C}}|} \tag{17}$$

where *GWSCC* represents climate-driven GWS changes; *GWSCGRACE* represents the annual GWS changes estimated by GRACE data; and *GWSCH* denotes the part of GWS changes induced by human factors. *η<sup>H</sup>* and *η<sup>C</sup>* denote the contribution of human and climatic factors to GWS changes, respectively. If *η* 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–63].

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

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( X(t) - \mathbf{Y}(t) \right)^{2}} \tag{19}$$

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

where *n* denotes the total number of observations; *X*(*t*) and *Y*(*t*) denote measurements and simulated values, respectively; and *Xmean* 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 GRACEbased 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). 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 GWFMbased 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.

**Figure 3.** Comparisons of GWSA time series from GRACE−GLDAS, GRACE−ERA5, GRACE−WGHM, and GWFM and verified against in situ groundwater-level measurements.

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


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.

**Figure 4.** Comparison of GWSA from GRACE−GLDAS, GRACE−ERA5, GRACE−WGHM, and GWFM estimates with in situ groundwater-level measurements after setting the lag time. (**a**) S1; (**b**) S2; (**c**) S3; (**d**) S4; (**e**) S5; (**f**) average of S1–S5.

**Figure 5.** Comparison of the evaluation index between GWSA and in situ groundwaterlevel measurements.

### *3.2. Comparison of GWSA*

Figure 6 shows the annual, monthly, and seasonal scales of GRACE-based and GWFMbased 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.

**Figure 6.** Comparison of GRACE-based and GWFM-based GWSA on different time scales over the HC during the period from 2003 to 2016. (**a**) Monthly; (**b**) annual; (**c**) seasonal.

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 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].

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


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 GWFMbased 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.

**Figure 7.** Vector diagram of the GWSA amplitude and phase change from 2003 to 2016. (**a**) Annual; (**b**) semi-annual.

#### *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.

**Figure 8.** The spatial distribution of GRACE-based and GWFM-based GWSA over the HC during the period from 2003 to 2016. (**a**) GRACE−GLDAS; (**b**) GRACE−ERA5; (**c**) GRACE−WGHM; (**d**) GWFM.

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.

**Figure 9.** Comparison of the spatial-trend distribution in the Hexi corridor between the average results of original results (GRACE−GLDAS, GRACE−ERA5, and GRACE−WGHM) and GWFMbased GWSA. (**a**) GWFM; (**b**) average.

#### *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.

**Figure 10.** Cross-wavelet transforms between GWSA and climatic factors at monthly scale in the HC. (**a**) Precipitation; (**b**) evapotranspiration; (**c**) temperature.

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.

**Figure 11.** Comparison between P–ET and GWFM-based GWSA over the HC during the period from 2003 to 2016.

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 m3/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.

**Figure 12.** Time series of groundwater withdrawal in the HC and its three subregions.

#### **4. Discussion**

#### *4.1. Spatial Distribution of Weight Index*

The weighted fusion model is presented in this study, which can merge three GRACEbased 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.

**Figure 13.** Weights of three GRACE-based GWSA for the fused GWSA. (**a**) GRACE−GLDAS; (**b**) GRACE−ERA5; (**c**) GRACE−WGHM.

#### *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.

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


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.
