Spatiotemporal Characteristics of Drought and Driving Factors Based on the GRACE-Derived Total Storage Deﬁcit Index: A Case Study in Southwest China

: Drought monitoring is useful to minimize the impact of drought on human production and the natural environment. Gravity Recovery and Climate Experiment (GRACE) satellites can directly capture terrestrial water storage anomalies (TWSA) in the large basin, which represents a new source of hydrological information. In this study, the GRACE-based total storage deﬁcit index (TSDI) is employed to investigate the temporal evolution and spatial distribution of drought in Southwest China from 2003 to 2016. The comparison results of TSDI with the standardized precipitation index (SPI), the standardized precipitation evapotranspiration index (SPEI), and the self-calibrating Palmer drought severity index (SC-PDSI) show that TSDI has signiﬁcant consistency with them, which veriﬁes the reliability of TSDI. The spatial distribution of TSDI was more consistent with the governmental drought reports than SC-PDSI in the most severe drought event from September 2009 to April 2010. Finally, the links between drought and climate indicators are investigated using the partial least square regression (PLSR) model. The results show that insufﬁcient precipitation has the most signiﬁcant impact on drought in Southwest China, followed by excessive evaporation. Although Southwest China is selected as a case study in this paper, the method can be applied in other regions as well.


Introduction
Drought is one of the most devastating natural disasters in the world, which will bring damage to agricultural production, social economy, and ecological environment [1,2]. The Intergovernmental Panel on Climate Change 2018 report warned that human activities have caused global temperatures to rise by~1 • C since industrialization, and global warming could reach 1.5 • C between 2030 and 2050 [3]. Due to global climate change, drought events in many areas are expected to intensify in the 21st century [4,5]. Therefore, it is urgent to carry out effective drought monitoring for the sake of reducing the damage of drought. Traditional drought monitoring approaches rely on meteorological and hydrological data observed by stations, and the spatial distribution of drought is realized through interpolation [6]. However, the drought monitoring result obtained by interpolation is inaccurate in areas with complex terrain or scarcity of stations [7,8].
With the development of remote sensing technology, drought monitoring has got rid of the dependence on traditional station-observed data, and can obtain the critical drought-related variables over a larger spatial and temporal scale. Gravity Recovery and Climate Experiment (GRACE) gravity satellite mission, launched in March 2002, is cooperatively developed by the National Aeronautics and Space Administration (NASA) and the German Aerospace Center. The GRACE twin satellites can capture terrestrial water storage anomalies (TWSA) by monitoring the changes in Earth's gravity field, which represents an alternative source of hydrological information [9,10]. GRACE-derived TWSA is vertically divided into five parts: surface water, ice and snow, soil moisture, groundwater, and water contained in biomass, which reflects all water storage information [11,12]. The development of GRACE realizes the dynamic monitoring of water storage information in the large-scale region, and shifts drought monitoring away from reliance on traditional station-based measurements. Therefore, GRACE has excellent potential in drought monitoring and has been applied in regional drought investigation. Using the GRACE-based total storage deficit index (TSDI), Yirdaw et al. [13] generated a pictorial representation of the long-term dryness and wetness in the Canadian Prairie. Yi and Wen [14] introduced a GRACE-based hydrological drought index for drought monitoring in the continental United States from 2003 to 2012. Thomas et al. [15] established the GRACE groundwater drought index for drought monitoring in the Central Valley of California, and captured the characteristics of groundwater drought due to complex human activities and natural changes. Yu et al. [16] evaluated the drought using GRACE water storage deficit index over Mongolia, and demonstrated the effectiveness of GRACE in analyzing and assessing the severity of drought in areas lacking hydrological observations.
Under the background of global climate change, the frequency of drought events has increased in Southwest China in the past few decades [17]. Extreme drought often leads to a large area of crop failure and water shortage for people and animals, causing severe damage to the local area [8]. Specifically, Sichuan province and Chongqing municipality in Southwest China suffered the worst drought in the summer of 2006 since records began in 1891, resulting in drinking water difficulties for more than 15 million people [18]. Besides, Southwest China suffered a "once-in-a-century" extreme drought from autumn 2009 to spring 2010, causing economic damage equal to 2.13% of the gross domestic product in 2010 [19]. Therefore, it is of considerable significance and practical value to study the drought in Southwest China.
This present study aims to explore the drought situation by the GRACE-derived TSDI in Southwest China during 2003-2016, and to reveal the relationship between the drought identified by GRACE and climate factors. The primary objectives of this work are: (1) to verify the reliability of TSDI for evaluating drought events in Southwest China; (2) to explore the temporal evolution and spatial distribution of drought using GRACE-based TSDI; and (3) to analyze the links between TSDI and climate indicators using the partial least square regression (PLSR) model.

