Spatiotemporal Variation Characteristics of Droughts and Their Connection to Climate Variability and Human Activity in the Pearl River Basin, South China

Abstract

Droughts have damaging impacts on human society and ecological environments. Therefore, studying the impacts of climate variability and human activity on droughts has very important scientific value and social significance in order to understand drought warnings and weaken the adverse impacts of droughts. In this study, we used a combined drought index based on five Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On solutions to characterize droughts in the Pearl River basin (PRB) and its sub-basins during 2003 and 2020. Then, we accurately quantified the impact of climate variability and human activity on droughts in the PRB and seven sub-basins by combining the hydrometeorological climate index and in situ human activity data. The results show that 14 droughts were identified in the PRB, particularly the North River basin with the most drought months (52.78%). The El Niño-Southern Oscillation and the Indian Ocean Dipole were found to have important impacts on droughts in the PRB. They affect the operation of the atmospheric circulation, as well as the East Asia summer monsoon, resulting in a decrease in precipitation in the PRB. This impact shows a significant east–west difference on the spatial scale. The middle and upper reaches of the PRB were found to be dominated by SM, while the lower reaches were found to be dominated by GW. Human activity was found to mainly exacerbate droughts in the PRB, but also plays a significant role in reducing peak magnitude. The sub-basins with a higher proportion of total water consumption experienced more droughts (more than 11), and vice versa. The Pearl River Delta showed the highest drought intensification. Reservoir storage significantly reduces the drought peak and severity, but the impact effect depends on its application and balance with the total water consumption. Our study provides a reference for analyzing the drought characteristics, causes, and impacts of sub-basins on a global scale.

1. Introduction

A drought is a very destructive hydrometeorological disaster, which is an extreme manifestation of the paradox between water supply and demand [1]. It can seriously threaten the sustainable development of the social economy [2,3]. Since 2000, the number and duration of global droughts have increased by 29%. Between 1998 and 2017, global economic losses due to drought were high as USD 124 billion, and the number of people affected by drought reached 1.4 billion [4]. Therefore, realizing the early warning signs of drought is of great significance in order to reduce the adverse effects of droughts. The premise of establishing a drought-based early warning model is to fully understand and master the formation and development mechanism of droughts, especially the impact of various factors on droughts.

Drought research requires various hydrological data, such as data on precipitation (PPT), soil moisture (SM), runoff, evapotranspiration (ET), etc. Although the ground stations can directly obtain these hydrological data, the insufficient site leads to low data coverage, and the method needs high costs [5]. Satellite remote sensing achieves all-weather tracking and full coverage of data monitoring, but it obtains hydrological data indirectly. Therefore, there are large uncertainties in the data [6,7]. The emergence of Gravity Recovery and Climate Experiment (GRACE) and its Follow-On (GRACE-FO) satellites has provided us with a geodetic approach that can be used to monitor all the terrestrial hydrological data, including surface and subsurface data, and terrestrial water storage changes (TWSCs) estimated by these data included the influence of natural factors and human activity [8,9]. Therefore, it has attracted significant attention among scientists from all over the world. Moreover, some scholars have used them in terrestrial hydrological research, e.g., for the Yangtze River basin [10,11] and the Amazon River basin [12,13], as well as in Africa [14,15] and Australia [16,17].

A drought index is a quantitative indicator often used in drought research, and it has the advantage of being easy to understand. Commonly used drought indices mainly include the Palmer drought severity index (PDSI), the standardized PPT ET index (SPEI), and the standardized PPT index (SPI) [18]. However, the above three drought indices are only based on one or two pieces of hydrometeorological data, so it cannot fully reflect the real situation of terrestrial hydrology [7,19,20]. The current causes of droughts are not only natural factors, but also human activity. The traditional drought indices cannot accurately reflect the impact of human activity, bringing great uncertainty within drought assessments [21]. Moreover, GRACE has provided an effective solution to the above problem. Sinha et al. [22] used GRACE TWSC data to construct a water storage deficit index (WSDI), and applied it to drought monitoring in the Indian Continent during 2002 and 2015. Subsequently, the WSDI was used in drought research in the Yangtze River basin and the Sao Francisco River basin, as well as China, and Turkey [7,23,24,25]. The WSDI has become a preferable and effective tool for monitoring large-scale droughts.

The Pearl River basin (PRB) is the largest river system in South China and the third longest river in China. Its annual runoff is more than 330 million m³, ranking second in China’s river systems [26]. GRACE data were applied to successfully detect severe drought events in the PRB in 2011 [27], and PPT, ET, runoff, and climate variability were found to play an important role in the drought propagation [28,29]. The El Niño-Southern Oscillation (ENSO) alters the drought duration and severity by affecting PPT levels [30]. Except for natural factors, human activity has a certain interference effect on droughts. Forest change has the ability to disturb regional climate states and thus contribute to hydrological droughts. In Baker Creek and Yaak River basins, forest change has a greater impact on droughts than climate change [31,32]. Reservoir storage also disrupts drought duration and severity [30]. The above studies did not consider the divergence between different sub-basins in the PRB, and precisely quantify the contribution of human activity to droughts. In this study, we discuss the impact of climate variability and human activity on each PRB sub-basin during the drought events from 2003 to 2021, and precisely quantified the impact extent of human activity to drought duration and severity. We briefly introduced the study area, data, and methods in Section 2, Section 3 and Section 4, respectively. Section 5 presents an analysis of drought characteristics and its influencing factors in the PRB. The discussion and conclusion are provided in Section 6 and Section 7, respectively.

