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

Abstract

Traditional drought monitoring is based on observed data from both meteorological and hydrological stations. Due to the scarcity of station observation data, it is difficult 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 significant 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 significantly 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.

1. 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 characteristics. PDSI, SPI, and SPEI have been the most widely used drought indices in recent decades [12,13]. Despite their broad applicability, these indices have their drawbacks and may not reproduce the conditions of drought accurately [14,15].

The aforementioned drought indices are related to land surface conditions and only a few indices are interrelated with remote sensing methods. Over the past years, the capabilities of remote sensing products amplified many folds, assessing the spatial and temporal variations at local and global scales. Gravity Recovery and Climate Experiment (GRACE) satellite measures the earth’s gravity field to evaluate water storage changes globally [16,17]. The GRACE terrestrial water storage (TWS) dataset comprises all forms of water stored above and below the land surface that constitutes groundwater, surface water, soil moisture, and snow water equivalent [18,19]. GRACE-based TWS variations are commonly used in the evaluation of different hydrological assessments such as groundwater storage changes [20,21,22,23,24,25,26], terrestrial water storage variations [27,28,29,30], and hydrologic drought characterization [17,18,31,32,33].

Recently, many studies have been focused on monitoring and characterizing regional droughts using GRACE TWS [34,35,36]. Thomas et al. [32] developed a framework for assessing drought characterization, focusing on the total water storage deficits using GRACE TWS. Cao et al. [37] utilized GRACE TWS to develop a total storage deficit index for China. The Hydrological Drought Index was introduced by Yi and Wen [38] using GRACE data on the continental United States. The water storage deficit index for drought quantification was developed by Sinha et al. [18]. Zhao et al. [36] proposed a GRACE-based global gridded drought severity index named GRACE-DSI. Thomas et al. [39] evaluated the GRACE Groundwater Drought Index (GGDI) over the central valley of California, a regional aquifer subjected to significant drought periods and intensive human activities. Sinha et al. [17], with GRACE TWS and rainfall, developed the combined climatologic deviation index over India.

Due to large-scale climate variations in India, the spatio-temporal availability of surface and groundwater is very diverse and affects the agricultural and industrial productivity of the country [40]. A 2016 drought in India affected 330 million people with a more than USD 100 billion loss in the economy [41]. From this perspective, for the conservation of water resources, it is crucial to understand the variations of surface water and groundwater and its association with the teleconnection in India.

Utilizing the GRACE data, many studies have investigated drought characteristics throughout the globe. These studies only verified the capabilities of drought using GRACE data but not the associations between GRACE-based droughts and teleconnection factors. It is clear from earlier studies that telecommunication factors have a major effect on drought [42,43]. Many worldwide attempts have been made over past years to establish the relationship between climate variability and GRACE TWS changes, with most studies focused on El Niño–Southern Oscillation (ENSO). To evaluate the associations between the Multivariate ENSO Index (MEI) and GRACE mass anomalies, Phillips et al. [44] utilized the monthly GRACE TWS. Huang et al. [45] concluded that hydrological drought over the Columbia River basin was greatly influenced by ENSO and Arctic Oscillation (AO). Over the entire globe, Ni et al. [46] examined the links between ENSO and GRACE TWS. Vissa et al. [47] evaluated the relationship between ENSO-induced groundwater changes derived from GRACE and GLDAS over India. Liu et al. [48] explored the role of teleconnections over TWS variations within the Asian and eastern European regions. In the same way, few studies related to linkages between GRACE and teleconnections include Ni et al. [46], Chen et al. [49], Zhang et al. [50], Ndehedehe et al. [51], Anyah et al. [52], Han et al. [53], and Wang et al. [54]. Thus, climate variables influence the drought directly or indirectly that results in detailed investigation between them.

It is worth mentioning that, in India, approximately 50% of the population depends on agriculture that relies on surface and groundwater resources. As the availability of surface water is not uniform through space and time, groundwater resources have emerged as the primary source for agriculture, domestic application, and industry. Increased water demand 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.

2. Materials and Methods

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

2.2. Data

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

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

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

2.2.4. Climate Data

