Monthly and Seasonal Drought Characterization Using GRACE-Based Groundwater Drought Index and Its Link to Teleconnections across South Indian River Basins

: Traditional drought monitoring is based on observed data from both meteorological and hydrological stations. Due to the scarcity of station observation data, it is difﬁcult to obtain accurate drought distribution characteristics, and also tedious to replicate the large-scale information of drought. Thus, Gravity Recovery and Climate Experiment (GRACE) data are utilized in monitoring and characterizing regional droughts where ground station data is limited. In this study, we analyzed and assessed the drought characteristics utilizing the GRACE Groundwater Drought Index (GGDI) over four major river basins in India during the period of 2003–2016. The spatial distribution, temporal evolution of drought, and trend characteristics were analyzed using GGDI. Then, the relationship between GGDI and climate factors were evaluated by the method of wavelet coherence. The results indicate the following points: GRACE’s quantitative results were consistent and robust for drought assessment; out of the four basins, severe drought was noticed in the Cauvery river basin between 2012 and 2015, with severity of − 27 and duration of 42 months; other than Godavari river basin, the remaining three basins displayed signiﬁcant negative trends at monthly and seasonal scales; the wavelet coherence method revealed that climate factors had a substantial effect on GGDI, and the impact of Southern Oscillation Index (SOI) on drought was signiﬁcantly high, followed by Sea Surface Temperature (SST) Index (namely, NINO3.4) and Multivariate El Niño–Southern Oscillation Index (MEI) in all the basins. This study provides reliable and robust quantitative result of GRACE water storage variations that shares new insights for further drought investigation.


Introduction
Drought is a dynamic natural phenomenon with high frequency and long duration characteristics that impact ecosystems and society in many ways [1][2][3]. Drought is a common natural disaster that has serious influence on water resources, agriculture, and socio-economic development due to its long-term persistence and frequent occurrence [4][5][6]. Thus, effective evaluation and monitoring of droughts are extremely necessary. Monitoring of traditional drought depends on ground data from meteorological and hydrological stations. Due to the scarcity of station observation data and the spatial heterogeneity of the regional environment, it is difficult to obtain accurate drought distribution characteristics, and also tough to replicate the large-scale information of drought [7].
To precisely monitor and assess drought characteristics, researchers have developed various drought indices such as the Palmer Drought Severity Index (PDSI) [8], Standardized Precipitation Index (SPI) [9], Standardized Precipitation Evapotranspiration Index (SPEI) [10], and Effective Drought Index (EDI) [11], which can be applied to explore drought mand contributes to the overexploitation of surface water and groundwater during drought times. Successful river basin scale drought monitoring is important for water resource management and drought mitigation plans. Several studies have shown that anthropogenic activities and climate change have exacerbated drought-related calamities [13,55,56].
To the best of our knowledge, previous studies have focused on the relationship of several atmospheric variables such as precipitation, temperature, vapor pressure, and humidity with teleconnections in India [57][58][59][60]. Nonetheless, a comprehensive and systematic analysis between GRACE and teleconnections is vague, particularly for India. Therefore, this novel study addresses this research gap by exploring the drought situation over south Indian river basins with GGDI and identifying the linkages between drought and large-scale climate oscillations during 2003-2016. In the present study, we assessed the effect of four major climate oscillations, namely, Multivariate ENSO Index (MEI), Southern Oscillation Index (SOI), Dipole Mode Index (DMI), and NINO3.4 Sea Surface Temperature (SST) on GGDI over Indian river basins using the GRACE TWS dataset for the period of 2003-2016. The detailed analysis was accomplished over Godavari (GRB), Krishna (KRB), Pennar and east flowing rivers between Pennar and Cauvery (PCRB), and Cauvery (CRB) river basins in order to illustrate the linkages between GGDI and climate oscillations. Moreover, gridded monthly and seasonal drought trends were evaluated using the modified Mann-Kendall (MMK) trend test over four river basins (combined). Seasonal trends were evaluated for four seasons of India, namely, (i) post-monsoon rabi (PMon-R -January to March), (ii) pre-monsoon (PMon-April to June), (iii) monsoon (Mon-July to September), and (iv) post-monsoon kharif (PMon-K-October to December).
The objectives of this study were (1) to determine the changing characteristics of terrestrial water storage anomaly (TWSA) over South Indian river basins at monthly, seasonal, and annual timescales; (2) to analyze the spatial distribution and temporal evolution of drought using GGDI; (3) to evaluate the trend characteristics of GGDI over South Indian river basins during 2003-2015 with the MMK trend test; and (4) apply wavelet coherence method to evaluate the relationships between GGDI and climate oscillations.