Study Area
Southwest China between 97-112 • E and 21-35 • N, with approximately 1.36 × 10 6 km 2 ( Figure 1a). This region consists of the municipality of Chongqing and the four provinces Sichuan, Yunnan, Guizhou, and Guangxi ( Figure 1b). There are 138 meteorological stations in the study area, and the station density in the east is significantly higher than that in the west (Figure 1c). Southwest China is mainly covered by grassland and forestland, and the agriculture land is concentrated in the Sichuan Basin (Figure 1d). The region has a subtropical monsoon climate, and the annual average temperature and precipitation arẽ 16 • C and~1100 mm, respectively [20]. Water resources are abundant in Southwest China, and the Yangtze River, the Pearl River, and the Lancang River pass through the region. However, the spatial and temporal distribution of precipitation is extremely uneven, and 85% of the precipitation is concentrated in April to October and decreases from southeast to northwest, which is conducive to the formation of drought [21]. Affected by climate change, Southwest China has suffered a frequent drought in recent years, which leads to a large area of crop failure and water shortage for people [22].
precipitation are ~16 °C and ~1100 mm, respectively [20]. Water resources are abundant in Southwest China, and the Yangtze River, the Pearl River, and the Lancang River pass through the region. However, the spatial and temporal distribution of precipitation is extremely uneven, and 85% of the precipitation is concentrated in April to October and decreases from southeast to northwest, which is conducive to the formation of drought [21]. Affected by climate change, Southwest China has suffered a frequent drought in recent years, which leads to a large area of crop failure and water shortage for people [22].

