Assessment of Variability and Attribution of Drought Based on GRACE in China from Three Perspectives: Water Storage Component, Climate Change, Water Balance

Abstract

Understanding how drought is impacted by both natural and human influences is crucial to the sustainable utilization and protection of water resources. We established a drought severity index (DSI) based on the terrestrial water storage anomaly (TWSA) derived from the GRACE satellite to detect drought characteristics and trends over ten major river basins in China from 2002 to 2017. The influence of natural factors (terrestrial water storage components, precipitation, evapotranspiration, runoff, NDVI, and teleconnection factors (ENSO, PDO, NAO, and AO)) and a human factor (LULC) on drought were investigated and quantified from the perspective of water storage components based on the Theil–Sen trend and Mann–Kendall test method, the perspective of climate change based on cross wavelet transforms, and the perspective of water balance based on Random Forest. The results indicated that (1) almost all humid and arid basins experienced major drought periods during 2002–2006 and 2014–2017, respectively. The southern IRB and central YZRB regions exhibited notable declines in DSI trends, while the majority of the HLRB, IRB, LRB, YRB, HRB, and SWRB experienced significant increases in DSI trends; (2) abnormal groundwater decreases were the main cause of drought triggered by insufficient terrestrial water storage in most basins; (3) ENSO was the strongest teleconnection factor in most humid basins, and NAO, PDO, and AO were the strongest teleconnection factors in the arid basins and PRB. Most significant resonance cycles lasted 12–64 months in 2005–2014; and (4) the influence of an anthropogenic driver (LULC) has become as important as, or more important than, natural factors (runoff and teleconnection factors) on hydrological drought.

1. Introduction

Approximately 80% of the world’s population experiences acute water scarcity or water insecurity as a result of climate change and anthropogenic activities [1,2,3]. As a result, drought has become the most dreaded natural disaster [4]. Characterized by slow and widespread development, recurrence, and high cost [5,6], drought is a frequent and serious natural disaster in China [7,8]. Severe and extreme droughts have become more serious since the late 1990s [9] and pose a challenge to the sustainable utilization and protection of water resources. Therefore, it is crucial to find the dominant factors behind China’s droughts.

Due to the multiple hydrometeorological drivers of drought at different spatial and temporal scales, it is challenging to comprehensively characterize droughts [10]. The GRACE (Gravity Recovery and Climate Experiment) satellite provides an alternative approach to monitor drought from an integrated perspective, which can detect vertically integrated terrestrial water storage (TWS) changes from the land surface to the deepest aquifers [11]. Since its launch in 2002, the GRACE satellite has been used in numerous drought studies. On one hand, GRACE data were used for drought monitoring by creating drought indices. For example, Zhao et al. [12] proposed the DSI (Drought Severity Index), which can compare drought characteristics in different regions and periods, and incorporates changes in water storage caused by human influence. Liu et al. [11] successfully used DSI based on detrended GRACE TWS to reflect droughts driven only by climate factors in major river basins in China. Sinha et al. [13] developed the CCDI (Combined Climatologic Deviation Index), which covers all aspects of meteorological, agricultural, hydrological, and anthropogenic drought occurrence, and was successfully applied to drought monitoring in China [14]. The drought-monitoring abilities of other indices, such as the TSDI (Total Storage Deficit Index) [15] and WSDI (Water Storage Deficit Index) [6], have been demonstrated in regions of China [16,17]. In addition, GRACE-based groundwater and soil moisture drought indicators were published by NASA’s Goddard Space Flight Center at https://nasagrace.unl.edu/About.aspx (accessed on 22 August 2023), which were derived from GRACE satellite data and other observations using a sophisticated numerical model of land-surface water and energy processes [18,19,20].

On the other hand, many studies have analyzed the potential factors of drought based on GRACE TWS from the perspective of water storage components. The TWSA (terrestrial water storage anomaly) derived from GRACE is the sum of water storage anomalies in groundwater, soil moisture, surface, ice and snow, plant canopy, and biological storage [21]. If one component significantly contributes to the reduction of TWS compared to others, it can be considered the dominant factor in causing drought due to insufficient TWS. By utilizing the soil moisture, canopy water storage, total runoff, and snow water equivalent provided by a hydrological model such as GLDAS, it is possible to separate groundwater storage from GRACE TWSA [22]. Excessive groundwater extraction in the North China Plain [23], Northwest India [24], and the northern Middle East [25] has resulted in a considerable depletion of groundwater storage, consequently inducing drought conditions. Long et al. [26] presented that the changes of soil moisture storage can explain 70–80% of TWS depletion during a drought in 2011 in Texas. Meng et al. [27] conducted TWS components in the Tibetan Plateau and observed that the increase in TWS in the upstream of Yangtze and Yellow River Basins was primarily attributed to an augmentation in soil moisture. Furthermore, they found that the decrease in TWS in the Brahmaputra Basin was predominantly associated with glacier mass loss.