Monthly climate oscillations, namely, NINO 3.4, Multivariate ENSO Index (MEI), Southern Oscillation Index (SOI), and Dipole Mode Index (DMI) for the period of 2003–2016 were utilized in the study. Monthly SST anomaly data (NINO3.4) were obtained from http://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Data/nino34.long.anom.data (accessed on 20 March 2021). For ENSO, the MEI was selected and obtained from https://www.esrl.noaa.gov/psd/enso/mei (accessed on 20 March 2021). The SOI data were obtained from the NOAA Earth System Research Laboratory (https://psl.noaa.gov/gcos_wgsp/Timeseries/SOI/ (accessed on 20 March 2021)). The Indian Ocean Dipole (IOD) is measured as DMI due to the dipole mode in the tropical Indian Ocean, and the DMI data were obtained from http://www.jamstec.go.jp/frcgc/research/d1/iod/DATA/dmi (accessed on 20 March 2021).

2.3. Method

2.3.1. Retrieval of Groundwater Storage Change

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:

$$SMS{A}_t=SM{S}_t-\overline{SM{S}_{2004-2009}}$$

(1)

wre $SMS{A}_t=$ soil moisture anomaly with respect to time $t$; $SM{S}_t=$ soil moisture at time $t$; and $\overline{SM{S}_{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:

$$CWS{A}_t=CW{S}_t-\overline{CW{S}_{2004-2009}}$$

(2)

wre $CWS{A}_t=$ canopy water storage anomaly w.r.t time t; $CW{S}_t=$ canopy water storage at time t; and $\overline{CW{S}_{2004-2009}}=$ average canopy water storage w.r.t the base line period from January 2004 to December 2009, same as that of GRACE terrestrial water storage.

The GWSA is calculated as

$$GWS{A}_t=TWS{A}_t-SMS{A}_t-CWS{A}_t$$

(3)

2.3.2. 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,\dots 12)$), was calculated as follows:

$$C_i=\frac{1}{n_i}{{\displaystyle \sum }}_1^{n_i}GWS{A}_i$$

(4)

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:

$$GGDI=\frac{GWS{D}_t-{\overline{x}}_{GWSD}}{S_{GWSD}}$$

(5)