GRACE Data
GRACE data provide the basis for understanding the terrestrial water cycle, ice sheet, and glacier mass balance, as well as their links to global climate change [23]. The mascon solutions can directly calculate the equivalent water height corresponding to the regional mass change [24]. The TWSA is the anomaly of the equivalent water height relative to the mean baseline from 2004 through 2009. In contrast to the traditional spherical harmonic solution, the mascon solution can reduce the signal leakage errors between the ocean and the land, and no additional decorrelation filtering is required to remove the north-south stripe error [25]. CSR RL06 mascon products (http://www2.csr.utexas.edu/grace/RL06_mascons.html) are the latest solutions provided by the university of Texas at Austin, Center for Space Research (UTCSR). The CSR RL06 mascon products, with a 0.25° × 0.25° equiangular grid and a monthly temporal resolution,

GRACE Data
GRACE data provide the basis for understanding the terrestrial water cycle, ice sheet, and glacier mass balance, as well as their links to global climate change [23]. The mascon solutions can directly calculate the equivalent water height corresponding to the regional mass change [24]. The TWSA is the anomaly of the equivalent water height relative to the mean baseline from 2004 through 2009. In contrast to the traditional spherical harmonic solution, the mascon solution can reduce the signal leakage errors between the ocean and the land, and no additional decorrelation filtering is required to remove the north-south stripe error [25]. CSR RL06 mascon products (http://www2.csr.utexas.edu/grace/RL0 6_mascons.html) are the latest solutions provided by the university of Texas at Austin, Center for Space Research (UTCSR). The CSR RL06 mascon products, with a 0.25 • × 0.

Meteorological Data
Daily meteorological data are provided by the National Meteorological Information Center of the China Meteorological Administration (http://data.cma.cn). Meteorological data used in this study include evaporation (EVP), precipitation (PRE), air temperature (TEM), air pressure (PRS), relative humidity (RHU), sunshine duration(SSD), wind speed (WIN), and ground temperature (0 cm) (GST). The comparison of EVP data provided by the station with evapotranspiration data from Moderate-resolution Imaging Spectroradiometer (MODIS), Global Land Data Assimilation System (GLDAS), and Global Land Evaporation Amsterdam Model (GLEAM) confirms the reliability of site-based EVP data ( Figure S1). Daily evaporation and precipitation data are accumulated to obtain the monthly precipitation records, and other daily data are averaged to get the monthly data. The stations established after 1966 are excluded for ensuring the time length of meteorological data as complete as possible, and a total of 138 meteorological stations remain in the study area. The data time series is from 1966 to 2016 and provided by meteorological stations established before 1966.

Methods
A schematic illustration of the methodology steps is provided in Figure 2. Firstly, the TSDI is calculated using GRACE data to assess the drought in Southwest China. The standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI) are calculated from meteorological data, and the self-calibrating Palmer drought severity index (SC-PDSI) datasets are provided by the Climatic Research Unit (CRU) (https://crudata.uea.ac.uk/cru/data/drought/). Secondly, the reliability of TSDI is verified by spatiotemporal comparison of TSDI with standardized drought indexes (SPI, SPEI, and SC-PDSI). Finally, the links between the TSDI and climate indicators are established by the PLSR model, and the main factors affecting drought in Southwest China are analyzed based on the variable importance of the projection (VIP) from the PLSR model.

Meteorological Data
Daily meteorological data are provided by the National Meteorological Information Center of the China Meteorological Administration (http://data.cma.cn). Meteorological data used in this study include evaporation (EVP), precipitation (PRE), air temperature (TEM), air pressure (PRS), relative humidity (RHU), sunshine duration(SSD), wind speed (WIN), and ground temperature (0 cm) (GST). The comparison of EVP data provided by the station with evapotranspiration data from Moderate-resolution Imaging Spectroradiometer (MODIS), Global Land Data Assimilation System (GLDAS), and Global Land Evaporation Amsterdam Model (GLEAM) confirms the reliability of sitebased EVP data ( Figure S1). Daily evaporation and precipitation data are accumulated to obtain the monthly precipitation records, and other daily data are averaged to get the monthly data. The stations established after 1966 are excluded for ensuring the time length of meteorological data as complete as possible, and a total of 138 meteorological stations remain in the study area. The data time series is from 1966 to 2016 and provided by meteorological stations established before 1966.

Methods
A schematic illustration of the methodology steps is provided in Figure 2. Firstly, the TSDI is calculated using GRACE data to assess the drought in Southwest China. The standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI) are calculated from meteorological data, and the selfcalibrating Palmer drought severity index (SC-PDSI) datasets are provided by the Climatic Research Unit (CRU) (https://crudata.uea.ac.uk/cru/data/drought/). Secondly, the reliability of TSDI is verified by spatiotemporal comparison of TSDI with standardized drought indexes (SPI, SPEI, and SC-PDSI). Finally, the links between the TSDI and climate indicators are established by the PLSR model, and the main factors affecting drought in Southwest China are analyzed based on the variable importance of the projection (VIP) from the PLSR model.

GRACE-based TSDI
Yirdaw et al. [13] renamed the soil moisture deficit index developed by Narasimhan and Srinivasan [26] as TSDI. Drought events may occur when an area is continuously dry, and the GRACE-based TSDI is employed to describe this long-term dryness [27]. Firstly,

GRACE-Based TSDI
Yirdaw et al. [13] renamed the soil moisture deficit index developed by Narasimhan and Srinivasan [26] as TSDI. Drought events may occur when an area is continuously dry, and the GRACE-based TSDI is employed to describe this long-term dryness [27]. Firstly, the total storage deficit (TSD, %) is calculated using TWSA, as shown in the following formula [13]: where TSA i,j denotes the GRACE-derived TWSA in month j of the year i. MTSA j , MaxTSA j , and MinTSA j represent the long-term mean, maximum and minimum values of total storage anomaly (TSA) in month j. Drought is a long-term process, and the current month's drought is not only related to the current TSD but also previous drought conditions. Therefore, TSDI is calculated based on past drought conditions and current TSD, as follows [13]: According to the cumulative TSD map during drought, parameters p and q can be determined as described in the following formula [13]: where m and b stand for the slope and intercept of the best-fitting line of cumulative TSD curve in the dry period, respectively. Parameter C is an important parameter for calculating TSDI, and represents the drought intensity during the selected dryness period (the period corresponding to the bestfitting line of the cumulative TSD). Unlike the traditional TSDI method regarding parameter C as an empirical parameter, parameter C is determined by SPI and SPEI calculated from the meteorological data of ground stations in this study, as follows: where C dsi stands for the drought intensity calculated by the drought severity index, D i represents the drought intensity obtained from the standardized drought index calculated from the meteorological data of the station, and P i is the percentage of stations with D i drought intensity. According to the definition of Palmer, the drought intensity D i is divided into five categories: -4 for extreme drought, -3 for severe drought, -2 for moderate drought, -1 for mild drought and 0 for no drought.
In previous studies, Parameter C is determined by previous research records, limiting the application of TSDI in areas lacking drought research records. In this study, the parameter C is determined by SPI and SPEI calculated from the data of meteorological stations, which makes the application of TSDI is no longer limited to the areas with abundant research records.

Standardized Drought Indices
Mckee et al. [28] proposed the meteorological drought index SPI, which is based on long-term precipitation records to characterize the probability of precipitation in a certain period. SPI is one of the most widely used drought indicators in the world because of its simple calculation and reliable results. In this study, the monthly precipitation data of 51 years were used to calculate the SPI of 3 months, 6 months, and 12 months, respectively. Firstly, the long-term precipitation records are first fitted into the gamma probability density function to calculate the SPI. Then the probability density function is transformed into a cumulative probability function. Finally, the drought grade is classified according to the cumulative frequency distribution of standardized precipitation, as follows [29]:  [29]. Vicente-Serrano et al. [30] introduced the standardized drought index SPEI, which considers the effect of temperature based on the SPI. The method to calculate SPEI is the same as SPI, but the cumulative distribution function for precipitation minus the potential evapotranspiration values. Both SPI and SPEI calculate three-time scales of 3month, 6-month, and 12-month. The 3-month time scale is mainly related to soil moisture (meteorological drought), 6-month time scale is mainly related to reservoir storages and river discharge (agricultural drought), and 12-month time scale is mainly related to the change of the groundwater storage (hydrological drought). Therefore, SPI and SPEI on different time scales can be useful for monitoring different types of drought.
The Palmer Drought Severity Index is a standardized drought index that comprehensively considers the impact of hydrological and meteorological variables on drought [31]. The original PDSI is limited by fixed climate weighting factors and duration factors, which makes it not comparable in different climate regions [32]. SC-PDSI replaces the fixed parameters with the parameters calculated by the climate data of each region, which realizes the comparison of the drought assessment results in different regions [33]. In this study, the SC-PDSI datasets, with a spatial resolution of 0.5 • , are provided by the CRU.
A unified drought classification standard should be established to compare different drought indices. The severity of drought is divided into four levels (D1, D2, D3, D4), while D0 indicates that the area has not suffered from drought, and D4 stands for extreme drought [34]. The drought severity category of four drought indices, including TSDI, SPI, SPEI, and SC-PDSI, is shown in Table 1.

Category
Description

Seasonal-trend Decomposition by Loess Method
The seasonal-trend decomposition by loess (STL) method, which is based on locally weighted regression, is considered to be robust and computationally efficient. For detailed principles of STL, the reader can refer to Cleveland et al. [35]. The STL can be used to decompose the original time series (Y) into trend (T), seasonal (S), and residual components (R), as follows [35]: In the previous study, the seasonal-trend decomposition method was used to decompose GRACE data to obtain long-term trend signals and seasonal changes [36,37]. The season has a significant influence on climatic factors, while drought is usually not related to the season. The correlation between all original climatic factors and TSDI is very low, bringing significant uncertainty to construct the PLSR model by using original climatic factors and TSDI. The climatic factors are decomposed by the STL, and then the seasonal and residual terms are removed, and only the annual trend is retained, which can greatly improve the correlation between climatic factors and TSDI (Table S1). Previous studies have reported the feasibility of using the STL method to decompose variables and then building a model to reconstruct GRACE data [38,39]. In this work, the STL method is used to extract the residuals from GRACE mascon time series for uncertainty estimation, and to extract interannual trends from TSDI and climate indicators to reflect the relationship between them.

The Partial Least Square Regression Model
In this study, the partial least squares regression (PLSR) model was used to analyze the contribution of different climate indicators to drought in Southwest China. PLSR is a linear model illustrating the relationship between dependent variable Y and a set of independent variables X, as follows [40]: where a 0 is the intercept, and a n is the regression coefficients calculated based on the data. The PLSR model finds out the independent variables which have the most significant influence on the dependent variable by extracting the data and screening the components [41,42]. The variable importance of the projection (VIP) from the PLSR model can quantitatively explain the impact of the independent variables on the dependent variable [43]. It is generally believed that the independent variable with VIP > 0.8 is meant to explain the dependent variable [40]. In this study, the independent variables of PLSR are EVP, PRE, TEM, PRS, RHU, SSD, WIN, and GST data, and the dependent variable is the drought index. All variables are decomposed by the STL method, and the seasonal component and residual component are removed. The VIP is used to determine which climate forcing factors have the most significant contribution to the development of drought.

Drought Events Detected by TSDI
Natural influences and human activities can lead to a shortage of groundwater, which is often a sign of a drought event. During the study period, the TWSA reveals apparent seasonal variation, falling to the trough from March to May, and rising to the peak from August to October ( seasonal and residual terms are removed, and only the annual trend is retained, which can greatly improve the correlation between climatic factors and TSDI (Table S1). Previous studies have reported the feasibility of using the STL method to decompose variables and then building a model to reconstruct GRACE data [38,39]. In this work, the STL method is used to extract the residuals from GRACE mascon time series for uncertainty estimation, and to extract interannual trends from TSDI and climate indicators to reflect the relationship between them.

The Partial Least Square Regression Model
In this study, the partial least squares regression (PLSR) model was used to analyze the contribution of different climate indicators to drought in Southwest China. PLSR is a linear model illustrating the relationship between dependent variable Y and a set of independent variables X, as follows [40]: where a0 is the intercept, and an is the regression coefficients calculated based on the data. The PLSR model finds out the independent variables which have the most significant influence on the dependent variable by extracting the data and screening the components [41,42]. The variable importance of the projection (VIP) from the PLSR model can quantitatively explain the impact of the independent variables on the dependent variable [43]. It is generally believed that the independent variable with VIP > 0.8 is meant to explain the dependent variable [40]. In this study, the independent variables of PLSR are EVP, PRE, TEM, PRS, RHU, SSD, WIN, and GST data, and the dependent variable is the drought index. All variables are decomposed by the STL method, and the seasonal component and residual component are removed. The VIP is used to determine which climate forcing factors have the most significant contribution to the development of drought.

Drought Events Detected by TSDI
Natural influences and human activities can lead to a shortage of groundwater, which is often a sign of a drought event. During the study period, the TWSA reveals apparent seasonal variation, falling to the trough from March to May, and rising to the peak from August to October (Figure 3). According to the TWSA trend line in Figure    After calculating the TSDI time series, the TDSI results are compared with SPI, SPEI, and SC-PDSI. As shown in Figure 5, TSDI is consistent with SPI, SPEI, and SC-PDSI in the description of drought in Southwest China, with similar peaks and troughs, and all the drought indices have a large trough in 2009-2010 when the GRACE-derived TWSA have the minimum values. However, there are also differences in the assessment of drought by these drought indices because they are developed using different variables and approaches. For example, TSDI is consistent with the standardized drought indices during 2003-2012, but is generally higher than other indices after 2012. Figure 6 demonstrates the significant correlation between TSDI and other drought indices over Southwest China from 2003 to 2016 (at a 95% confidence interval). TSDI has the highest Parameter C is determined jointly by SPI-6 (SPI with a 6-month time scale) and SPEI-6 (SPEI with a 6-month time scale) from September 2009 to April 2010. The SPI-6 indicates that 92.75% of the stations suffer drought, while SPEI-6 shows that 97.83% of the stations are attacked by drought (Figure 4c). Two drought indices illustrate that nearly half of the stations have suffered extreme drought events (D4). The SPI-6 indicates that 20.29% of the stations suffer severe drought (D3), and this proportion reaches 28.26% in terms of SPEI-6. According to Equation (4), the values of C SPI-6 and C SPEI-6 are −2.90 and −3.10, respectively, so the value of parameter C is considered to be −3. The best-fitting line of the cumulative TSD (solid black line in Figure 4b) is defined as the upper limit of severe drought. The horizontal zero line in the drought monograph is defined as no drought, and then the interval from no drought to severe drought is divided into three equal intervals. The two dashed lines marked "-1" and "-2" above the best-fitting line represent the upper limits of mild drought (red line) and moderate drought (blue line), respectively. The fourth line is added below the best-fitting line at equal intervals to represent the upper limit of extreme drought (yellow line).
According to the values of parameters m, b, and C, it can be calculated that the parameter p is -0.15 and q is 0.15. By substituting the values of p and q into formula Equation (2), the TSDI of the study area was obtained, as shown in follow: where the initial TSDI is obtained by multiplying TSD 1 by 2% [26] As shown in Figure 4d, the TSDI in Southwest China fluctuated from -4.56 to 5.41 from 2003 to 2016. When the TSDI value of an area was less than -1 for three consecutive months or more, the area is considered to have a drought event [44]. Therefore, seven drought events were detected between 2003 and 2016, as shown in Table 2. The drought Remote Sens. 2021, 13, 79 9 of 19 degree of each drought event can be obtained by calculating the best-fitting line of the cumulative TSDI. After calculating the TSDI time series, the TDSI results are compared with SPI, SPEI, and SC-PDSI. As shown in Figure 5, TSDI is consistent with SPI, SPEI, and SC-PDSI in the description of drought in Southwest China, with similar peaks and troughs, and all the drought indices have a large trough in 2009-2010 when the GRACE-derived TWSA have the minimum values. However, there are also differences in the assessment of drought by these drought indices because they are developed using different variables and approaches. For example, TSDI is consistent with the standardized drought indices during 2003-2012, but is generally higher than other indices after 2012. Figure 6 demonstrates the significant correlation between TSDI and other drought indices over Southwest China from 2003 to 2016 (at a 95% confidence interval). TSDI has the highest correlation coefficient of 0.74 with SPI-6 and the lowest correlation coefficient of 0.59 with SPEI-3. SPI and SPEI have a high correlation at the same time scale, and SC-PDSI is better correlated with SPI and SPEI on a 12-month time scale, with correlation coefficients of 0.84 and 0.81, respectively. Overall, TSDI is consistent with SPI, SPEI, and SC-PDSI in drought monitoring over Southwest China, which verifies the reliability of GRACE-based TSDI.
correlation coefficient of 0.74 with SPI-6 and the lowest correlation coefficient of 0.59 with SPEI-3. SPI and SPEI have a high correlation at the same time scale, and SC-PDSI is better correlated with SPI and SPEI on a 12-month time scale, with correlation coefficients of 0.84 and 0.81, respectively. Overall, TSDI is consistent with SPI, SPEI, and SC-PDSI in drought monitoring over Southwest China, which verifies the reliability of GRACE-based TSDI.

Spatial Distribution of Drought
The spatial distribution of drought obtained by TSDI is also depicted to reflect the initiation and termination of the drought event during 2009.09-2010.04 (Figure 7a). During the drought event, Yunnan is covered by severe and extreme drought throughout, while Chongqing suffered little damage from drought. In November 2009, when the

Spatial Distribution of Drought
The spatial distribution of drought obtained by TSDI is also depicted to reflect the initiation and termination of the drought event during 2009.09-2010.04 (Figure 7a). During the drought event, Yunnan is covered by severe and extreme drought throughout, while Chongqing suffered little damage from drought. In November 2009, when the drought is the worst, the average TSDI values of Sichuan, Chongqing, Guizhou, Guangxi, and Yunnan are -3.22, -1.22, -2.85, -4.34, and -5.98, respectively. Besides, Figure S2

The Links between TSDI and Climate Factors
Drought is closely related to climatic factors, so it is necessary to assess the relationship between drought and climatic factors for better understanding the causes and development of drought from the perspective of climate change. In this study, the PLSR model is used to simulate the links between climate factors and drought, and the impact of various climate indicators on drought development is evaluated based on the VIP obtained by the PLSR model. Before constructing the PLSR model, both climate factors and TSDI were decomposed by STL to obtain the inter-annual trend, and the results are shown in the Figure S3. The simulation result of the PLSR model has sufficiently high goodness of fit with the determination coefficient of 0.80, which indicates that meteorological factors are the dominant factor of drought in Southwest China. As shown in Figure 8, the values of VIP greater than 0.8 are PRE (1.52), EVP (1.25), RHU (1.14), PRS (1.06), and WIN (0.87), while the VIP of TEM and GST are lower, with the values of 0.41 and 0.40, respectively. Besides, at 95% confidence level, there is a significant positive correlation between precipitation and TSDI series (correlation coefficient = 0.78) and a significant negative correlation between evaporation and TSDI series (correlation As shown in Figure 7, the spatial distribution of TSDI and SC-PDSI is also compared. The spatial distribution of TSDI reveals the process of drought spreading from Yunnan to other areas during 2009.09-2010.04, which is consistent with the record of the government-issued bulletin of flood and drought disasters in China [19]. However, SC-PDSI underestimates the severity of drought at the beginning in Yunnan and Guizhou and overestimates the severity of drought at the end in Guangxi. TSDI shows a more intuitive and accurate scope and severity of drought than SC-PDSI, and more consistent with government reports.

The Links between TSDI and Climate Factors
Drought is closely related to climatic factors, so it is necessary to assess the relationship between drought and climatic factors for better understanding the causes and development of drought from the perspective of climate change. In this study, the PLSR model is used to simulate the links between climate factors and drought, and the impact of various climate indicators on drought development is evaluated based on the VIP obtained by the PLSR model. Before constructing the PLSR model, both climate factors and TSDI were decomposed by STL to obtain the inter-annual trend, and the results are shown in the Figure S3. The simulation result of the PLSR model has sufficiently high goodness of fit with the determination coefficient of 0.80, which indicates that meteorological factors are the dominant factor of drought in Southwest China. As shown in Figure 8,

Drought Severity Evaluation
Southwest China has suffered more frequent extreme disasters over the past few decades [45]. The drought events detected by the TSDI in 2006 and 2009-2010 reach extreme drought levels in the worst months (Table 2). Zhang and Zhou [46] also confirmed that extreme drought events occurred in 2006 and 2009-2010 in their retrospective analysis. It is worth noting that the estimated overall drought level in 2006 is only moderate drought due to the overestimation of the duration of drought in 2006 by TSDI, which is lower than the actual drought severity. The deviation of the drought investigate may be due to the deviation caused by insufficient time series of GRACE data [44]. One crucial measure of drought events is severity, and the other is spatial distribution. The GRACE-based TSDI detects that the drought from 2009 to 2010 is mainly concentrated in Yunnan, Guizhou, and the northwest of Guangxi, which is consistent with the spatial distribution of this severe drought event analyzed by Zhao et al. [22]. TSDI clearly shows the worst areas of the drought, and the spatial distribution detected is consistent with other studies and government reports.
To better assess the severity of drought in different periods, Figure 9 demonstrates the percentage of areas affected by the drought of different severity in Southwest China each month. The area most affected by drought is the drought event from September 2009 to April 2010, and 33% of Southwest China area suffer from extreme drought (D4), 13% experience severe drought (D3), 14% suffer from moderate drought (D2), 19% experience

Drought Severity Evaluation
Southwest China has suffered more frequent extreme disasters over the past few decades [45]. The drought events detected by the TSDI in 2006 and 2009-2010 reach extreme drought levels in the worst months (Table 2). Zhang and Zhou [46] also confirmed that extreme drought events occurred in 2006 and 2009-2010 in their retrospective analysis. It is worth noting that the estimated overall drought level in 2006 is only moderate drought due to the overestimation of the duration of drought in 2006 by TSDI, which is lower than the actual drought severity. The deviation of the drought investigate may be due to the deviation caused by insufficient time series of GRACE data [44]. One crucial measure of drought events is severity, and the other is spatial distribution. The GRACE-based TSDI detects that the drought from 2009 to 2010 is mainly concentrated in Yunnan, Guizhou, and the northwest of Guangxi, which is consistent with the spatial distribution of this severe drought event analyzed by Zhao et al. [22]. TSDI clearly shows the worst areas of the drought, and the spatial distribution detected is consistent with other studies and government reports.
To better assess the severity of drought in different periods, Figure 9

Influence Factors of Drought
Under the background of the aggravation of global climate change, the accurate analysis of the influencing factors of drought is of great significance to the comprehensive management of drought [47]. In this study, the VIP results based on the PLSR model determine that the abnormal decrease of precipitation and the increase of evaporation are the most critical factors leading to drought in Southwest China. The shortage of precipitation will lead to insufficient water supply, which directly hinders the growth of vegetation and human production, and excessive evaporation will further accelerate the loss of available water resources and destroy the balance of the water cycle [48]. There are only three years in Southwest China with an average annual precipitation of less than 1000 mm since 1966, namely, 2006 (997.14 mm), 2009 (954.97 mm), and 2011 (910.28 mm), and TSDI has detected moderate or severe drought in these three years. Similarly, the evaporation in 2006, 2009, and 2011 is higher than the average value from 2000 to 2016 by 61.73 mm, 48.17 mm, and 22.38 mm, respectively. In addition, the STL method is used to extract the annual trends of TSDI and climate indicators, and the relationship between them is reflected by comparison. As shown in Figure 10, the abnormal decrease of precipitation and the increase of evaporation often correspond to the trough of TSDI, which further confirms that the accuracy of PLSR model simulation results.

Influence Factors of Drought
Under the background of the aggravation of global climate change, the accurate analysis of the influencing factors of drought is of great significance to the comprehensive management of drought [47]. In this study, the VIP results based on the PLSR model determine that the abnormal decrease of precipitation and the increase of evaporation are the most critical factors leading to drought in Southwest China. The shortage of precipitation will lead to insufficient water supply, which directly hinders the growth of vegetation and human production, and excessive evaporation will further accelerate the loss of available water resources and destroy the balance of the water cycle [48]. There are only three years in Southwest China with an average annual precipitation of less than 1000 mm since 1966, namely, 2006 (997.14 mm), 2009 (954.97 mm), and 2011 (910.28 mm), and TSDI has detected moderate or severe drought in these three years. Similarly, the evaporation in 2006, 2009, and 2011 is higher than the average value from 2000 to 2016 by 61.73 mm, 48.17 mm, and 22.38 mm, respectively. In addition, the STL method is used to extract the annual trends of TSDI and climate indicators, and the relationship between them is reflected by comparison. As shown in Figure 10, the abnormal decrease of precipitation and the increase of evaporation often correspond to the trough of TSDI, which further confirms that the accuracy of PLSR model simulation results. and TSDI has detected moderate or severe drought in these three years. Similarly, the evaporation in 2006, 2009, and 2011 is higher than the average value from 2000 to 2016 by 61.73 mm, 48.17 mm, and 22.38 mm, respectively. In addition, the STL method is used to extract the annual trends of TSDI and climate indicators, and the relationship between them is reflected by comparison. As shown in Figure 10, the abnormal decrease of precipitation and the increase of evaporation often correspond to the trough of TSDI, which further confirms that the accuracy of PLSR model simulation results.

The Sources of Uncertainties
The uncertainties of TSDI includes grace data and the development of TSDI. In this study, mascon solutions are used as the input parameter to calculate the drought index. Compared with the spherical harmonic solutions, the mascon solutions reduce the additional post-processing and obtain more accurate results. However, there remain multiple uncertainties in using GRACE mascon products to evaluate drought. GRACE data is recommended to assess water storage change information in large basins (> 200000 km 2 ), which limits its application in small scale areas [9]. Besides, long time series (at least 30 years) of GRACE data are needed to determine the baseline of the occurrence and severity of water storage deficit [44]. Although GRACE data can provide vertical integrated TWSA, it is difficult to reflect the contribution of various components of terrestrial water storage (i.e., soil moisture and groundwater) to drought. The GRACE-derived TWSA is an abnormal value of the monthly terrestrial water storage minus the average value from 2004 to 2009, and TWSA may have similar values in areas with abundant or scarce water storage, which cannot reflect the real situation of water availability. The root mean square (RMS) of residual components obtained by STL decomposition is used to evaluate the uncertainty of CSR RL06 mascon products in Southwest China ( Figure 11). The results indicate that the uncertainty in the northern region is lower than that in the southern part, which is consistent with Scanlon et al. [25], who evaluated the RMS of the CSR RL05 mascon products in the global land area.
average value from 2004 to 2009, and TWSA may have similar values in areas with abundant or scarce water storage, which cannot reflect the real situation of water availability. The root mean square (RMS) of residual components obtained by STL decomposition is used to evaluate the uncertainty of CSR RL06 mascon products in Southwest China (Figure 11). The results indicate that the uncertainty in the northern region is lower than that in the southern part, which is consistent with Scanlon et al. [25], who evaluated the RMS of the CSR RL05 mascon products in the global land area. Figure 11. RMS of the residuals in TWSA after the long-term trend and seasonal amplitudes have been removed.
The uncertainty of GRACE data is one of the sources of error, and the developed TSDI also introduces uncertainty. TSDI takes the negative value of TSD for three Figure 11. RMS of the residuals in TWSA after the long-term trend and seasonal amplitudes have been removed.
The uncertainty of GRACE data is one of the sources of error, and the developed TSDI also introduces uncertainty. TSDI takes the negative value of TSD for three consecutive months or more as the standard for drought detection, which may not be suitable for monitoring all types of droughts, especially short-term drought. The definition of drought events also introduces errors, for example, some droughts are just continuing the last drought, especially for the close two drought events, but they are divided into two droughts [34]. Additionally, TSDI only takes GRACE-derived TWSA as the input parameter, ignoring the impact of human activities (i.e., irrigation, livestock, groundwater abstraction) on drought [49].

Future Direction
With the operation of GRACE Follow-On satellites, GRACE time series data will continue to expand, which will provide sufficient data for the study of drought. However, there is a long period (2017.07-2018.05) of missing data between GRACE and GRACE-FO, which hinders GRACE-FO data from being used for drought monitoring. At present, there have been some researches on GRACE data reconstruction using machine/deep learning, which can effectively fill in long-term missing data in large basins [50]. The reconstruction of GRACE records with long time series is very important for drought monitoring and evaluation. However, drought evaluation requires gridded GRACE data (e.g., 0.25 degree), and gridded data reconstruction still has great uncertainty compare with basin scale [39]. Moreover, data assimilation of GRACE data with terrestrial hydrological models will provide more accurate data on changes in different water storage composition, which will contribute to the study of the causes and development of drought. Improving the definition of drought and incorporating human activity into TSDI would also provide more reliable drought assessment results.

Conclusions
GRACE gravity satellites have proven to be an alternative method of drought monitoring, especially in areas where meteorological and hydrological data from station-observed are insufficient. In this study, the GRACE-based TSDI is used to evaluate the drought events in Southwest China from 2003 to 2016, and compare with SPI, SPEI, and SC-PDSI in temporal evolution and spatial distribution. Subsequently, the links between TSDI and climate indicators are clarified using the PLSR model, and the main factors affecting drought are determined. The primary conclusions are as follows: The spatial distribution of TSDI was more consistent with the government report than SC-PDSI. (3) The PLSR model can reveal the links between drought and climate indicators. The VIP results based on the PLSR model indicate that insufficient precipitation has the most significant impact on drought in Southwest China, followed by excessive evaporation. There is a significant positive correlation between precipitation and TSDI (correlation coefficient = 0.78), while a significant negative correlation between evaporation and TSDI (correlation coefficient = -0.64), which further indicates that the decrease of precipitation and excessive evaporation are the causes of drought. In addition, the change trend comparison of precipitation and evaporation with TSDI also verified the results of the PLSR model.