Study Area
India is the seventh largest country in the world, covering 22 major river basins of which there are 4 river basins, namely, (i) Godavari river basin (GRB), (ii) Krishna river basin (KRB), (iii) Pennar and east flowing rivers between Pennar and Cauvery (PCRB), and (iv) Cauvery river basin (CRB), which were all considered for this study (India-WRIS, 2012). In Figure 1, the study area map is presented; further details regarding the study area can be found in Kumar et al. [33].

Gravity Recovery and Climate Experiment (GRACE)
GRACE monthly products are mainly generated by the Jet Propulsion Laboratory (JPL), Center for Space Research at the University of Texas (CSR), and the German Research Center for Geosciences (GFZ). In the present study, GRACE monthly mass grids (RL06 mascon solutions) processed at Jet Propulsion Laboratory (JPL) with a spatial resolution of 0.5 ° × 0.5° between 2003 and 2016 were used to estimate the terrestrial water storage anomalies (TWSA) (https://grace.jpl.nasa.gov (accessed on 20 March 2021)). The mascon solutions are associated with the baseline period from January 2004 to December 2009 [33,61]. Moreover, glacial isostatic adjustment (GIA) correction was removed, with no need for smoothing filter and north-south striping [61]. To fill the missing monthly GRACE datasets, we adopted the linear interpolation method, which is effective and is extensively used in handling the missing data [33,62].

Global Land Data Assimilation System (GLDAS)
In the present study, the latest release of GLDAS Noah model, i.e., NOAH10_M 2.1 products with the spatial resolution of 1 ° × 1° consistent with GRACE product was adopted for the period of 2003-2016 and used for the analysis (https://disc.gsfc.nasa.gov/ (accessed on 20 March 2021)). All the datasets were resampled to 1 ° × 1° using the bilinear interpolation method, and further analysis was performed in this study.

Global Land Data Assimilation System (GLDAS)
In the present study, the latest release of GLDAS Noah model, i.e., NOAH10_M 2.1 products with the spatial resolution of 1 • × 1 • consistent with GRACE product was adopted for the period of 2003-2016 and used for the analysis (https://disc.gsfc.nasa.gov/ (accessed on 20 March 2021)). All the datasets were resampled to 1 • × 1 • using the bilinear interpolation method, and further analysis was performed in this study.

Meteorological Data
In the present study, gridded precipitation data from the India Meteorological Department (IMD) was considered from 2003 to 2016 with a spatial resolution of 0.25 • × 0.25 • for the study region [63,64] Terrestrial water storage anomalies (TWSA) were derived from GRACE satellite observations. The groundwater storage anomalies (GWSA) at any time t are calculated by subtracting soil moisture storage anomalies (SMSA) and canopy water storage anomalies (CWSA) from the GRACE-based TWSA.
The soil moisture storage anomaly (SMSA) was calculated for NOAH land surface model using the following equation: wre SMSA t = soil moisture anomaly with respect to time t; SMS t = soil moisture at time t; and SMS 2004−2009 = average soil moisture w.r.t at the baseline period from January 2004 to December 2009, the same as that of GRACE terrestrial water storage. Similarly, the canopy water storage anomaly (CWSA) was calculated for NOAH land surface model using the following equation:

GRACE Groundwater Drought Index (GGDI)
In the present study, dimensionless GGDI was implemented to examine the drought characteristics related to groundwater. Firstly, monthly climatology, C i (climatology for month i, (i = 1, 2, . . . 12)), was calculated as follows: where i represents month (i = 1, 2, . . . 12) and n represents number of years. In the present study, GRACE TWS was considered from 2003 to 2016 with n = 14. The monthly climatology C i was calculated for each month individually using groundwater storage anomalies (GWSA). The effect of seasonality in groundwater storage changes is removed by using the monthly climatology [39]. Secondly, the monthly climatology was subtracted from GWSA to obtain groundwater storage deviation (GWSD), which signifies the net deviation in the volume of groundwater storage on the basis of seasonal variability. Finally, the GWSD was normalized by removing the mean and divided by standard deviation as follows: where x GWSD and S GWSD are the mean and standard deviation of GWSD, respectively. GGDI is the normalized net deviation in groundwater storage volume; the GGDI classification is given in Table 1. For detailed information regarding GGDI, refer to Thomas et al. [39].  (Wang et al., 2020).

Grade Classification GGDI
The run theory approach is used to determine the characteristics of drought, such as severity and duration using GGDI. Drought duration (D) is the period of time where the GGDI remains below the fixed threshold value (threshold value of −0.8). The minimum duration of drought is considered as 1 month, as the drought event is defined at aggregation of monthly time scale. Drought severity (S) is the cumulative values of GGDI within the drought duration.

Standardized Precipitation Evapotranspiration Index (SPEI)
Standardized Precipitation Evapotranspiration Index (SPEI) considers potential evapotranspiration (PET) along with precipitation, utilizing all the advantages of SPI. In this study, SPEI was estimated on a 12-month timescale using the IMD precipitation and temperature datasets from 2003 to 2016. The positive SPEI indicates a wet condition and negative SPEI indicates a dry condition. The index is reliable and flexible with respect to space and time in reproducing water deficiencies; therefore, drought characteristics are well assessed with SPEI at different timescales. The detailed procedure regarding SPEI is provided in the Supplementary Materials section.

Modified Mann-Kendall (MMK) Trend Test
The traditional nonparametric Mann-Kendall test is the most widely applied trend test all over the world. However, the persistence of hydrometeorological dataset will affect the Mann-Kendall test results. Therefore, Hamed and Rao [65] estimated the MMK test where it can remove the autocorrelation, which is consistent and robust in finding the trend of a time series [66]. This study implemented the MMK trend test to evaluate the spatial drought trend characteristics over 4 river basins of south India from 2003 to 2016. The detailed procedure regarding MMK test is provided in the Supplementary Materials section.

Wavelet Coherence
Within the time-frequency space, wavelet coherence can be used to determine the relationship between the 2 time series data by estimating the correlation between them that varies between 0 and 1. In accordance with Torrence and Webster [67] and Grinsted et al. [68], coefficient of wavelet coherence between the 2 sets of time series data can be denoted as follows: where R 2 (s, τ) = coherence coefficient minimum and maximum coherence at 0 and 1, and W xy (s, τ) = cross wavelet transforms between two series. Equation (6) resembles the coefficient of determination equation, and thus the wavelet coherence varies between 0 and 1 [69]. S = smoothing operator represented as given below [70]: Climate 2021, 9, 56 7 of 27 The smoothing along wavelet axis is represented as S scale and S time . In the present study, the wavelet coherence was examined at 5% significance level or at the confidence interval > 95%.     Overall, from the beginning of the 21st century, TWSA showed a downward trend over the four river basins on seasonal and annual scales. We can see from the GGDI time series that droughts have become more frequent across these river basins in recent years. For the drought characterization, dry spells of more than three months were considered for drought event analysis in this study [32]. The identified drought events were denoted as "DE", followed by event order for that particular basin.

Basin-Wise GGDI-Identified Drought Event Analysis
As shown in Figure 3a, with reference to GGDI, four major drought events were detected in GRB, i.      were two major drought events that were noticed (SPEI and GGDI). Variations in the drought duration were observed between GGDI-and SPEI-based droughts. More drought events were observed using GGDI when compared with the SPEI. As each drought index is different in terms of construct and variables involved, differences in characterizing drought events are expected among the indices. Thus, variations in the drought events were observed between GGDI and SPEI. Therefore, GGDI-based drought analysis is important and may offer additional insights in identifying the extreme droughts over the river basins in which 50% of the population depends on agriculture.

Gridded Monthly and Seasonal GGDI-Based Drought Trend Characteristics
Gridded monthly and seasonal drought trends were evaluated using MMK trend test over four river basins (combined) during 2003-2016, which are presented in Figure 5. Figure 6a represents the Z s values of GGDI over four river basins during 2003-2016. " ", " ", and " " denote significant (positive and negative) values at 0.1, 0.05, and 0.01 levels, respectively. Figure 6b represents the Kendall tau values.

Monthly Trends
As shown in Figure 5, GRB exhibited monthly significant positive and negative trends from January to December. From January to July, GRB exhibited positive trends at different significant levels (0.01, 0.05, and 0.1 percentiles). The highest significant positive trends were observed in the month of May at the upper part of the basin, whereas the bottom part of the basin showed a downward trend from January to July. From August to December, significant negative trends at 1, 5, and 10% were observed in the downward region of the basin. The highest significant negative trends were seen during August and September. Highly fluctuating positive and negative trends were observed in GRB compared to the other three basins. As shown in Figure 6a, most of the significant positive Z s values of GGDI were observed in GRB compared to other basins. Out of 12, 4 months (January, March, June, and July) exhibited significant positive values at 0.1 and 0.05 levels. In terms of Figure 6b, strong positive trends were observed in the months of March, June, and July.
As seen in Figure 5, KRB displayed monthly positive and negative (non-significant) trends from January to July. From August to December, significant negative trends were observed at 0.01, 0.05, and 0.1 levels. Most of the negative trends were observed in the month of September followed by August. No significant positive trends were observed in the Krishna basin. As shown in Figure 6a, significant negative Z s values of GGDI were observed during August, September, and October months at 0.1 level. The remaining months showed positive and negative Z s values (no significant). Figure 6b shows that strong negative trends were observed from August to October over KRB.  As we can see from Figure 5, PCRB and CRB demonstrated monthly positive and negative (non-significant) trends from January to July. From August to December, significant negative trends were observed at 0.01, 0.05, and 0.1 levels. A complete downward trend was observed over CRB (PCRB) in the month of November (December). In terms of Figure 6a, we can see that significant negative Z s values of GGDI were observed during August to December months at 0.01, 0.05, and 0.1 levels. The highest significant negative Z s values were observed for the month of November for both PCRB and CRB. Overall, KRB and CRB displayed most of the significant negative trends and Z s values compared to GRB and PCRB, and significant positive trends were observed only in GRB. As seen in Figure 6b, strong negative trends were observed in the month of November, followed by October and December in both CRB and PCRB.

Seasonal Trends
The four seasons considered for the study, namely, PMon-R, PMon, Mon, and PMon-K exhibited significant positive trends over GRB. In KRB, CRB, and PCRB, PMon-R and PMon seasons showed upward and downward (non-significant) trends. In comparison, KRB, CRB, and PCRB showed significant negative trends in the Mon and PMon-K seasons. As shown in Figure 5, most of the significant negative trends were observed in the PMon-K, followed by the Mon season. Overall, from season to season, significant positive trends were converted to significant negative trends, with highest significant negative trends shown in the PMon-K season. As seen in Figure 6a, GRB exhibited positive Z s values in all the seasons, with all remaining basins having positive and negative Z s values (significant and insignificant). PMon-K season displayed the highest significant negative Z s values in all the basins, followed by Mon season. In terms of Figure 6b, we see that PMon-R and PMon exhibited strong positive trends in GRB followed by CRB. On the other hand, in PMon-K season, strong negative trends were observed in PCRB and CRB. Decrease in precipitation was observed during 2002-2016, which led to nearly four major drought events. Figure 7 represents the monthly and seasonal precipitation trends using the MMK trend test over the four river basins (combined) during the period of 2003-2016. From Figure 7, we can see that GRB exhibited monthly significant downward trends in the lower and middle parts of the basin during January to December, except in March and June. However, the upper part of the basin showed a significant upward trend, except in May and August. In January, and from August to December, GRB exhibited significant downward trends (0.05, and 0.1 percentiles, respectively). The highest significant positive trends were observed in the months of June and July at 0.01 and 0.05 percentiles, respectively, whereas downward trends were observed during November, December, and January at 1, 5, and 10 %, respectively. The highest significant negative trends were shown during August and September. Highly fluctuating positive and negative trends were observed in GRB compared to other three basins. From Figure 6a From Figure 7, we see that monthly significant positive trends were observed in June, August, and December over PCRB and in the months of June and December over CRB at 0.01 and 0.05 levels, respectively. Highly significant negative trends were observed during October in PCRB and CRB at 0.01 and 0.05 levels. Significant positive and negative trends were observed during January to May and September to November in PCRB and CRB at 0.01, 0.05, and 0.1 levels.

Seasonal Trends
Seasonally (Figure 7), significant positive trends were observed during the PMon and Mon seasons at 0.01 and 0.05 levels in GRB and KRB. Significant negative and positive trends were observed during PMon-R and PMon-K over GRB, KRB, CRB, and PCRB at 0.01, 0.05, and 0.1 levels. However, in case of PMon and Mon seasons, PCRB and CRB exhibited the most significant negative trends compared to positive trends at 0.05 and 0.01 levels. Highly significant positive trends were observed in CRB compared to all other basins.
Overall, from the analysis, GGDI was found to be strongly influenced by variability of precipitation in the study region. Results stated that study region experienced significant decreasing trend in precipitation and GGDI. Assessment of GGDI and precipitation variability showed a significant linear trend both monthly and seasonally.

The Correlation between GGDI and Climate Factors
Previous studies have shown that droughts were closely related to climate variables [47,48,55,60]. In the present study, MEI, NINO3.4, SOI, and DMI were chosen to describe the influences of teleconnections over droughts. Moreover, wavelet coherence was employed to evaluate the link between GGDI and climate factors over South Indian river basins during 2003-2016.
The wavelet coherence between monthly GGDI and climate factors were presented in For KRB, the wavelet coherence between GGDI and climate factors are provided in Figure 9. It can be noted from the figure that high wavelet coherence is noticed at an annual scale, characterizing the dominant effect of groundwater for MEI, SOI, NINO 3.4, and DMI. Interannual variability was detected at time scales of 2 to 14 months for MEI, 4-10 months for SOI, and 2-10 months for NINO3.4. However, for DMI, interannual variability highly varied in the overall period.
For CRB, the wavelet coherence between monthly GGDI and climate factors are presented in Figure 10. From Figure 10     For PCRB, the wavelet coherence between monthly GGDI and climate factors were presented in Figure 11. As seen in Figure 11, the annual variability of MEI was significantly dominant for all the years between the time scale of 16-32 months. Interannual variability was also observed between the time scale of 4-10 months at different years. The influence of interannual variability was observed at different years between the time scale of 2-14 months. However, annual variability was seen for all the years between 16 and 32 months. Annual variability was comparatively less in NINO3. 4

Influence Factors of Drought
Droughts in India pose extraordinary challenges to the food production, socio-economic aspects, livelihood, and gross domestic product. India has a long history of droughts with lasting effects on crops, surface, and subsurface water resources, and rural livelihoods [55,71]. Water availability and crop production were affected by the recent drought of 2015-2016 over large parts of southern India with lower reservoir storage. Additionally, the 2015-2016 drought affected around 330 million people and caused groundwater depletion in the South Indian states [72,73]. Failure of monsoon rainfall or its receipt in smaller quantities may often result in drought over major parts of India. The effect of anthropogenic factors on drought, along with natural factors, cannot be overlooked. Due to the uneven distribution of rainfall, spatially and temporally, surface and subsurface water resources are scarce over India [33]. Mishra [73] stated that the 2015-2018 drought affected groundwater and surface water availability in the southern part of India and was linked to climate indices. Farmers disproportionately use electricity and fossil fuels to pump groundwater to compensate for the lack of rainfall. In particular, cultivation costs for rice and other rainy season (kharif) crops are also rising due to increased use of energy and diesel for the collection of groundwater [74]. Excessive withdrawal of groundwater to save crops in drought conditions has drained groundwater in most parts of the world, therefore ultimately triggering a drought crisis.
Researchers have established that climate factors play a major role in the process of drought formation [42,75]. Additionally, the wavelet coherence results from Section 3.5 have shown that climate factors (MEI, SOI, DMI, and NINO3.4) have had an extreme influence on drought evolution. In particular, in the Indian regions, MEI, SOI, and NINO3.4 had the greatest influence on drought (Figures 8-11). There are several teleconnections that influence the variability of TWS along with its components over India. Although earlier studies primarily focused on the effect of ENSO over TWS [46,47], it is unlikely to be able find a single indicator to represent all climatic variability features over large regions [76]. In the current study, four widely used climate factors were considered, and their links with GGDI were evaluated; the results show that for each river basin, the teleconnections differed considerably with GGDI. Therefore, TWS attributions and predictions or indices calculated using TWS (GGDI) centered on a single teleconnection should be treated with caution, and multiple teleconnections are suggested for the assessment of TWS changes and their components.
As shown in Figures 8-11, the coherence of GGDI with teleconnections (MEI, SOI, DMI, and NINO3.4) at ≈32 months period may be due to the correlation between climate indices (correlation o MEI with SOI/DMI/NINO3.4). Therefore, analyzing the standalone effect of teleconnection factors on the GGDI series may provide better correlation between GGDI and teleconnection after removing the effect of other influential time series [77,78].

Uncertainty Analysis
Numerous areas of uncertainty were encountered in the present study. First, to avoid the uncertainty induced by observational and data processing errors, we considered modern mascon solution data instead of the spherical harmonic coefficients data [79]. Still, the mascon solutions developed by various organizations still have ambiguities due to diverse background models as well as data processing approaches. In addition, the JPL mascon solutions considered different hydrological models to adjust the scale factors, which eventually contributes to the presence of uncertainty. Second, to minimize the GLDAS uncertainty, the ensemble mean of several hydrological models (Noah, VIC, CLM, Mosaic) was suggested [37]. However, similar to GRACE data, the Noah model has the same spatial resolution. Thus, we adopted the GLDAS Noah model outputs in the present study to reduce the uncertainty related to spatial resolution in evaluating water storage calculations. Finally, using linear interpolation techniques, the missing GRACE data were filled out, which may have triggered some uncertainties. However, since the approach is prevalent, easy, and extensively used in the handling of missing data, we followed similar techniques in filling the GRACE dataset gaps [33,80], and the results suggest that the linear interpolation approach was appropriate.

Advantages and Limitations
GRACE satellite gravimetry plays an important role in the identification of drought in regions where data related to water storage variations is inadequate [81]. To reduce the influence of various errors in the GRACE spherical harmonic coefficients, one can apply various filtering processes, yet results suggest the possibility of weak signals in the derived product. Therefore, scale factors are applied to recover the signal leakage caused by the filtering processes [82]. In order to resolve these data processing errors, we adopted GRACE mascon solutions, which are equal to or superior to traditional GRACE spherical harmonic coefficients, in the current research [65]. Therefore, using GRACE mascon solutions, changing characteristics of TWS, in identifying the teleconnections with GGDI, and hence the drought situation over river basins in southern India were explored. In regions where hydrometeorological data are minimal, GRACE data are an important means of estimating and managing the drought. The GGDI drought index was evaluated in this analysis to identify drought events. The GGDI is a normalized index that can be used to objectively compare spatio-temporal drought, providing a strong evidence for evaluating surplus and deficit groundwater availability [39]. This study positively established the drought events between 2003 and 2005, 2008 and 2010, and 2013 and 2015, which were consistent with the results of Sinha et al. [17] and Kumar et al. [33]. Severe drought events between 2008 and 2010 for GRB, 2015 and 2016 for KRB, and 2003 and 2004 for PCRB were also reported by Kumar et al. [33], which are consistent with the drought events obtained using GGDI in this study ( Figure 3, Table 2). Moreover, the interaction of GGDI with climate indices exhibited the fact that teleconnections had a substantial effect on drought across southern India's river basins.
Although GGDI based on GRACE dataset can be used to categorize drought characteristics effectively and expediently, there are still some limitations. GRACE datasets have only been available since 2002, covering only a few years of data; with longer temporal datasets, the results will therefore be more accurate [32]. Although GRACE data resolution is relatively low, it considers the changes in water storage (including soil water and surface snow) that comes from rainfall, evapotranspiration, river transportation, and deep underground infiltration. GRACE is the key technology in gravity satellite sensors to improve accuracy and monitoring gravity field in terrestrial hydrology. GRACE provides realistic spatiotemporal variations of vertically integrated measurement of water storage at precision of tens of millimeters of equivalent thickness. The new GRACE Follow On (GRACE FO) satellite datasets will provide a valuable solution for the long-term evaluation of TWS variations and their associated studies, resulting in significant improvements in our knowledge of GRACE-related studies. Moreover, the effect of anthropogenic activities (groundwater extraction, regional water division, mining) over mass changes of the earth surface cannot be overlooked [83]. These effects (influence of human activities) are generally ignored due to the lack of observed datasets, as well as difficulties in the collection and measurement of relevant information. Therefore, analyzing the influence of anthropogenic activities on variations in water storage may provide a fresh insight into the future of science.

Conclusions
In the present study, during 2003-2016, the drought characteristics were examined and evaluated using the GRACE-based groundwater storage as a metric over four major river basins in India. The spatial distribution, temporal evolution of drought, and trend characteristics were analyzed using GGDI during 2003-2016 over KRB, CRB, GRB, and PCRB. Then, using the wavelet coherence method, we evaluated the relationship between GGDI and climate factors. GRACE datasets provide significant benefits in detecting droughts and revealing information about large-scale groundwater depletion, where hydrometeorologi-cal data are limited and data related to water storage variations are insufficient. This study has provided reliable and robust quantitative results of GRACE water storage variations that provide a new approach to link surface and subsurface condition while investigating droughts, and the methodology can be applied to any other region. The key findings from this research are as follows: The monthly and seasonal trends were evaluated using the MMK test. Significant monthly negative trends were observed during August to December in KRB, CRB, and PCRB. Seasonal negative trends were also significant in Mon and PMon-K in CRB, KRB, and PCRB.

•
The wavelet coherence analysis effectively demonstrated the teleconnections between climate indices and drought events. The influence of SOI on drought was significantly high, followed by NINO3. 4 and MEI in all the basins. SOI has the strongest impact in detecting the progression of drought compared to other climate indices in these river basins. Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. We acknowledge the partial support received by the corresponding author from the Virginia Agricultural Experiment Station (Blacksburg) and the Hatch Program of the National Institute of Food and Agriculture, the U.S. Department of Agriculture (Washington, DC, USA).