Numerous studies have explored dominant drought drivers from the perspective of water balance, which consider both natural climate variability and human activity [28,29]. According to the water balance equation P–E–R = TWSA (i)–TWSA (i–1) = TWSC (terrestrial water storage change, obtained by the difference in TWSA between month (i) and (i–1)), TWSA is influenced by a combination of precipitation (P), evapotranspiration (E), and runoff (R) [28,30,31]. Moreover, with the increasing human activity influence, the direct impact of land use and land cover (LULC) change on water resources cannot be ignored [29,32,33]. LULC has an important impact on water cycle processes by changing the properties of the underlying surface and affecting the exchange of heat and water between the land surface and the atmosphere [34]. Social and economic development have led to the conversion of large areas of forest and grassland into agricultural land, and population growth has led to increased urbanization [33]. Cui et al. [35] found that drought events in temperate monsoon and subtropical monsoon climates were mainly influenced by precipitation, while drought events in temperate continental and plateau mountain climates were mainly influenced by evapotranspiration. Deng et al. [36] found that human water use, especially for agriculture, was highly correlated with WSDI in the arid and semi-arid regions of North China and the Loess Plateau. Zhu and Zhang [37] identified that natural net recharge (P–E) and anthropogenic factors (LULC) are crucial factors affecting GRACE-based groundwater drought in the Yangtze and Yellow River basins in China, where groundwater exploitation plays a major role. Lv et al. [38] quantitatively analyzed TWS in the Yellow River basin of China under the Grain for Green project with a consideration for irrigation, which revealed that the increasing vegetation coverage and intensive irrigation was of great importance to TWS.

Teleconnection factors are strong signals derived from large-scale ocean–atmosphere interactions in global change [23,29,39], which are closely related to drought caused by the abnormal reduction in TWS. Hence, many researchers focus on the causes of drought from the perspective of climate change. The El Niño–Southern Oscillation (ENSO) and North Atlantic Oscillation (NAO) cause different degrees of interference to the normal operation of atmospheric circulation, and Cui et al. [35] detected that the affected area is prone to drought in China. Zhong et al. [40] proposed that the subsequent ENSO after 1997 induced a decade-long drought period in the west Liaohe River Basins. Yu et al. [9] proposed that the Arctic Oscillation (AO) and Pacific Decadal Oscillation (PDO) are important factors affecting groundwater drought in the Yangtze River Basin and Yellow River Basin, respectively. Wang et al. [23] pointed out that ENSO, AO, and PDO are associated with groundwater drought in North China Plain.

At present, there are limited comprehensive studies that have simultaneously analyzed the variability and attribution of drought using GRACE TWS data in China, considering the three perspectives mentioned above. Moreover, few studies have compared the differences among major river basins in China. These aspects are crucial for scientific and rational water resource management, decision-making, and mitigation of drought risks. In order to fill these gaps, this study had the following main objectives: (1) identify drought characteristic and trends using the Drought Severity Index (DSI) for the period from April 2002 to June 2017; (2) find the main components affecting drought from the perspective of water storages components; (3) analyze the relationship between teleconnection factors (ENSO, PDO, NAO, and AO) and drought from the perspective of climate change; and (4) quantitatively compare the attribution of natural factors and human factors to the development of drought, from the perspective of water balance.

2. Materials and Methods

2.1. Study Area

China (73°33′E–135°05′E, 3°51′N–53°33′N) is in the eastern part of the Asian and European continent, west of the Pacific Ocean. The spatial distribution of annual precipitation decreases from southeast to northwest, and the average annual temperature decreases from south to north, with obvious seasonal changes.