2. Study Area

The PRB (Figure 1), approximately located at 21° N–27° N and 102° E–116° E, covers five provinces in South China and northeastern Vietnam, with a total area of 442,000 km². The terrain is high in the northwest and low in the southeast. The region has the subtropical monsoon climate, with a mild and rainy climate. The annual average temperature is 14~22 °C and the annual average PPT is 1200~2200 mm. PPT distribution obviously decreases gradually from east to west. The spatial and temporal distribution of PPT is uneven, and more than 80% of PPT is concentrated in the rainy season [33,34]. The sub-basins of the PRB are the Nan-Bei Pan River basin (NBPRB), the Hongliu River basin (HRB), the Yu River basin (YRB), the West River basin (WRB), the North River basin (NRB), the East River basin (ERB), and the Pearl River Delta (PRD).

3. Data

3.1. GRACE/GRACE-FO Data

At present, GRACE/GRACE-FO solutions are mainly published by three official agencies, which are the Center for Space Research at the University of Texas at Austin (CSR), the Helmholtz-Centre Potsdam-German Research Centre for Geosciences (GFZ), and the Jet Propulsion Laboratory (JPL). GRACE/GRACE-FO solutions are divided into spherical harmonic (SH) and Mascon solutions, and gridded TWSC data can be obtained from them. Among them, only SH solutions require pretreatment. Firstly, first-order and C₂₀ term coefficients needed to be replaced [35,36]. Then, SH coefficients were processed using a 250 km fan filter and a P3M8 polynomial filter. Finally, the scale factor method was used to repair the signal [37]. Since GRACE and GRACE-FO data are exactly the same, GRACE and GRACE-FO are collectively referred to as GRACE. For convenience, three GRACE SH solutions and two Mascon solutions are termed as CSR-SH, GFZ-SH, JPL-SH, CSR-M, and JPL-M, respectively.

3.2. Reconstructed TWSC Data

The data gap between GRACE and GRACE-FO missions (from July 2017 to May 2018) was filled using the dataset of reconstructed TWSC data in China based on PPT (2002–2019), derived from the National Tibetan Plateau Data Center. The dataset was calculated using the CSR GRACE/GRACE-FO RL06 Mascon solution, China’s daily gridded PPT real-time analysis system (version 1.0), and CN05.1 temperature data and other datasets, according to the PPT reconstruction model [38,39].

3.3. Traditional Drought Index

We used the two traditional drought indices, the self-calibrating PSDI (SCPDSI) and the SPEI, to evaluate the quality of WSDIs. The SCPDSI is an upgraded version of the PDSI, which is more suitable for drought assessment on a global scale and calculated based on PPT, temperature, and parameters relating to soil/surface characteristics [40]. The monthly 0.5° × 0.5° global gridded SCPDSI data were taken from the Climate Research Unit at University of East Anglia. The SPEI was estimated based on PPT and potential ET. In our study, the monthly SPEI gridded data with a spatial solution 0.25° × 0.25° had three different time scales (3, 6, and 12 months). The 3-month scale represents an agricultural drought; the 6-month scale represents a hydrological drought; and the 12-month scale represents a hydrogeological drought [41].

3.4. In Situ Hydrometeorological Data

In situ hydrometeorological data, including PPT, soil moisture (SM), runoff, and ET, were derived from China Global Land Surface Reanalysis Products (CMA-RA/Land). The products were provided by the China Meteorological Administration, with a spatial solution of 0.5° × 0.5° and monthly temporal solutions.

3.5. Groundwater Data

Since groundwater (GW) also has a certain impact on regional drought, we discussed GW as a factor of drought. In our study, the monthly 0.25° × 0.25° GW gridded data were taken from the Global Land Data Assimilation System 2.1 Catchment Land Surface Model, which is obtained by subtracting soil moisture, plant canopy surface water, and snow water equivalent from terrestrial water storage.

3.6. Climate Teleconnection

Droughts have a certain teleconnection with global climate events [42]. Climate change in China is vulnerable to sea surface temperature anomalies in the Pacific and Indian Oceans [43]. The ENSO is a persistent anomalous sea surface temperature in the equatorial Pacific that affects global and regional PPT levels and water cycles, especially in Pacific coast regions [44,45]. Previous studies [30,46] indicate that PPT and TWSC in the PRB are mainly influenced by the ENSO. PPT in the PRB is closely related to the Indian Ocean Dipole (IOD) [47,48]. Therefore, it is necessary to study the impact of the ENSO and the IOD on climate change in the PRB. In our study, the monthly Niño 3.4 index and the Indian Ocean Dipole model index (DMI) during 2003 and 2020, respectively, were provided by the National Oceanic and Atmospheric Administration, which has represents the state of the ENSO and the IOD over time.

3.7. In Situ Human-Induced TWSC Data

In situ human-induced TWSC data were provided by the PRB Water Resources Bulletins, and the bulletins were published by the Pearl River Water Resources Commission of the Ministry of Water Resources. These data comprise PPT, surface water, groundwater, reservoir storage (RS), total water supply, total water consumption (TWC), etc., both in the PRB and its sub-basins on an annual scale.