where ${\overline{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. Classification of the Gravity Recovery and Climate Experiment (GRACE) Groundwater Drought Index (GGDI) used in the study (Wang et al., 2020).

Grade	Classification	GGDI
I	No drought	−0.5 < GGDI
II	Mild drought	−1.0 < GGDI ≤ −0.5
III	Moderate drought	−1.5 < GGDI ≤ −1.0
IV	Severe drought	−2.0 < GGDI ≤ −1.5
V	Extreme drought	GGDI ≤ −2.0

. For detailed information regarding GGDI, refer to Thomas et al. [39].

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.

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

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

2.3.5. 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:

$$R^2\left(s,\tau \right)=\frac{{\left|S\left(s^{-1}{W}_{xy}\left(s,\tau \right)\right)\right|}^2}{S\left(s^{-1}{\left|W_x\left(s,\tau \right)\right|}^2\right).S\left(s^{-1}{\left|W_y\left(s,\tau \right)\right|}^2\right)}$$

(6)

where $R^2\left(s,\tau \right)=$ coherence coefficient minimum and maximum coherence at 0 and 1, and $W_{xy}\left(s,\tau \right)=$ 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]:

$$S\left(W\right)=S_{scale}\left(S_{time}\left(W\left(s,\tau \right)\right)\right)$$

(7)

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

3. Results

3.1. Changing Characteristics of TWSA over River Basins

Seasonal and annual scale GRACE TWSA analysis was performed over four basins. The results showed distinct seasonal and annual TWSA characteristics over these river basins. In the Mon and PMon-K seasons, the GRB showed positive TWSA values, while negative TWSA values were detected in the PMon season. The positive and negative TWSA values were shown by the PMon-R season. On annual basis, TWSA displayed a significant upward trend with the highest values observed in 2013 and 2014. The TWSA in the GRB tended to decrease from 2003 to 2005, leading to a significant drought between 2003 and 2005, as shown in Figure 2. As seen in Figure 2, except for the year of 2015, the KRB showed positive TWSA values in the Mon and PMon-K seasons. Except for the PMon-R season in 2011, the PMon and PMon-R seasons displayed negative TWSA values. Annual TWSA values showed a significant upward and a downward trend of TWSA values and most of the negative trends were observed between 2003 and 2005, and between 2012 and 2016, leading to severe droughts.

From Figure 2, the CRB and PCRB exhibited positive and negative TWSA for PMon-K and PMon seasons, respectively. Positive and negative TWSA values were observed in the Mon and PMon-R seasons, with most of the negative values between 2003 and 2005, and between 2012 and 2016. Annual TWSA values exhibited a downward trend between 2012 and 2016, leading to severe drought, followed by the 2003–2005 event. Overall, from the beginning of the 21st century, TWSA showed a downward trend over the four river basins on seasonal and annual scales.

3.2. Basin-Wise GGDI-Identified Drought Event Analysis

Figure 3 represents GGDI-based temporal evolution characteristics of drought and major drought events for the four river basins from 2003 to 2016. The solid red line represents the GGDI, and the green shaded region indicates the period of drought events. The GGDI exhibited upward and downward trends over each river basin with different change characteristics, demonstrating that the droughts reported using GGDI were increasing in these river basins during the period of 2003–2016. For all the four basins, major decreasing trends were observed between 2003 and 2005, 2011 and 2013, and 2014 and 2016. 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.

As shown in Figure 3a, with reference to GGDI, four major drought events were detected in GRB, i.e., during (i) DE1—January to December 2003, (ii) DE2—March 2004 to July 2005, (iii) DE3—June 2008 to July 2010, and (iv) DE4—August 2015 to February 2016. The longest drought event over GRB extended for 26 months between 2008 and 2010 (DE3). Among four drought events, DE1, 2, and 4 were characterized as moderate drought, and DE3 as severe drought.

As shown in Figure 3b, five major drought events were observed in KRB: (i) DE1—January 2003 to May 2004, (ii) DE2—July 2004 to July 2005, (iii) DE3—June 2012 to May 2013, (iv) DE4—June to November 2014, and (v) DE5—July 2015 to June 2016. Among these drought events, DE1, 3, and 4 were characterized as mild drought; DE2 as moderate drought; and DE5 as extreme drought. For KRB, the longest drought was observed for 17 months between 2003 and 2004 (DE1).

As shown in Figure 3c, four major drought events were identified over CRB: (i) DE1—January 2003 to May 2004, (ii) DE2—July 2004 to August 2005, (iii) DE3—June 2012 to November 2015, and (iv) DE4—August to December 2016. DE1, 3, and 4 were characterized as severe drought, and DE2 as moderate drought, with the longest drought extended for 42 months between 2012 and 2015 (DE3).

As shown in Figure 3d, five major drought events were observed in PCRB: (i) DE1—January 2003 to May 2004, (ii) DE2—July 2004 to August 2005, (iii) DE3—June 2012 to July 2013, (iv) DE4—June 2014 to February 2016, and (v) DE5—August to December 2016. DE1 and 3 were characterized as moderate drought, and DE2, 4, and 5 were characterized as severe drought, with the longest drought extended for 21 months between 2014 and 2016 (DE4). Overall, for the four river basins of South India, GGDI can be considered a strong indicator of drought.

3.3. SPEI-Based Drought Event Analysis

Figure 4 represents SPEI-12 timeseries with major drought events represented in red bands for the four river basins from 2003 to 2016. The SPEI-12 exhibited upward and downward trends over each river basin during the period of 2003–2016. For all the four basins, major decreasing trends were observed between 2003 and 2004, 2012 and 2013, and 2014 and 2016. From the SPEI-12 time series, it is evident that droughts have become more recurrent and prolonged across four river basins in recent years. The identified drought events were denoted as “DE”, followed by event order for that particular basin.

As seen in Figure 4a, four major drought events were evident in GRB: (i) DE1—January 2003 to August 2003, (ii) DE2—February 2009 to September 2010, (iii) DE3—November 2011 to April 2013, and (iv) DE4—April 2014 to December 2016. DE1 and 3 were characterized as moderate drought, and DE2 and 4 were characterized as severe drought, with the longest drought extended for 33 months between 2014 to 2016 (DE4).

As seen in Figure 4b, four major drought events were observed in KRB: (i) DE1—January 2003 to April 2004, (ii) DE2—February 2009 to January 2010, (iii) DE3—November 2011 to June 2013, and (iv) DE4—February 2014 to December 2016. DE1, 2, and 3 were characterized as moderate drought, and DE4 was characterized as severe drought, with the longest drought extended for 35 months between 2014 to 2016 (DE4).

As seen in Figure 4c, four major drought events were observed in CRB: (i) DE1—January 2003 to May 2004, (ii) DE2—June 2012 to October 2013, (iii) DE3—January 2014 to December 2015, and (iv) DE4—May 2016 to December 2016. DE1, 2, and 4 were characterized as moderate drought, and DE3 was characterized as severe drought, with the longest drought extended for 24 months between 2014 to 2015 (DE3).

As seen in Figure 4d, four major drought events were observed in PCRB: (i) DE1—January 2003 to April 2004, (ii) DE2—February 2009 to May 2010, (iii) DE3—November 2011 to September 2013, and (iv) DE4—January 2014 to December 2015. DE1, 2, and 3 were characterized as moderate drought, and DE4 was characterized as severe drought, with the longest drought extended for 24 months between 2014 to 2015 (DE4).

Comparisons of GRACE-based drought with SPEI from Figure 3 and Figure 4 were analyzed. In GRB, drought events were commonly observed during 2008–2010 and 2015–2016 for GGDI as well as SPEI. In KRB, three major drought events were observed between 2003 and 2004, 2012 and 2013, and 2015 and 2016 for both SPEI and GGDI. In CRB, for both SPEI and GGDI, two major drought events were observed during 2003–2005 and 2012–2015. IN PCRB, 2003–2004 and 2012–2016 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.

3.4. Basin-Wise Drought Characteristics Using GGDI

Table 2. Drought characteristics (severity and duration) calculated from GGDI for the drought events identified for GRB, KRB, CRB, and PCRB.

Time Period	Severity	Duration (No. of Months)
Godavari (GRB)
Jan 2003 to Dec 2003	−8.72	12
Mar 2004 to Jul 2005	−12.91	17
Jun 2008 to Jul 2010	−14.64	26
Aug 2015 to Feb 2016	−4.87	7
Krishna (KRB)
Jan 2003 to May 2004	−11.56	17
Jul 2004 to Jul 2005	−10.08	13
Jun 2012 to May 2013	−8.47	12
Jun to Nov 2014	−1.56	6
Jul 2015 to Jun 2016	−15.72	12
Aug to Dec 2016	−5.41	5
Cauvery (CRB)
Jan 2003 to May 2004	−19.27	17
Jul 2004 to Aug 2005	−10.56	14
Jun 2012 to Nov 2015	−27.02	42
Aug to Dec 2016	−6.71	5
Pennar and east flowing rivers between Pennar and Cauvery (PCRB)
Jan 2003 to May 2004	−16.33	17
Jul 2004 to Aug 2005	−12.52	14
Jun 2012 to Jul 2013	−9.71	14
Jun 2014 to Feb 2016	−13.38	21
Aug to Dec 2016	−7.14	5

represents the drought characteristics (severity and duration) calculated from GGDI over four river basins. For GRB, the highest severity of −14.64 was observed for a duration of 26 months between June 2008 and July 2010. In KRB, the highest severity of −15.72 was observed for a duration of 12 months (July 2015 to June 2016), followed by −11.56 severity with a duration of 17 months (January 2003 to May 2004). The CRB experienced the height drought period among all the four basins, with severity of −27.02 observed for 42 months from June 2012 to November 2015. PCRB observed the highest severity of −16.33 for a duration of 17 months (January 2003 to May 2004), followed by severity of −13.38 with the highest duration of 21 months (June 2014 to February 2016) in this basin. All the four basins experienced droughts during 2003–2005 and 2015–2016.

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

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

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

3.6. Gridded Monthly and Seasonal Trend Characteristics in Terms of Precipitation

3.6.1. Monthly Trends

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, we can see that most of the significant positive $Z_s$ values of GGDI were observed in GRB in comparison to other basins. Out of the 12 months, 3 months (March, April, and June) exhibited significant positive values at 0.01 and 0.05 levels.

From Figure 7, we see that during February to April and December, KRB displayed monthly significant positive trends at 0.01, 0.05, and 0.1 levels in the upper part of the basin, whereas the other portion of the basin showed significant negative trends at the 0.05 level. The basin was covered with significant negative trends in the months of May, September, October, and November at the 0.1 level. No significant positive trends were observed in KRB. As shown in Figure 6a, significant negative $Z_s$ values of GGDI were observed during the months of August, September, and October at the 0.1 level. Both insignificant negative and positive trends were negotiable in KRB.

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.

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

3.7. 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 Figure 8 over GRB. From Figure 8, the coherence at interannual variability was observed continuously between 2007 and 2016 on time scales of 20–36 months, and intermittency was observed between 4 and 12 months at different years for MEI. In the case of SOI, coherence at interannual variability was observed between 2002 and 2005, and 2007 and 2016 periods on time scales varying from 16 to 40 months. For NINO3.4, intermittency was reduced between the 2002 and 2007 period for time scales of 10–34 months, and intermittency was observed at different years between 4 and 10 months. Compared to other teleconnections at different scales during different years, the effect of DMI was weak.

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, we can see that the influence of MEI is observable over the time scales of 10–30 months for the period 2002–2011. Interannual variability of SOI was observed within the 10 to 12-month time scale, whereas annual variability was observed for the time scale of 16–32 months over the period of 2002–2012. Annual variability of NINO3.4 was observed between 2002 and 2011 for an 18–32-month time scale, and interannual variability was observed for the time scale of 6–16 months. Influence of DMI was highly significant for the time scale of 14–40 months; for the period 2002–2009, the annual variability was significantly high compared to NINO3.4, SOI, and MEI.

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 when compared with MEI, SOI, and DMI. Interannual variability was observed between the time scales of 4–14 months over the years of 2009–2016. Annual variability of DMI was highly significant for 16–32-month time scales between the years 2002 and 2010, continuing for 33–40-month time scales for all the years. Moreover, interannual variability was observed between the time scale of 4–14 months for different years.

3.8. Spatial Distribution of Drought

Generally, most of the drought events were seen during 2003–2005, 2008–2010, 2012–2013, and 2014–2016. Spatial distribution of the yearly average GGDI timeseries was plotted (Figure 12) to depict the spatial variation of drought throughout the basins. From Figure 12, we can see that 2003, 2015, and 2016 were the most drought-affected years in all the river basins. KRB exhibited the severe drought throughout the basin in 2016 followed by PCRB and CRB. For GRB, the GGDI varied between −0.5 and −2 during 2003–2005, 2009, 2012, 2015, and 2016. For KRB, it varied between −0.5 and −2 for the years 2003–2005, 2009, and 2012–2016. For CRB and PCRB, a range from −0.5 to −2 for the years 2003–2005 and 2012–2016 were noticed. KRB was the severely affected basin during the past decade compared to the other basins. Drought had been considerably relieved during 2006, 2007, and 2011. Therefore, appropriate drought-resistant measures should be implemented in these river basins to reduce the impact of drought disaster and improve the capability of drought resistance.

4. Discussion

4.1. 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 (Figure 8, Figure 9, Figure 10 and Figure 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 Figure 8, Figure 9, Figure 10 and Figure 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].

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

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

5. 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 hydrometeorological 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 distinct seasonal and annual variations of TWSA were observed in four river basins. The PMon and PMon-R seasons exhibited negative TWSA values, while Mon and PMon-K seasons showed positive TWSA variations in all the river basins. Annually, TWSA values showed significant upward and downward trends, with most of the negative trends observed between 2003 and 2005, and 2012 and 2016, indicating severe droughts.
• The GGDI-identified drought events exhibited different temporal change characteristics in all the river basins. The most severe drought event was observed in CRB between 2012 and 2016, followed by GRB between 2008 and 2010. All the four basins exhibited drought events between 2003 and 2005, and KRB, CRB, and PCRB experienced droughts between 2012 and 2016.
• Drought severity and duration were evaluated using GGDI for four river basins. The CRB experienced the longest drought period among all the four basins, with a severity of −27.02 observed for 42 months during June 2012 to November 2015.
• 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.