For the convenience of the study, only the main land areas of China were covered in this study, excluding areas such as islands and reefs in the South China Sea. According to the spatial division scheme for water resource basins from the Ministry of Water Resources of the People’s Republic of China (http://www.mwr.gov.cn/) (accessed on 10 August 2021), the study area was divided into Songhua River Basin (SRB), Liao River Basin (LRB), Haihe and Luanhe Rivers Basin (HLRB), Huaihe River Basin (HRB), Yellow River Basin (YRB), Yangtze River Basin (YZRB), Southeast Rivers Basin (SERB), Pearl River Basin (PRB), Southwest Rivers Basin (SWRB), and Inland Rivers Basin (IRB) (Figure 1). In the following study, basins were divided into arid (SRB, LRB, HLRB, HRB, YRB, and IRB) and humid (YZRB, SERB, PRB, and SWRB) basins based on annual average precipitation, characterized by red and blue contour of basins or font color for basins’ names in Figure 1.

2.2. Data

2.2.1. GRACE Data

GRACE can monitor TWSA (terrestrial water storage anomaly) through the inversion of changes in the Earth’s gravity field. TWSA is the equivalent height of water (mm) relative to the long-term baseline average (January 2004–December 2009), which is available from the GRACE RL06 mascon from the University of Texas Space Research Center (http://www2.csr.utexas.edu/grace/) (accessed on 20 August 2022) for the period from April 2002 to June 2017, with a spatial resolution of 0.25° × 0.25° and a time resolution of 1 month. The mascon solutions have gone through C20 replacement, degree 1 correction and GIA (glacial isostatic adjustment) correction. Compared with spherical harmonies (SH) method, mascon-derived TWS have improved resolution and signal-to-noise ratio, and reduced leakage error [41]. The linear interpolation method was used to compensate for missing values due to technical issues [42].

2.2.2. GLDAS Data

The National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA) jointly developed the Global Land Data Assimilation System (GLDAS), which combines satellites and ground observations data [43]. The capabilities of Noah model driven by GLDAS for obtaining groundwater has been successfully verified, and Noah model has high spatial resolution (0.25° × 0.25°) consistent with GRACE RL06 CSR mascon solutions [23]. Thus, this paper adopted GLDAS Noah model L4 monthly V2.1 products (http://disc.sci.gsfc.nasa.gov/hydrology/data-holdings) (accessed on 25 August 2022) to estimated soil moisture at four layer depths of 10, 40, 100, and 200 cm, canopy water storage, total runoff including surface runoff, subsurface runoff, and snowmelt runoff, and snow water equivalent [44].

2.2.3. ERA5-Land Data

ERA5 is the fifth-generation reanalysis from ECMWF (European Centre for Medium-Range Weather Forecasts), with an updated physics model and improved data assimilation process compared to ERA-Interim data [45,46]. ERA5-Land is a reanalysis dataset providing a consistent view of the evolution of land variables over several decades at an enhanced resolution, which is produced by replaying the land component of the ECMWF ERA5 climate reanalysis [47]. This study used total precipitation (accumulated liquid and frozen water, including rain and snow), total evapotranspiration (accumulated amount of water evaporated from the Earth’s surface), and runoff (the sum of surface and sub-surface runoff) from “ERA5-Land monthly averaged data from 1950 to present” as meteorological variables to investigate the attribution of natural factors to drought (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=form) (accessed on 25 August 2022).

2.2.4. NDVI

A change in vegetation can indicate the occurrence of drought [48]. NDVI (Normalized Difference Vegetation Index) is the normalized reflectance difference between the near-infrared (NIR) and visible-red bands. Higher NDVI values generally represent greater vigor and photosynthetic capacity (or greenness) of vegetation canopy [49]. Monthly NDVI dataset with 1 km resolution of China from 2001 to 2022 was used in this study (available at National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://www.geodata.cn) (accessed on 20 June 2023)). The dataset was generated through the monthly synthesis, mosaicking, and cropping from MODIS MOD13A2 data.

2.2.5. Teleconnection Factor

According to previous studies [50,51,52,53] ENSO, PDO, NAO, and AO are closely related to drought in different regions of China in recent years. This study used Multivariate El Niño/Southern Oscillation Index version 2 (MEI.v2) provided by NOAA Earth System Research Laboratory (https://psl.noaa.gov/enso/mei/) (accessed on 20 June 2023), Pacific Decadal Oscillation, North Atlantic Oscillation, and Arctic Oscillation supplied by NOAA National Centers for Environmental Information (http://www.ncdc.noaa.gov/ teleconnections/) (accessed on 20 June 2023) to study the effects of teleconnection factors on drought.

2.2.6. Land Use and Land Cover

The Landsat-derived annual land cover product of China (CLCD) from 1985 to 2019 built by Yang and Huang [54] was utilized in this study. CLCD reflects China’s rapid urbanization and a series of ecological projects, revealing the anthropogenic impact on land cover under climate change conditions (https://zenodo.org/record/5816591) (accessed on 20 June 2023). The CLCD product has 9 types of land use and land cover: cropland, forest, shrub, grassland, water, snow, barren, urban and wetland.

2.3. Methods

2.3.1. Drought Index and Drought Characteristic

GRACE-DSI (drought severity index) proposed by Zhao et al. [12] was obtained using the following algorithm, which was based solely on GRACE TWSA and is therefore a hydrological drought index.

$${DSI}_{i,j}=\frac{{TWSA}_{i,j}-{\overline{T}WSA}_j}{std\left({TWSA}_j\right)} .$$

(1)

where $i$ represents the $i$-th year from 2002 to 2017, $j$ represents the $j$-th month from 1 to 12, ${\overline{T}WSA}_j$ and $std\left({TWSA}_j\right)$ represent the average and deviations of $TWSA$ for month $j$, respectively. DSI was classified into five drought (wet) categories by the United States Drought Monitor (USDM) (

Table 1. Classification of dry/wet conditions of drought severity index.

Category	Description	DSI
W5	Exceptionally wet	[2.00, $+\mathrm{\infty }$)
W4	Extremely wet	[1.60, 2.00)
W3	Severely wet	[1.30, 1.60)
W2	Moderately wet	[0.80, 1.30)
W1	Slightly wet	[0.50, 0.80)
N0	Near normal	(−0.50, 0.50)
D1	Abnormally drought	(−0.80, −0.50]
D2	Moderate drought	(−1.30, −0.80]
D3	Severe drought	(−1.60, −1.30]
D4	Extreme drought	(−2.00, −1.60]
D5	Exceptional drought	(−$\mathrm{\infty }$, −2.00]

). Drought characteristics were extracted using run theory (duration, severity, etc.) and pixel-based drought evaluation (affected area) [11].

Drought event was defined as a period with DSI values less than the abnormal drought threshold (−0.50) over three or more consecutive months, and two adjacent drought events separated by a month were merged into a single drought event. Based on run theory, drought duration from onset to recovery was set as drought duration. Drought affected area was expressed as the ratio of drought pixels to total pixels [55].

2.3.2. Retrieval of Groundwater Storage Anomalies

As shown in Equation (2), groundwater storage anomalies (GWSA) can be estimated by deducting the contribution of soil moisture storage anomalies (SMSA), snow water equivalent anomalies (SWEA), canopy water storage anomalies (CWSA), and surface water storage anomalies (SWSA) provided by the GLDAS Noah model from GRACE TWSA, where the contribution of biological water reserves are negligible [23]. For consistency with the GRACE data, these water storages were deducted from the mean values for the period January 2004–December 2009 to obtain the corresponding monthly anomalies.

$$GWSA=TWSA-SMSA-SWEA-CWSA-SWSA .$$

(2)

2.3.3. Theil–Sen Trend Analysis and Mann–Kendall Trend Test

Theil–Sen trend analysis is a robust non-parametric statistical trend calculation method for estimating the slope of a linear trend in a time series, without requiring the sample to follow a certain distribution, and it is not be disturbed by outliers [56,57].

$$\beta =median\left(\frac{x_j-x_i}{j-i}\right), \hspace{1em}\forall j>i .$$

(3)

where $x_i$ and $x_j$ are the time series monitoring data in sequential order. If $\beta$> 0 reflects an upward trend, $\beta$< 0 reflects a downward trend.

Mann–Kendall trend test (M-K trend test) is a non-parametric statistical method, which is used to determine the significance of trends [56,57]. The calculations are performed as follows:

$$S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}sign(x_j-x_i) ,$$

(4)

$$sign(x_j-x_i)=\left\{\begin{array}{c}1,\hspace{1em} x_j-x_i>0\\ 0,\hspace{1em} x_j-x_i=0\\ -1,\hspace{1em} x_j-x_i<0\end{array}\right. ,$$

(5)

$$var(S)=\frac{n(n-1)(2n+5)-{\sum }_{k=1}^m{t}_k(t_k-1)(2t_k+5)}{18} ,$$

(6)

$$Z_c\left\{\begin{array}{c}\frac{S-1}{\sqrt{var\left(S\right)}},\hspace{1em}\hspace{1em}S>0\\ 0,\hspace{1em} S=0\\ \frac{S+1}{\sqrt{var}},\hspace{1em}\hspace{1em}S<0\end{array}\right. .$$

(7)

where $n$ is the length of the time series, $t_k$ is the number of the tie in the sample in the kth value, and m is the number of tied variables. At a given significance level $\mathsf{\alpha }$, a value of|$Z_c$| > $Z_{1-\alpha /2}$ indicates that the trend has passed the significance test. The significant level $\mathsf{\alpha }$ is usually 0.01 and 0.05, and the corresponding $Z_{1-\alpha /2}$ is 2.5758 and 1.9600, respectively.

2.3.4. Random Forest (RF)

Random Forest (RF), an ensemble classifier method based on decision trees, can be used to solve classification and regression problems [58]. Random Forest (RF) is capable of generating multiple independent decision trees simultaneously and randomly selecting the input training sample variables, as well as choosing the variables at each tree node with good accuracy and generalization ability [59]. The number of decision trees and the nodes that split the variables are two important parameters of the Random Forest, and the optimal result will be obtained when the optimal number of trees and leaf nodes are determined [60].

2.3.5. Cross Wavelet Transforms

The cross wavelet transforms, a method for analyzing the correlation of two time series, is often used to reveal the correlation between hydrological series and climatic series. This method reflects regions where two time series periodically have high common power between the time and frequency domains [61]. In the result of cross wavelet transforms, the arrows indicate the phase difference, with left (right) indicating that the two sequences are negatively (positively) correlated; and up (down) indicating that the change in former lags (exceeds) the change in the latter by 3 months.

Define the continue wavelet transform of two time series $X=\left\{x_1,x_2,\dots ,x_n\right\}$ and $Y=\left\{y_1,y_2,\dots ,y_n\right\}$ is $W_n^X\left(S\right)$ and $W_{\mathrm{n}}^{Y*}\left(S\right)$, and the cross wavelet transforms between them is

$$W_n^{XY}\left(S\right)={W_n^X\left(S\right)W}_n^{Y*}\left(S\right) .$$

(8)

where $W_{\mathrm{n}}^{Y*}\left(S\right)$ denotes the complex conjugation of $W_{\mathrm{n}}^Y\left(S\right)$, $S$ is time delay, and $W_n^{XY}\left(S\right)$ is cross wavelet power.

The theoretical distribution of the cross wavelet power of two time series with background power spectra $P_k^X$ and $P_k^Y$ is described as follows [62]:

$$D(\frac{\left|W_n^X\left(S\right)W_{\mathrm{n}}^{Y*}\left(S\right)\right|}{\sigma X\sigma Y}<p)=\frac{Z_v\left(p\right)}{v}\sqrt{P_k^X{P}_k^Y} .$$

(9)

where $\sigma X$ and $\sigma Y$ represent the standard deviation of $X$ and $Y$, respectively. $Z_v\left(p\right)$ is the confidence level associated with the probability $p$ for a pdf defined by the square root of the product of two ${\chi }^2$ distributions, $v$ is the degree of freedom, $v$ is equal to 1 for real ($Z_1=2.182$, with a significance level of 0.05) and 2 ($Z_2=3.999$, with a significance level of 0.05) for complex wavelets.

3. Results

3.1. Drought Characteristics and Trends

Based on the Drought Severity Index (DSI) derived from the GRACE data, the drought duration and affected area calculated using the run theory and pixel-based drought evaluation approach are revealed in Figure 2. The results showed that almost all humid and arid basins experienced major drought periods of 2002–2006 and 2014–2017, respectively. The heaviest drought with the longest drought duration and highest drought severity identified using the DSI occurred in YZRB, SERB, and PRB (humid basins) in 2002/9–2003/5, 2003/5–2004/7, and 2004/10–2006/5, respectively. Similarly, LRB, HLRB, HRB, YRB, and IRB (arid basins) experienced the most serious drought period during 2014/7–2016/9, 2014/7–2016/9, 2013/8–2014/9, 2015/8–2017/6, and 2014/9–2017/6, respectively. In addition, the arid basin SRB and humid basin SWRB were detected as experiencing long-term drought from 2007/7 to 2008/12 and from 2014/4 to 2017/6, respectively.

Figure 3 demonstrated the spatial characteristics of hydrological drought trends derived using the DSI over ten basins in China via the Theil–Sen trend analysis and M–K trend test. There was obvious spatial heterogeneity in the variation of DSI over the study period. Significant decrease trends were observed in HLRB, the northern and eastern parts of IRB, central and western parts of LRB, northern part of YRB, northern part of HRB, and northern part of SWRB, with a slope ranging from −0.0173 to −0.0036. Moreover, southern IRB (Tibetan Plateau) and central YZRB (Szechuan Basin and east Yunnan and Guizhou Plateau) revealed significant increase trends, with a slope ranging from 0.0033 to 0.0171. It is important to investigate the dominant divers of drought spatiotemporal heterogeneity of China’s major river basins.

3.2. Attribution Analysis of Drought

To comprehensively assess the dominant factors influencing hydrological drought in China, the effect of natural factors (including terrestrial water storage components, precipitation, evapotranspiration, runoff, NDVI, and teleconnection factors (ENSO, PDO, NAO and AO)) and anthropogenic influences (LULC) on drought were identified and quantified in terms of water storage components, climate change, and water balance.

3.2.1. From the Perspective of Water Reserve Component

To identify the main water storage components causing drought, the spatial distribution (Figure 4) and temporal evolution (Figure 5) of the trends of various components were investigated using an M–K trend test and Theil–Sen trend analysis. To clearly see the trend of CWSA, the results of CWSA in Figure 5 were enlarged by a factor of 300. As shown in Figure 4, the GWSA had a significant downward trend in most of the LRB, HLRB, HRB, YRB, SWRB, and IRB, with a trend of <−4–0 mm/month, especially in central SWRB (Southeastern Tibet), southern HLRB (North China Plain), and northwestern IRB (Tianshan Mountains), where it was highly significant. On the contrary, the GWSA in central SRB, southwestern IRB, and most of YZRB showed a prominent upward trend of 0->1.5 mm/month, especially with respect to the remarkable increasing trend in the southwestern part of IRB (Qaidam basin). It was obvious that the trends of GWSA were consistent with those of TWSA. Moreover, as shown in Figure 5, it should be noted that the changes in GWSA for SERB and PRB were not significant, while the prominent upward trends in SWSA (0.18 mm/month) and SMSA (0.133 mm/month) explained the significant increase in TWSA in SERB and the middle of PRB; the significant decrease in CWSA (−1.76 × 10⁻⁴ mm/month) in the central part of YZRB neutralized the significant increase in GWSA (0.300 mm/month) and slowed down the rate of increase in TWSA; SWEA showed a significant downward trend only in the northern SWRB (−0.0525 mm/month) and southwestern IRB (−8.65 × 10⁻⁵ mm/month), contributing to the significant decline in TWSA; and the SMSA showed a clear downward trend in most parts of SWRB (−0.326 mm/month) and HRB (−0.263 mm/month), and contributed to the decline in TWSA. In addition, the rise of SMSA in most parts of IRB (0.0198 mm/month) and SRB (0.221 mm/month) contributed to the rise and decline in TWSA, respectively. The specific results are summarized in

Table 2. Trends of TWSA and TWSA components of China’s major river basins (mm/month) (where * represents regions where water storage is statistically significant ($\mathsf{\alpha }$ = 0.05, |$Z_c$| > 1.96). Red and blue font color for basins’ name highlight arid and humid basins respectively).

Basins	TWSA	GWSA	SWSA	SWEA	CWSA	SMSA
SRB	0.085	−0.162 *	0.0060	8.13 × 10−5	6.43 × 10−5	0.221 *
LRB	−0.453 *	−0.474 *	0.0042	2.72 × 10−6	4.43 × 10−5	0.0147
HLRB	−1.241 *	−1.150 *	−0.0006	0	−7.19 × 10−6	−0.004
HRB	−0.486 *	−0.170 *	−0.0142	0	−9.75 × 10−5	−0.263 *
YRB	−0.513 *	−0.457 *	−0.0015	4.10 × 10−6	−3.24 × 10−5	0.003
YZRB	0.264 *	0.300 *	0.0020	−4.25 × 10−4 *	−1.76 × 10−4 *	−0.029
SERB	0.441 *	0.0572	0.1800 *	0	−7.17 × 10−5	0.133 *
PRB	0.414 *	0.170	0.0717	0	−1.37 × 10−4	0.133
SWRB	−0.867 *	−0.431 *	−0.0230	−0.0525 *	−1.28 × 10−4	−0.326 *
IRB	−0.166 *	−0.184 *	−0.0028	−8.65 × 10−5	4.05 × 10−6	0.0198 *

.

3.2.2. From the Perspective of Climate Change

In order to analyze the driving forces on hydrological drought of China’s major river basins from the perspective of climate change, this section used cross wavelet transforms to investigate the relationship between teleconnection factors (ENSO, PDO, NAO, and AO) and the DSI in the time and frequency domains. The results of the YRB (typical arid basin) and the YZRB (typical humid basin) are presented below as examples.

For the YRB, as shown in Figure 6a, the DSI exhibited unstable intermittent resonance cycles with ENSO, and a similar situation occurred for the NAO (Figure 6e). The correlations between DSI and PDO passed the significance test from 2008 to 2015 with a 14–32 month-long resonance cycle (Figure 6c). There was a prominent and negative correlation between DSI and AO from 2007 to 2011, with a 24–32 month-long signal, and a positive correlation between DSI and AO from 2012 to 2014, with a 15–20 month-long signal (Figure 6g). For the YZRB, as shown in Figure 6b, the DSI had a significant negative correlation with ENSO, and they had a 14–28 month-long resonance cycle from 2005 to 2011. The DSI exhibited a significant lag of 3 months with PDO, and there was a resonance cycle of 8–12 months from 2011 to 2013 (Figure 6d). The NAO showed a significant negative correlation with DSI for 24–32 months from 2008 to 2014 (Figure 6f). There were unstable intermittent resonance cycles between DSI and AO (Figure 6h).

By comparing the strength of significant resonance cycles, it can be concluded that there was a statistically significant correlation between DSI and teleconnection factors, and the teleconnection factors had a great influence on the evolution of hydrological drought derived from the DSI. Furthermore, the impacts of the PDO and ENSO were stronger than other factors in the YRB and YZRB, respectively. Other results can be found in the Supplementary Material Figure S1(1)–(4). Through a similar analysis,

Table 3. The most important teleconnection factors in different basins (red and blue font colors for basins’ names highlight arid and humid basins, respectively).

Basins	SRB	LRB	HLRB	HRB	YRB	YZRB	SERB	PRB	SWRB	IRB
Factors	AO	NAO	NAO	PDO	PDO	ENSO	ENSO	PDO, NAO, AO	ENSO	AO
Cycles	35–64 months in 2007–2012	20–35 months in 2007–2014	24–32 months in 2008–2013	20–32 months in 2009–2014	14–32 months in 2008–2015	14–28 months in 2005–2011	12–48 months in 2008–2014	12–16 months in 2005–2009	20–64 months in 2006–2014	25–45 months in 2008–2014

summarizes the most important teleconnection factors affecting drought and the main significant resonance cycles in all basins. It can be concluded that most of the significant resonance cycles were 12–64 month-long signals in the period of 2005–2014, and ENSO was the strongest teleconnection factor in most humid basins, while NAO, PDO, and AO were the strongest teleconnection factors in arid basins and the PRB.

3.2.3. From the Perspective of Water Balance

According to the water balance, TWSC was caused by natural factors of precipitation, evaporation, runoff, NDVI, teleconnection factors, and anthropogenic factors such as LULC change. This section used a Pearson correlation analysis and Random Forest (RF) to attribute hydrological droughts characterized by insufficient TWS in China’s major river basins. The following section takes YRB and YZRB as examples to select the dominant drought factors.

The variable importance scores obtained using the RF were used to further quantify the contribution of different influencing factors on hydrological drought. As shown in Figure 7, the variable importance scores of LULC were equal to or even higher than that of the natural factors for most basins. For instance, the importance of cropland (0.49), forest (0.36), shrub (0.55), grassland (0.35), and water (0.40) was high in the YZRB (Figure 7f), and had a negative correlation (−0.54–0.1) with the DSI (Figure S2f); and similarly, the importance of urban (0.58) was high, but positively correlated (0.52) with the DSI (Figure S2f). Due to the dramatic increase in urban coverage and the reduction in cropland and grassland [37], the YZRB tended to become wetter (Figure 3), which suggested that the negative impact of urban expansion is compensated by cropland and grassland reduction. Runoff and ENSO had relatively higher variable importance scores (0.45, 0.45) and higher positive correlation coefficients (0.38, 0.32) than other natural factors. In the YRB, forest (0.78), grassland (0.45), urban (0.78) and cropland (0.43) received relatively high variable importance scores (Figure 7e), and the area of forest, grassland and urban land gradually increased while cropland decreased [37]. Figure 3 revealed that most parts of the YRB had a drying trend, which was negatively related to forest (−0.53), urban (−0.68), and grassland (−0.12) cover, and had a positive relationship with cropland (0.30) (Figure S2e). Therefore, the expansion in forest, urban, and grassland cover induced drought and the decrease in cropland failed to reverse the drying trend in the YRB. Runoff and PDO had relatively higher variable importance scores (0.40, 0.30) and higher positive correlation coefficients (0.23, 0.14) than other natural factors.

It was found that in many basins (Figure 7 and Figure S2), LULC had an almost equal importance with natural factors in the SRB, HRB, YRB, YZRB, SERB, and PRB, and an even higher importance than natural factors in the LRB, HLRB, SWRB, and IRB. Moreover, the influence of teleconnection factors could not be ignored in all basins, and the influence of runoff had a relatively high importance in the SRB, HRB, YRB, YZRB, SERB, and PRB. In summary, the dominant natural or human factors affecting drought formation in other river basins could be derived from the attribution and correlation of meteorological variables, teleconnection factors, and LULC with the DSI.

4. Discussion

TWSA is a crucial indicator of regional water resources changes triggered by the natural water cycle and human exploitation [63]. In this study, a Drought Severity Index (DSI) based on GRACE TWS data was employed for monitoring hydrological droughts and investigating spatiotemporal changes in drought. Additionally, the study analyzed dominant factors contributing to hydrological droughts by considering water storage components, climate change, and water balance.

Changes in terrestrial water storage are a comprehensive reflection of changes in snow and ice, runoff, canopy water storage, soil water, and groundwater. Hydrologic drought caused by insufficient TWS is influenced by water storage components contributing to the decreasing trend in TWSA. These components will exacerbate the drought. Conversely, the water storage components that contribute to the increasing trend of TWSA will mitigate or suppress the drought. The results of the trend analysis revealed a spatially and temporally significant consistency between the groundwater storage anomaly (GWSA) and terrestrial water storage anomaly (TWSA) in most of the basins, indicating both upward and downward trends. For instance, in arid basins, groundwater storage is the main source of water supply compared with other components. For example, in the southern HLRB (North China Plain), intense human activities (agricultural irrigation, industrial and domestic water demand) affect terrestrial water reserves mainly through the exploitation of groundwater [23]. SWSA and SMSA were attributed to the TWSA increase in SERB and PRB, and CWSA retarded the increase in the TWSA trend in the YZRB. SWEA caused a significant decrease in the TWSA of the SWRB and IRB, and SMSA caused a significant decrease in TWSA of the SWRB and HRB. Hence, a GWSA decrease was the dominant drought factor in the perspective of water storage components.

From the perspective of climate change, a statistically significant correlation between DSI and teleconnection factors could be derived from the cross wavelet transforms. The influences of PDO and ENSO were greater than those of other factors in the YRB and YZRB, respectively. From the Supplementary Material Figure S1, it could be concluded that the ENSO is the strongest teleconnection factor in most humid basins, and the NAO, PDO, and AO are the strongest teleconnection factors in arid basins and the PRB. Most of the significant resonance cycles were 12–64 month-long signals in the period of 2005–2014, which indicated that teleconnection factors have a great influence on hydrological drought.

As a result of the attribution analysis based on RF from the perspective of water balance, LULC had an equal or greater impact on hydrological drought than other drivers in most basins. LULC change will alter the water cycle and result in TWS deficits. The expansion of urban cover can create a higher terrestrial water demand for urban production and livelihood, which poses a risk of drought [64]. Moreover, the increase in cropland, forest, and grassland will lead to the increase of water demand and aggravate drought. The implementation of the “Grain-to-Green Program” on the Loess Plateau in northern China has increased water consumption, resulting in a runoff reduction and soil moisture decrease, and apparent widespread aridification [65,66,67]. In addition, some studies have found that irrigation return water can recharge groundwater on cropland and mitigate drought [28]. Lv et al. [38] illustrated that the increasing trend of evapotranspiration resulting from the expansion of vegetation cover and irrigation water played a dominant role in the decline of total water storage capacity (TWSC) in the Loess Plateau, considering the “Grain-to-Green Program” and irrigation. In summary, TWSC, as well as drought variability due to LULC changes, is complex, especially when anthropogenic factors except for LULC changes are included. These complex factors interact and balance each other, ultimately leading to drought or wet trends in a basin.

In addition, except for the IRB and SWRB, runoff of natural aspects had relatively high variable importance scores in most basins, especially in the SRB, HRB, YRB, YZRB, SERB, and PRB. This result was consistent with Yang et al. [68], who proposed that runoff made a prominent contribution to the TWS of the entire territory, except for the northwest river basins. This suggests that runoff (main output fluxes) plays a more important role than precipitation (main input fluxes) and evapotranspiration (the other main input fluxes) in drought derived from TWS deficits in most basins. Furthermore, teleconnection factors have a strong influence on drought variations. In summary, the results revealed the enormous impact of human activities on hydrological droughts in China’s major basins, and may inspire projects related to agricultural and urban development.

5. Conclusions

Climate change and human activities are the two main factors influencing changes in TWS [67]. In this study, the variability and attribution of drought derived from TWS deficits were analyzed from 2002–2017 in China. The main results are as follows:

• Almost all humid and arid basins experienced major drought periods during 2002–2006 and 2014–2017, respectively. The southern IRB and central YZRB had notable declines in DSI trends, while most parts of the HLRB, IRB, LRB, YRB, HRB, and SWRB experienced significant increases in DSI trends;
• Abnormal groundwater decreases were the main cause of drought triggered by TWS deficits in most basins;
• ENSO is the strongest teleconnection factor in most humid basins, and the NAO, PDO, and AO are the strongest teleconnection factors in arid basins and the PRB. Most of the significant resonance cycles were 12–64 month-long signals in the period of 2005–2014;
• The impact of human activities (LULC) has become equally or even more significant than natural factors such as runoff and teleconnection factors in influencing hydrological drought in most basins of China.

This study confirms the dominant factors from comprehensive perspective on hydrological drought in China’s major river basins, which has great significance for optimizing the allocation of regional water resources and inspiring development for agriculture and urbanization in China.