4. Methods

4.1. Data Fusion

Data fusion is mainly carried out in three steps: (1) the uncertainties of five TWSC results of different GRACE solutions are estimated using the generalized three-cornered hat method; (2) according to the uncertainties results, the weight of the five TWSC results are determined; (3) the five TWSC results are fused based on the least squares method. Technical details on the data fusion can be found in Refs. [49,50].

4.2. Natural-Induced TWSCs

To analyze the impact of human factors on the droughts in the PRB, we need to evaluate the natural-induced TWSC ($TWS{C}_c$). According to the water balance equation, $TWS{C}_c$ is expressed as follows:

$$TWS{C}_c=P-R-E$$

(1)

where $P$ is the PPT rate, $E$ is the ET rate, and $R$ is the runoff.

4.3. WSDI

The expression of the WSDI is as follows [22]:

$$WS{D}_{i,j}=TWS{C}_{i,j}-TWS{C}_j^{c\lim }$$

(2)

where $WS{D}_{i,j}$ and $TWS{C}_{i,j}$ are the water storage deficit and GRACE TWSC for the $j$th month in the year $i$, respectively; $TWS{C}_j^{c\lim }$ is the climatology of TWSC for the $j$th month. In our study, the climatology of TWSC was calculated by averaging each month of GRACE TWSC data during 2003 and 2020 in the PRB, which was used as the baseline for identifying the WSD [51]. A negative WSDI represents the drought conditions, whose magnitude indicates the drought intensity [22].

4.4. Definition and Characteristics of Drought Event

In our study, to determine the occurrence of a drought event in the PRB, the monthly WSDI must be less than 0 for three consecutive months. The drought duration is the number of months during the start and end month of the drought event; the peak magnitude indicates the maximum monthly WSDI during the drought event; the drought area percentage represents the ratio of the grid area with a WSDI less than 0 to the total area of the PRB. The drought severity (S) is expressed as [52]:

$$S=\overline{M}\times D$$

(3)

where $\overline{M}$ is the monthly average of the WSDI during the drought event and $D$ is the drought duration calculated up to the current month.

4.5. SPI

The SPI is estimated based on long-term PPT data and is widely used in regional drought research [53]. Firstly, the probability density function of gamma distribution is used for long-term PPT data. Then, the function is converted to a normal distribution to obtain an SPI mean value at 0. The expression is as follows:

$$SPI=W-\frac{c_0+c_{1}W+c_2{W}^2}{1-d_{1}W+d_2{W}^2+d_3{W}^3}$$

(4)

$$W=\left\{\begin{array}{l}\sqrt{-2\ln P},P\le 0.5\hfill \\ \sqrt{-2\ln \left(1-P\right)},P>0.5\hfill \end{array}\right.$$

(5)

where $P$ is the cumulative probability of PPT exceeding the threshold value. $c_0$, $c_1$, $c_2$, $d_1$, $d_2$, and $d_3$ are constant, which are 2.52, 0.80, 0.01, 1.43, 0.19, and 0.0013, respectively [54]. In our study, the three different time scales of the SPI (SPI-3, SPI-6, and SPI-12) were calculated. The above three SPIs represent meteorological drought, agricultural drought, and hydrological drought [55].

4.6. Partial Least Squares Regression Model

The partial least squares regression model (PLSR) is a regression modeling method that uses multiple dependent variables and multiple independent variables. It solves the problem of multiple correlations between independent variables in the traditional least squares regression model [56,57]. In our study, we took the WSDI as the dependent variable and hydrogeological factors as independent variables to study the impact of different hydrometeorological factors on droughts in the PRB.

Suppose $Y$ and $X=\left\{X_1,X_2,\cdots ,X_n\right\}$ are the WSDI and hydrometeorological factors, respectively. $X_i\left(i=1,2,\cdots ,n\right)$ represents different hydrometeorological factors, and $n$ is the number of hydrometeorological factors. The time series of $Y$ and $X_i$ can be expressed as $\left\{y_1,y_2,\cdots ,y_p\right\}$ and $\left\{x_{i1},x_{i2},\cdots ,x_{ip}\right\}$, and $p$ is the length of the time series. Therefore, $X$ can be expressed as:

$$X=\left[\begin{array}{ccc}{x}_{11}& \cdots & x_{1n}\\ ⋮& \ddots & ⋮\\ x_{p1}& \cdots & x_{pn}\end{array}\right]$$

(6)

Firstly, $Y$ and $X$ are normalized. Normalized $X$ and $Y$ are denoted as $A$ and $B$, respectively.

Then, the principal components of $A$ and $B$ are extracted, respectively [58]. Assuming that the regression equation meets the accuracy requirements after the $m$th principal component extraction and regression, the extracted principal components can basically represent all independent variables.

Finally, the regression equation of $Y$ with respect to $\left\{X_1,X_2,\cdots ,X_n\right\}$ can be obtained.

The variable importance of projection (VIP) from the PLSR model can accurately quantify the impact of hydrometeorological factors on droughts. The expression of VIP is as follows:

$$VI{P}_j=\sqrt{\frac{n{\displaystyle \sum _{i=1}^{m}Rd\left(Y,v_i\right)w_{ij}^2}}{{\displaystyle \sum _{i=1}^{m}Rd\left(Y,v_i\right)}}}\left(j=1,2,\cdots ,n\right)$$

(7)

where $VI{P}_j$ is the VIP value of the $j$th hydrometeorological factor. $v_i\left(i=1,2,\cdots ,m\right)$ represents the $m$ principal components of $Y$, and $w_{ij}\left(i=1,2,\cdots ,n;j=1,2,\cdots ,p\right)$ represents the weight of each element in the matrix after the matrix $X$ is standardized. $Rd\left(Y,v_i\right)=r^2\left(Y,v_i\right)$ and $r\left(Y,v_i\right)$ are the correlation coefficients between $Y$ and $v_i$, respectively. When VIP is more than 0.8, the corresponding factor has a significant impact [59].

4.7. Correlation Coefficient and Delay Months

Suppose the two independent time series are $x_1$ and $x_2$, the correlation coefficient between the above time series is expressed as [60]:

$$\rho \left(\tau \right)=\frac{{\sigma }_{12}\left(\tau \right)}{\sqrt{{\sigma }_{11}{\sigma }_{22}}}$$

(8)

where $\rho \left(\tau \right)$ is the correlation coefficient, $\tau$ is the delay factor, ${\sigma }_{11}$ and ${\sigma }_{22}$ are the variance levels of $x_1$ and $x_2$, respectively. ${\sigma }_{12}$ is the covariance of $x_1$ and $x_2$. $\left|\rho \left(\tau \right)\right|$ is the maximum ($\left|\rho \left(\tau \right)\right|\le 1$), which means that $\tau$ is the corresponding delay in months ($\left|\tau \right|\le 12$).

5. Results

5.1. Spatiotemporal Distribution of the WSDI in the PRB

We fused the TWSC results from six GRACE solutions to improve the reliability of TWSC results. Figure 2 shows that fused results are in good agreement with TWSC results from six single solutions, and correlation coefficients between fused results and TWSC results from six single solutions are larger than 0.95. The fused results (3.183 mm) offer less uncertainty than the one of the results from six single solutions (CSR-SH, 13.596 mm; GFZ-SH, 15.434 mm, JPL-SH, 12.711 mm; CSR-M, 33.570; JPL-M, 23.902).

Figure 3 compares the WSDI with seven traditional drought indices and presents the corresponding scatter plots. The eight drought indices had similar change trends, and the performance of the WSDI in detecting droughts in the PRB had a high consistency with traditional drought indices. The WSDI had the highest correlation coefficient with the SPEI-06 (0.6917) and the lowest score with SPI-03 (0.5318). The correlation coefficients of SCPDSI, SPEI-03, SPEI-12, SPI-06, and SPI-12 were 0.5688, 0.5848, 06227, 0.6469, and 0.6120, respectively. Overall, the WSDI agreed well with the SCPSDI, the SPEI, and the SPI in drought monitoring in the PRB.

5.2. Drought Characteristics in the PRB

In our study, the WSDI was used to monitor and characterize droughts in the PRB and its sub-basins from 2003 to 2020 (

Table 1. Summary table of drought events identified by the WSDI.

No.	Time Span of Drought Events	Duration (Months)	WSDI			Drought Area Percentage (%)			Local Meteorological Data Validation (Y/N)
			Peak Magnitude	Average Magnitude	Severity	Peak Magnitude	Average Magnitude	Cumulative Magnitude
1	200307–200406	12	−1.7852	−0.9388	−11.2652	100.00	87.60	1051.15	-
2	200410–200505	8	−1.6718	−1.2597	−10.0776	100.00	99.68	797.44	Y
3	200507–200605	11	−1.6282	−0.9810	−10.7907	100.00	94.38	1038.22	Y
4	200609–200705	9	−0.9943	−0.5177	−4.6590	100.00	90.28	812.50	Y
5	200710–200805	8	−1.3530	−0.7105	−5.6841	100.00	90.43	723.45	Y
6	200901–200906	6	−0.8185	−0.3814	−2.2884	100.00	70.20	421.17	Y
7	200909–201005	9	−1.8509	−1.1458	−10.3120	100.00	97.15	874.39	Y
8	201011–201207	21	−1.7394	−0.5912	−12.4149	100.00	88.03	1848.56	Y
9	201305–201310	6	−0.6143	−0.3546	−2.1276	79.49	61.50	368.97	Y
10	201406–201408	3	−0.7469	−0.4754	−1.4262	69.18	60.63	181.90	Y
11	201504–201508	5	−1.8096	−0.8534	−4.2668	100.00	82.13	410.65	Y
12	201706–201711	6	−1.0703	−0.6612	−3.9671	100.00	74.32	445.95	Y
13	201807–201809	3	−0.5566	−0.4366	−1.3099	84.81	70.24	210.71	Y
14	201911–202004	10	−1.0355	−0.4898	−4.8976	100.00	84.99	849.89	Y

and Figure 4). The fourteen droughts occurred in the PRB during the study period. Among them, the worst drought occurred from November 2010 to July 2012, whose peak magnitude and drought severity values were −0.6143 and −12.4149, respectively, which is also the longest one (21 months). During this drought, the HRB was most affected (−14.6244). The seven droughts belonged to the basin-wide one (nos. 1, 2, 3, 5, 7, 8, and 11). During the seven droughts, the most severely affected sub-basins were the WRB (12.7318), the WRB (11.4342), the HRB (14.6244), the WRB (5.3541), the NBPRB (9.5206), the HRB (15.7994), and the ERB (6.8657). During the study period, the WRB and the HRB suffered the most droughts (12). In non-basin-wide droughts, the no. 10 drought (June 2014 to August 2014) only affected two sub-basins (the NBPRB and the YRB). Four droughts mainly occurred in the middle and upper reaches, while three ones mainly appeared in the lower reaches. Among the droughts that occurred in the middle and upper reaches, the most affected sub-basins were all NBPRBs, while the ones that occurred in the lower reaches were all ERBs during the droughts. Therefore, the NBPRB, the WRB, the HRB, and the ERB were found to be the most drought-prone. With 14 droughts, the drought duration varied across sub-basins (Figure 4b). For example, in the no. 8 drought (November 2011 to July 2012), the YRB had the longest drought months (12 months), followed by the NRB (19 months), the HRB (18 months), the WRB (17 months), the ERB (16 months), the NBPRB (15 months), and the PRD (15 months).

Figure 5 shows that all the sub-basins experienced drought nearly or more than half the time (43%). The NRB had the highest number of drought months (52.78%), while the ERB and the PRD were basically in drought half the time (50.46% and 50.93%). The NRB, the ERB, and the PRD belonged to the lower reaches. This suggests that the lower reaches are more susceptible to droughts.

5.3. Impact of Climate Variability on Droughts

Drought is closely related to hydrometeorological conditions. We selected PPT, runoff, SM, GW, and ET as the drought-influencing factors in the PRB for analysis.

Table 2. The correlation coefficient (CC) and VIP between the WSDI and PPT, runoff, SM, GW, and ET anomalies in the PRB.

Variance	PRB		NBPRB		HRB		YRB		WRB		NRB		ERB		PRD
	CC	VIP	CC	VIP	CC	VIP	CC	VIP	CC	VIP	CC	VIP	CC	VIP	CC	VIP
WSDI vs. PPT	0.2603	0.4018	0.1024	0.1438	0.2871	0.4897	0.1640	0.2609	0.2338	0.3969	0.1689	0.3604	0.1627	0.4018	0.1660	0.4359
WSDI vs. runoff	0.2873	0.4896	0.1482	0.3012	0.2713	0.4366	0.2097	0.4263	0.2320	0.3906	0.1635	0.3375	0.1986	0.5984	0.1583	0.3965
WSDI vs. SM	0.6477	2.4888	0.4561	2.8542	0.6750	2.7032	0.5222	2.6439	0.5745	2.3963	0.4128	2.1528	0.3473	1.8299	0.3740	2.2119
WSDI vs. GW	0.5607	1.8651	0.2758	1.0434	0.5621	1.8745	0.4462	1.9302	0.5733	2.3859	0.4670	2.7546	0.4356	2.8784	0.4150	2.7229
WSDI vs. ET	0.3566	0.7543	0.3379	1.5667	0.2890	0.4955	0.2759	0.7377	0.2432	0.4294	0.7166	0.3940	0.1346	0.2751	0.1197	0.2265

shows that SM and GW had significant effects on the drought in the PRB, whose correlation coefficients and VIPs with the WSDI were 0.6477 and 2.4888 and 0.5607 and 1.8651, respectively. However, the situation varied across sub-basins. The HRB, the WRB, and the YRB results were consistent with those of the PRB. Although SM and GW still had significant effects in the NRB, the ERB, and the PRD, and the correlation coefficients between the two and the WSDI were less than 0.5. In the NBPRB, expect for SM and GW, ET had a significant effect on the drought. In the upper reaches of the PRB (the NBPRB, the HRB, and Y the RB), SM was the dominant factor in drought, while GW was the one in the lower reaches of the PRB (the NRB, the ERB, and the PRD). In the WRB, the degree of influence of SM and GW was found to be basically the same.

Due to the delay in response, we calculated the maximum correlation coefficients and the corresponding time delay months between the WSDI and PPT, runoff, SM, GW, and ET anomalies (Figure 6). The spatial correlation distribution of PPT (Figure 6a) shows high distribution characteristics in the middle and low distribution characteristics on both sides, as well as high distribution characteristics in the north and low distribution characteristics in the south. The maximum correlation coefficient was only 0.28, which means that the direct relationship between PPT and the WSDI is not strong. The delay of PPT mainly lasted 1~2 months in most regions (Figure 6b). The spatial distribution of correlation and delay months of runoff was the same as that of PPT (Figure 6c,d). This suggests that PPT and runoff have a close relationship (CC = 0.8682). However, the correlation coefficients of runoff were less than PPT, and the maximum one was 0.25. The correlation between SM and the WSDI was stronger than other factors, and the spatial distribution of correlation of SM was the same as PPT (Figure 6e). This suggests that PPT has an important impact on SM (CC = 0.5663). The delay SM basically lasted 0 months~1 month (Figure 6f). The correlation of GW shows a significant high distribution in the east and low distribution in the west (Figure 6g). This shows that the regions with greater GW impacts are concentrated in the lower reaches, which may be related to the relatively flat terrain in the region, given that the flat terrain is conducive to the recharge and accumulation of GW. The delay of GW almost lasted 0 months. ET had a similar distribution with SM, attributed to SM belong the main source of ET (CC = 0.5274). ET had a strong correlation with the WSDI in the western regions of the NBPRB. Because the eastern NBPRB belongs to the Yunnan-Guizhou Plateau, this region has sufficient sunshine.

With global warming, extreme climates occur frequently, altering the regional climate. Therefore, we studied the teleconnection between the ENSO, the IOD, and droughts in the PRB. The ENSO and the IOD had a positive correlation with PPT in the PRB (CC = 0.2862 and 0.2112). Yuan et al. [61] showed that the western Pacific subtropical high is stronger and more southward in the ENSO year, so the rising movement in South China is strengthened, convection develops, and PPT is increased. Xiao et al. [62] indicated that when the IOD is in a positive phase, North China and the Jiang-Huai River basin are dry and less rainy during the flood reason, and the coastal regions of South China are rainy belts. The ENSO has a greater correlation with PPT than the IOD. This is attributed to the PRB being closer to the Pacific Ocean.

Figure 7 shows the maximum correlation coefficients and the corresponding time delay months between the ENSO and PPT, SM, and the WSDI. The correlations between the ENSO and PPT show high distribution characteristics in the west and low characteristics in the east (Figure 7a). The strongest one occurred in the WRB (0.3334) and the weakest one appeared in the PRD (0.1450). This shows that the ENSO has a greater impact on PPT in the middle and upper reaches than in the lower reaches. Based on the response time of PPT to the ENSO, the closer to the ocean, the faster the response (Figure 7b). The correlation between SM and the ENSO shows high values in the southwest and low values in the northeast (Figure 7c), and the maximum and minimum ones (0.3404 and 0.1737) appeared in the WRB and the YRB, respectively. SM instantly responded to the ENSO in most regions, and the response time in individual regions reached 11 months (Figure 7d). The spatial distribution of SM was quite significantly different from that of PPT. This may have been affected by human activity, such as agricultural irrigation, reservoir regulation, human water discharge, etc. The spatial distribution of the WSDI was consistent with that of PPT (Figure 7e). This indicates that PPT is the main meteorological factor affecting the occurrence of drought in the PRB. The basin most affected by the ENSO was the YRB (0.3181), while the basin with the smallest effect was the ERB (0.01). The response of the WSDI to the ENSO was slow in the northeast and fast in the southwest (Figure 7f). Because the ERB, the NRB, and the PRD belong to the densely populated and economically developed regions of the PRB, human activity may be the main reason for the differences in the responses of PPT and the WSDI to the ENSO.

Figure 8 estimates the maximum correlation coefficients and the corresponding time delay months between the DMI and PPT, SM, and the WSDI. Comparing Figure 8a,e, the spatial distribution of correlations between PPT and the WSDI were similar, i.e., indicating the high distribution characteristics in the northeast and the low distribution characteristics in the southwest. This explains that the impact of the IOD relates to distance in the Indian Ocean, which is the same as that of the ENSO. The response times of PPT and the WSDI to the IOD were significantly different (Figure 8b,f). The correlations between the DMI and SM were high in the west and low in the east, and the low correlations were mainly concentrated in the ERB, the NRB, and the PRD. The reason is the same as for the ENSO. Figure 8d shows that SM responses to the IOD were faster in the east and slower in the west.

5.4. Impact of Human Activity on Droughts

Since the impact of human activity on droughts is difficult to estimate, we used indirect methods for evaluation purposes. We used natural-induced TWSCs to monitor and characterize the droughts in the PRB (

Table A1. Summary table of drought events identified by the WSDI based on natural TWSC.

No.	Time Span of Drought Events	Duration (Months)	WSDI			Drought Area Percentage (%)
			Peak Magnitude	Average Magnitude	Severity	Peak Magnitude	Average Magnitude	Cumulative Manitude
1	200302–200312	11	−1.7283	−0.4644	−5.1088	100.00	79.25	871.72
2	200408–200501	6	−1.0604	−0.5037	−3.0220	100.00	86.00	516.02
3	200507–200511	5	−1.4033	−0.6229	−3.1146	99.36	88.49	442.46
4	200709–200711	3	−0.7144	−0.4736	−1.4208	100.00	84.12	252.35
5	200802–200804	3	−0.3451	−0.2004	−0.6011	98.07	79.02	237.06
6	200812–200902	3	−0.6211	−0.4518	−1.3555	100.00	97.61	292.84
7	200904–200906	3	−0.4907	−0.3164	−0.9491	78.13	67.72	203.16
8	200908–200911	4	−1.5915	−0.7640	−3.0561	98.69	86.89	347.55
9	201101–201108	8	−1.8329	−0.7901	−6.3209	100.00	80.04	640.30
10	201208–201210	3	−1.0076	−0.4514	−1.3541	95.45	73.29	219.86
11	201212–201302	3	−0.5513	−0.2662	−0.7985	99.36	79.14	237.41
12	201502–201504	3	−0.8139	−0.5444	−1.6332	98.73	87.61	262.82
13	201612–201702	3	−0.4576	−0.2658	−0.7974	98.73	80.93	242.79
14	201802–201807	6	−0.6643	−0.3944	−2.3667	99.34	73.35	440.13

). Figure 9 compares the characteristics of droughts caused by natural factors and actual occurrences (drought-nature and drought-GRACE) in the PRB and its sub-basins. Comparing Table 1 and Table A1, the number of drought natures (11) was significantly less than the number of drought GRACE indices (14) in the PRB. This shows that three droughts were caused by human activity. Expect for the ERB, the above situation existed in each sub-basin, and the NBPRB had the most human-induced droughts (3). Figure 9a shows that drought duration extensions appeared in 81.82% of droughts (9) due to human activity. The longest increase in drought duration reached 10 months (no. 8), and the average value was 4 months. The highest increase rate in drought duration was 125% (no. 7) and the average value was 75.37%. The duration of one drought (no. 13) was shortened by 3 months due to human activity, while there was one drought (no. 6) whose duration did not change due to human activity. As shown in Figure 9, due to human activity, the drought with the longest prolonged drought (18 months) occurred in the YRB. The greatest increase in the average drought duration was found in the NBPRB and the ERB (both 6.5 months), and the smallest one was found in the NRB (3.3 months). The drought with the greatest increase rate of the drought duration appeared in the YRB (600.00%). The greatest average increase rate in drought duration was found in the NBPRB (186.31%) and the smallest one was found in the NRB (81.75%). Across the study regions, time extension was the main effect of human activity on drought duration, but duration was also shortened in sporadic droughts.

In the PRB (Figure 9b), the peak magnitude of 81.82% of droughts (9) increased. The largest increase in the peak magnitude was 0.6386 (no. 5), and the maximum growth rate was 133.89% (no. 12) and the minimum growth rate was 3.29% (no. 1). There were two droughts with declines in the peak magnitude, and the reduction rates were 5.10% and 16.21%. However, there are large differences in the results affected by human activity between the sub-basin and the PRB. Except for the WRB and the PRD, at least 50% of drought peaks were weakened in the other five sub-basins. In the NBPRB and the HRB, the peak magnitudes of 55.56% droughts were cut. Droughts with the largest peak reductions occurred in the HRB (91.13%). The seven sub-basins were sorted in ascending order of the average peak magnitude reduction levels, and their arrangement was as follows: the NRB (39.29%), the HRB (37.01%), the WRB (35.37%), the NBPRB (26.71%), the ERB (26.24%), the PRD (20.84%), and the YRB (15.13%). Combining the probability and procedures of peak magnitude reduction, the HRB was found to be the sub-basin with the largest peak magnitude reduction. Human activity had an effect on reducing the peak magnitude and the ability to exacerbate them. According to statistics, the PRD had the highest probability of the exacerbated peak magnitude (70%). The average peak magnitude intensification levels in each sub-basin were as follows: the NBPRB (349.58%), the YRB (118.05%), the HRB (105.75%), the WRB (36.59%), the NRB (56.20%), the ERB (53.40%), and the PRD (140.52%). Although the average level in the NBPRB was greater than that in the PRD, only two droughts in the NBPRB had peak magnitude intensification. In the PRD, there were seven occurrences of this situation. Therefore, the PRD was found to be a sub-basin with the largest peak magnitude intensification.

Except for drought duration and peak magnitude, human activity had a significant impact on drought severity. Figure 9c shows that 90.91% of droughts experienced increased severity in the PRB, and the average growth rate reached 163.77%. The maximum and minimum growth rates were 397.50% (no. 12) and 16.96% (no. 6), respectively. Only one drought occurred with reduced severity (44.65%). Only the HRB, the WRB, the NRB, and the ERB had reduced drought severity, and there were one, one, two, and one droughts in the sub-basins, respectively. The vast majority of droughts increased in drought severity. Among them, the greatest increasing rate was 1928.68% (the NBPRB), while the smallest one was 13.92% (the NRB). The seven sub-basins were sorted in ascending order of the average level of drought severity intensification, and their arrangement was as follows: the NBPRB (385.37%), the YRB (326.16%), the PRD (283.75%), the HRB (252.62%), the NRB (205.59%), the WRB (151.78%), and the ERB (148.54%).

Margariti et al. [63] indicated that human activity can affect drought via water storage and consumption. The water storage was mainly RS, while the water consumption includes industry, agriculture, and domestic water use. In our study, the above water consumptions were considered as a factor, called TWC. Figure 10a shows that TWC (56.30%) was significantly greater than RS (43.70%) during 2003–2013, and the greater proportion of TWC was found in 2013 (72.04%). Since 2011, TWC shows a continuous downward trend, which may be attributed to China’s energy conservation and emission reduction policies implemented since 2010 [64]. However, the average proportion of TWC still reached 51.02% during 2011–2022. TWC mainly plays a role in exacerbating the drought, so human activity mainly aggravates droughts in the PRB, which is consistent with the analysis results shown in Figure 9a–c.

The purpose of RS is not only flood control and drought relief, but also water supply and power generation, so RS has both positive and negative effects on droughts. Due to human activity, the drought severity of no. 13 in the PRB decreased by 44.65% (Figure 9c). This shows that RS played a role in mitigating drought in this case and offset the negation effect of TWC. Previous studies indicate [65,66] that when terrestrial water storage is abundant, the reservoir increases the water storage capacity, and when the drought occurs, stored water is released to alleviate the drought. During the study period, RS showed an increasing trend, from 23.23 billion cm³ to 70.65 billion cm³. However, the drought did not appear to significantly ease. On the one hand, this is because the proportion of RS was always smaller than that of TWC; on the other hand, RS was not all used for drought protection.

Figure 10f shows that the average proportion of RS in the PRD was the smallest (7.29%); thus, drought exacerbations caused by human activity were the most severe (drought duration, peak magnitude, and drought severity). Figure 10b,c show that the NBPRB and the HRB had the greater proportion of RS (63.20%) and the greatest reservoir storage (12.09 billion m³). Although there was no significant performance in terms of drought duration and severity, the above two sub-basins were found to be the most effective in the terms of peak magnitude reduction. This suggests that relying solely on RS to mitigate the effects of drought is not feasible in the PRB, and various measures must be integrated to reduce the adverse impacts of drought on the social economy.

6. Discussion

The PRB belongs to the Subtropical Monsoon Climate Zone, and the Indian Ocean and the South China Sea are the main sources of water vapor transportation in the region. Therefore, the East Asian Summer Monsoon (EASM) and the Indian Summer Monsoon (ISM), which originate from the Pacific and Indian Oceans, respectively, are bound to affect climate anomalies in the PRB. Huang et al. [67] indicate that the spatiotemporal evolution characteristics of PPT in the PRB are significantly coupled with anomalous changes in the sea surface temperature (SST) in the Pacific and Indian Oceans at the equator. Previous studies [68,69,70] show that the ENSO affects anomalous changes in atmospheric circulation by affecting the Walker and Hadley circulations. During the El Niño event, Walker circulation weakens and Hadley circulation increases, which causes the EASM to weaken, and PPT is mainly concentrated in South China. When a positive IOD appears, the equatorial region of the Indian Ocean is anomalously easterly, and there is an obvious westerly response from the eastern Bay of Bengal to the Indochina Peninsula, which strengthens the ISM. However, there is a significant westerly wind response in the western Pacific region, which weakens the water vapor transportation from the tropical ocean. This causes the PPT belt to recede to South China. In most cases, the ENSO and the IOD are concomitant [53]. At this time, there is an obvious easterly response from the eastern Bay of Bengal to the Indochina Peninsula. This hinders the northward movement of the high-altitude easterly belt. As a result, the southwest monsoon in East Asia continues to weaken, and the water vapor cannot be transported to North China, so a large amount of water vapor stays in South China [68]. Therefore, during El Niño and positive IOD events, the PRB is prone to floods and vice versa; thus, it is prone to droughts. In our study, PPT is more affected by the ENSO, which is consistent with the results of Huang et al. [30] and Deng et al. [46].

With rapid increases in the population and the industrialization process, people endlessly demand a large amount of water from nature, which breaks the water balance in nature, disrupts the terrestrial water cycle, and increases the intensity and frequency of extreme droughts. As shown in Figure 10, the YRB (57.56%), the WRB (81.12%), the NRB (72.26%), and the PRD (92.71%) are all regions with a large proportion of TWC. Therefore, the above four sub-basins experienced relatively more droughts (more than 11) during the study period. The NBPRB (36.80%) and the ERB (38.05%) were found to have a smaller proportion of TWC, so the drought numbers in these two sub-basins were less than 10. This shows that the TWC exacerbates droughts, which has good agreement with the findings of Wang et al. [71] centered around the fact that irrigation projects and agricultural production in the PRB increased the risk of droughts. In addition, RS can change the spatiotemporal distribution of a TWSC region, which may affect the peak magnitude, severity, and duration of droughts. Yang et al. [72] indicate that if RS is large and the proportion of TWC is small, the droughts will be alleviated, and the opposite will aggravate the droughts. Zhang et al. [65] suggest that RS plays a role in smoothing runoff change, effectively reducing the drought severity, but prolonging the drought duration.

Our results comprehensively quantify the extent to which climate variability and human activity have influenced the drought in the PRB, especially a precise quantification of the impact of RS and TWSC on drought characteristics. We used an indirect method to quantify the impact of human activity on droughts, which has definite uncertainty. In the follow-up research, it is hoped that high-precision monthly human activity data can be obtained to accurately quantify the impact of human activity on droughts.

7. Conclusions

In our study, a combined drought index based on GRACE and GRACE-FO solutions was used to analyze the drought characteristics and influencing factors in the PRB and its sub-basins during 2003 and 2020. We also combined multiple hydroclimate data, climate indexes, and human activity data to study the impact of climate variability and human activity on droughts. Our conclusions are as follows:

(1)

Although 14 droughts occurred in the PRB, not all droughts are basin-wide. Among them, the worst drought occurred from Nov 2010 to July 2012. There are significant differences in the drought characteristics of each sub-basin. The NRB is the basin that was affected by droughts for the longest durations.

(2)

The middle and upper reaches of the PRB were mainly dominated by SM, while the lower reaches were mainly dominated by GW. The ENSO had a greater impact on PPT in the PRB than the IOD. The impacts of the ENSO and the IOD on PPT present significant spatial east–west differences.

(3)

In the PRB, human activity leads to longer drought durations and more severe drought, but can weaken the peak magnitude. The PRD is the basin most affected by human activity. Due to TWC, the drought frequency increased. RS can play a role in reducing the peak magnitude and drought severity, but this depends on the use of a reservoir.

Our study helps to establish a connection mechanism between climate variability, human activity, and regional droughts, which is of great scientific value and social significance for the realization of early warning signs of regional droughts and the assessment of drought impacts.