Unraveling the pertinence of drought indices in the changing climate

The increasing divergence/disagreement between the meteorological drought indicators, majorly due to the changing climate conditions, has disputed their pertinency in the monitoring and modeling of drought events. In the present study, we attempted to quantify the divergence between widely used drought indices such as the standardized precipitation index (SPI), standardized precipitation evapotranspiration index (SPEI), and the newly developed standardized net precipitation index (SNEPI) and traced the evolution of this divergence/disagreement. A persistent presence of enhanced dry extremes in the annual divergence of SPEI with SPI and SNEPI was observed, which raises doubts about the utility of traditional drought indices in capturing the changing climatic characteristics. The seasonal dispersion of this annual divergence revealed a strong, spatiotemporally dominant, signal of dryness in the monsoon season by SNEPI, which SPEI failed to unravel. Furthermore, the attribution analysis revealed that the SPI-SPEI disagreements evolved under the influence of the shortwave radiative fluxes. In addition to that, the divergence between SPEI and SNEPI, is driven by the characteristics of wet spells, with the relationship strengthening in the monsoon season and tropical climate zones. All these findings caution us to evaluate the need to incorporate the changing precipitation characteristics before deploying drought indices in operational uses, especially in tropical climate zones.


Introduction
A realistic assessment of drought has always been a subject of bubbling ambiguity, primarily generated by the presence of a multitude of drought indicators. These indicators impart enhanced uncertainty in drought analysis owing to their mutually exclusive methodologies, changing variable behavior, and the data sources employed These disagreements have been primarily attributed to enormous evaporative losses induced by the increasing global temperatures in the changing climate (Carlson and Benjamin 1980, Estournel et al 1983, Houghton et al 1990, 1992, Arnel and Reynard 1996, Malek et al 2018. However, the climate change impacts propagating to drought behavior are not solely explained by the increased global temperatures, but also by the changing precipitation characteristics. Reportedly, precipitation variance is increasing with no significant change in the overall magnitude (Goswami et al 2006, Trenberth 2008 Xu et al 2018, Giorgi et al 2019. Such precipitation behavior causes long intra-monthly dry spells, catalyzing drought impacts on agriculture and social dynamics, which traditional indices like standardized precipitation evapotranspiration index (SPEI) may overlook. As a result, drought spells during a wet month are mischaracterized (Engelbrecht et al 2006, Udmale et al 2014, Lemma et al 2016, Konapala et al 2017, Ma et al 2017. This misrepresentation of drought can impact its preparedness, especially in regions where the agrarian economy is dependent on seasonal rainfall. For instance, despite the occurrence of monsoon floods in South Asia, the agricultural sector suffered from drought owing to highly skewed distribution of extreme wet spells (Madhusoodhanan and Sreeja 2019). Considering these shortcomings of the contemporary drought indicators, a new index, standardized net precipitation distribution index (SNEPI), was developed (Singh et al 2021). With SNEPI addressing the intraperiod distribution characteristics of net precipitation, an expansion of its inherent divergence with other contemporary indices can be hypothesized, especially in the monsoon season.
This study attempts to analyze the pertinence of drought indicators in monitoring dry extremes over regions experiencing strong seasonality, through analysis of divergence. The authors focus on identifying seasons/climate zones experiencing erratic precipitation behavior, wherein the operational value of SNEPI could be maximum. Further, the principal component analysis (PCA) is used to isolate the various physical processes driving the inter-index divergences in the changing climate. A brief overview of the study area, data and methods is discussed in the next section followed by the results and discussion subsequently.
The study utilizes daily gridded precipitation (0.25 • × 0.25 • ) (Pai et al 2014) data from January 1951 to December 2016, released by IMD. Daily gridded (0.25 • × 0.25 • ) potential evapotranspiration (PET), estimated from the Penman-Monteith method (Monteith 1964, Shuttleworth 1993, downward longwave radiation (R LW ), downward shortwave radiation (R SW ), specific humidity (SH), wind speed (W s ) and air temperature (T air ) data was obtained from version3, global meteorological forcing data set, available at the National Centre for Atmospheric Research using URL https://rda.ucar. edu/datasets/ds314.0/, released by Sheffield et al (2006). It is to be noted that though the sensitivity of PET-estimation method has little/no impact on the final drought index value (Singh et al 2021), employing a physical-based PET-estimation method could provide more realistic identification of various physical processes governing the evolution of the interindex divergence.

Drought indices
Standardized precipitation index (SPI) (McKee et al 1993), and SPEI (Vicente-Serrano et al 2010), utilize long-term monthly precipitation (P) and net precipitation (D = P − PET) data, respectively, to evaluate droughts at multiple timescales (minimum 1 month). The data is fitted to a suitable probability distribution function (PDF). From the fitted PDF, the cumulative distribution function (CDF) is obtained, which is later converted to the standard normal distribution, enabling spatio-temporal intercomparison.
SNEPI (Singh et al 2021) evaluates drought at multiple timescales through the refined monthly aggregate (D m ). D m , at a particular timescale, inherently accounts the distribution of D over excess/ deficit periods, identified by a suitable threshold. The authors recommend the use of D = 0 as a universal threshold, as it acts as a perfect balance between atmospheric demand (PET) and supply (P). The intraperiod magnitude distribution is also considered through the computation of a uniformity coefficient (U c ). D m is fitted to a suitable established/empirical PDF, followed by estimation of the CDF and standardization. A detailed methodology on the evaluation of SPI, SPEI and SNEPI can be obtained from McKee et al (1993), Vicente-Serrano et al (2010) andSingh et al (2021), respectively.
In the present study, SPI, SPEI and SNEPI are estimated at the 3 month timescale as this timescale acts as a proxy to represent meteorological droughts. Further, a 3 month timescale facilitates seasonal drought analysis as the average seasonal length Figure 1. The spatial extent of the climate zones over India (tropical monsoon (Am); humid subtropical (CWa); hot desert (BWh); hot semi-arid (BSh); tropical wet and dry (Aw); subtropical highland (CWb); cold desert (BWk)). The grey region, representing the cold desert climate is not considered due to unreliable data. observed over India amounts to 3 months. This does not reduce the relevance of other timescales, that may be required for different applications.

Divergence between drought indices 2.3.1. Divergence description
The inherent divergence present among various drought indices SPI 3 , SPEI 3 and SNEPI 3 is evaluated as, (1) The terms D 1 , and D 2 , are consistently used throughout the study to evaluate the index divergences. Since the possible impacts of evaporative losses and precipitation characteristics on drought behavior would be explored through D 1 and D 2 , respectively, a comparison between SPI 3 and SNEPI 3 , is not considered in this study to avoid redundancy. A negative divergence (−D 1 or −D 2 ) would imply more severe drought conditions represented by SPEI 3 , for both D 1 and D 2 . A positive divergence (+D 1 or +D 2 ) would imply dry conditions by SPI 3 and SNEPI 3 for D 1 and D 2 scenarios, respectively.

Differential divergence index (DDI)
Different divergence characteristics like magnitude, duration, and frequency are amalgamated for the first time into a single novel index called Divergence Index (DI). The annual maximum (minimum) series of the magnitude (M i ) of positive (negative) divergence is evaluated along with its respective duration (D i ). A decade-wise moving window, with a one-year lag is considered to evaluate the frequency (F i ) of ±D 1 , ±D 2 . The entire range of divergence characteristics were split into four quartiles and each quartile was given a score (s i ) in ascending order based on the magnitude range. For instance, a grid observing a divergence magnitude, duration, or frequency from the first, second, third and fourth quartile, would be given a score of 1, 2, 3 and 4, respectively. Weights (w i ) are assigned to define the importance of each divergence characteristic in the formulation of DI. Since unequal distribution of weights may lead to a bias towards a particular divergence characteristic, w i = 1 is considered. Further a weighted average of these scores would give DI, as seen in equation (3): where, i = number of elements in a time series of interest. Further, DDI is obtained by subtracting the DIs for the negative case from the positive one, to determine the dominant side.

Extraction of spell characteristics
Spell characteristics are obtained by dividing the daily rainfall time-series of a given location into wet (consecutive wet days) and dry (consecutive dry days) periods.
where, p i = daily rainfall in a wet spell, WS dur = maximum length of wet spell. DS dur implies the maximum length of dry spell. WS mag considers the maximum depth of rainfall in a wet spell. WS int is the maximum depth of wet spell divided by its respective length. over a given spatial domain. It works by decomposing the time-series (X t ) over a grid into various oscillatory components known as intrinsic mode functions (IMFs), as described in equation (6). The subsequent elimination of IMFs (Ij(t)) results in the formation of the final residue (R(t)), which is either monotonic or contains one extremum, and is free from the presence of oscillations,

Non-linear trend analysis and test for significance
Trend at a given temporal realization can be obtained as the difference in the residual time series from the present time (R t ) to the reference time (R 0 ) (Ji et al 2014), given by equation (7), The significance of these MEEMD trends was evaluated using Monte Carlo Simulation (Ji et al 2014). Here we compute 10 000 samples of red noise time series having the same temporal length as the time series of interest. The slow varying/residual component of the red noise time series is evaluated along with its empirical PDF. The trend at any point is considered significant if it falls within the 10% significance level of the generated PDF. More details on the Monte Carlo Simulation can be obtained from Vinnarasi et al (2017).

Attribution analysis
PCA (Jollifie 2002) is a multivariate analysis method that produces uncorrelated principal components (PCs) as linear combinations of correlated variables. The candidate variables are standardized to reduce the sensitivity of PCA to larger magnitude ranges. Furthermore, the covariance matrix (C), eigen values (E i ) and the percentage variance explained by each PC (V PCi ) are evaluated. More details on the PCA methodology can be obtained from supplementary section S1. While the first PC explains the maximum variance, the last explains the least. Further, E i is utilized to obtain the respective eigen vectors, also called the loading matrix (L), representing coefficients of the linear equations, relating the PCs with the candidate variables V 1 , V 2 , …, V n , given as, While the load magnitudes represent the strength of the relationship of a variable with the corresponding PC, the sign represents the direction. It is to be noted that if the column-wise variation of loads for certain variables is similar, the variables tend to provide similar contributions to the total variance that the PCs collectively explain. Such variables tend to cluster together in the PC co-ordinate system (Kovač-Andrić et al 2009) (see figure S1, supplementary). Using this concept, the candidate variables in the study are divided into two variable groups, (1) D 1 , R LW , R SW , SH, W S , T air, and (2) D 2 , DS int , WS mag , WS dur , WS int . Each group undergoes PCA and variables behaving similarly with D 1 and D 2 in the group process represented by the PCs, are identified. With this exercise we would obtain variables behaving similarly to D 1 and the ones behaving similarly with D 2 . More details on the clustering mechanism can be obtained in section S2 of the supplementary.

Evolutionary characteristics of divergence
3.1.1. DDI Variation of spatially averaged annual maximum (SAAM) magnitude (blue line) and duration (red markers) for +D 1 and +D 2 is illustrated in figures 2(a) and (b), respectively. Similarly, the SAAM magnitude (green line) and duration (pink markers) for −D 1 and −D 2 are presented in figures 2(c) and (d), respectively. Divergence frequency (spatially averaged at decadal time window with one-year lag) for ±D 1 and ±D 2 is presented in figures 2(e) and (f), respectively. Cyan (red) bars represent variation of positive (negative) divergence frequency, for both D 1 and D 2 . While the divergence characteristics (magnitude, duration, and frequency) of the positive sides have significantly declined (as seen in figures 2(a), (b), (e) and (f)), vice versa is observed for the negative side (as seen in figures 2(c)-(f)). The significant Mann Kendall's trend (Mann 1945, Kendall 1975 and Sen's slope (Sen 1968) for different divergence characteristics is tabulated in supplementary table S1. DDI displays a declining (shades of blue) trend with spatial majority for both D 1 (figure 2(g)) and D 2 (figure 2(h)). This declining trend implies a possibility of stronger negative divergence for both D 1 and D 2 , in the near future. The presence of a negative divergence, by definition in section 2.3.1., implies the projection of more intense droughts by SPEI. Thus, SPEI in general negatively disagrees from reality, assuming that SNEPI, by methodological definition, is a more realistic drought index (Singh et al 2021).

Trend analysis and seasonal behavior of divergence
Non-linear MEEMD trend exhibited by D 1 and D 2 from 1951 to 2016 is presented in figures 3(a) and (b), respectively. Shades of yellow-orange (blue) highlight the positive (negative) trend. There is an evident reduction in the trends exhibited by D 1 and D 2 , implying more intense droughts being projected by SPEI, annually. However, seasonally the trend behavior in D 1 and D 2 is quite disparate. The spatial agreement of the annual non-linear trend for the window of 1951-2016, with different seasons, in terms of the percentage match area, is presented in figures 3(c) and (d), respectively. Percentage area match implies what percentage of grids showing negative trend annually are also showing negative trend seasonally, and vice versa. For D 1 , the overall trend behavior largely matches (up to 70%) with all the seasons, represented in green color. However, in case of D 2 , the SM season experiences a significant mismatch (72%), highlighted in red. The temporal evolution of the spatial match of the annual trend behavior of D 1 and D 2 with those of the respective seasons is illustrated in figures 3(e) and (f), respectively. Both figures highlight a disparate divergence behavior of the SM season (blue line). Since the match percentage during the SM season for D 2 is the lowest, it signifies that a major portion of India is experiencing a positive D 2 trend during the SM season. Such trends points towards the persistence of drought like conditions in the SM season, which is aptly captured by SNEPI.

Evolution of climate variables
The evolutionary characteristics of D 1 and D 2 share close functionality with various climatic variables like PET and precipitation standard deviation (P SD : represents precipitation variability. A high standard deviation implies more spaced-out values from the mean both on the higher and the lower end, thus covering in this case skewed precipitation events), respectively. The non-linear MEEMD trend of PET and P SD is examined from 1951 to 2016. While PET shows an increasing trend at all grid points ( figure 4(a)), there is an increase in P SD in over 65% of the area ( figure 4(b)). Furthermore, the seasonal dispersion of PET trend, illustrated in figure 4(c), reveals that all seasons (SM, WM, W and S) collectively experience an increasing trend. P SD (figure 4(d)), on the other hand possesses a strong enhanced trend in the SM season, when compared to other months, thereby explaining the low percentage of match of the SM season with other seasons and annual divergence trend as illustrated in figure 3. Similar patterns are noticed across India, wherein the enhanced trends of standardized PET and P SD during the SM season are highlighted ( figures S3 and S4, supplementary).
Further, the functional dependency of PET on variables like R LW , R SW , SH, W s and T air , is not unknown, thereby identifying them as the potential independent variables causing D 1 . D 2 , on the other hand is driven by the seasonal (3 month) maxima of DS dur , WS mag , WS dur , and WS int . The long-term MEEMD trend, for the time window of 1951-2016, for each of these formative variables is illustrated in figures 4(e)-(m). While there is a countrywide increase (shades of red) in the trend exhibited by R SW (figure 4(f)), SH (figure 4(g)), T air (figure 4(i)), WS int (figure 4(l)) and DS dur (figure 4(m)), a reduction (shades of blue) is visible for R LW (figure 4(e)), W s (figure 4(h)), WS mag (figure 4(j)), and WS dur (figure 4(k)). Increased trend in R SW , T air, WS int, and DS dur points towards the aggravated impacts that climate change has recently brought to various meteorological variables, leading a warmer future which is likely to experience erratic precipitation events characterized by high intensity, immediately superseded by long dry spans. However, which variable amongst the group actually explains the evolution of D 1 and D 2 , is what drives our interest.

Explaining the evolution of divergence
The first three PCs are selected, based on the description of a collective variance threshold of 90%. The load clusters on the PC co-ordinate system are observed using six representative grids, randomly sampled over India, each representing a distinct climate zone. The variables clustering around D 1 , and D 2 are presented in figures 5(a) and (b), respectively. It is observed that R SW or T air , cluster around D 1 , and variables WS int or WS dur cluster around D 2 . This means that the contributions of D 1 , R SW and T air , and D 2 , WS int and WS dur , to PC 1 , PC 2 and PC 3 is of a similar kind, implying similar behavioral evolution. Thus, among several contributing variables considered, D 1 is most likely being influenced by R SW and T air , and D 2 by WS int and WS dur .
Similar patterns were obtained when the clustering procedure was extrapolated across India. It is observed that R SW , clusters the closest around D 1 , unanimously across India ( figure 5(c)). In case of D 2 , while WS int clusters in over 65% of the region, patches of WS dur also dominate in the extremely wet climatic zones ( figure 5(d)). The seasonal dispersion of the divergence clusters across India are illustrated in figures S5 and S6 of the supplementary, with their percentage area-wise summary being described in figures 5(e) and (f), respectively. While D 1 unanimously clusters around R SW across all the seasons, three out of four seasons experience the formation of D 2 -WS int clusters. The SM season, however, experiences a D 2 -WS dur clustering, with spatial domination. Further, in the WM season, the D 2 -WS dur cluster is confined to the southeastern part of India which experiences monsoon during the winter

Discussion
The spatiotemporal variation of divergence characteristics and DDI revealed an increasing trend of the negative divergence, implying more intense droughts being projected by SPEI. Similar results were also obtained by Stagge et al (2017) over Europe. The seasonal decomposition of annual divergence trend, revealed a disparate trend in D 2 during the SM season. While most seasons reciprocated the country wide annual negative D 2 trend, the SM season witnessed a spatial majority of positive D 2 trend, implying more intense monsoon droughts captured by SNEPI. A possible reason for such a strong signal of monsoon droughts could be a significantly increasing trend in the precipitation standard deviation during the SM season. A higher precipitation standard deviation is indicative of a greater number of intense wet spells in a given time frame. With a higher concentration of annual precipitation in the SM season, and a distinct divergence behavior, accurate detection of droughts in this season has become paramount, especially for a country like India, where the agriculture is largely rain-fed.
The absence of a dense network of observed drought data over India limits validation of results with ground truth. Also, the drought evaluation mechanism over India shares close proximity with the SPI methodology (Shewale and Kumar 2005). Thus, disagreements concluded among indices, in this study, do not imply that one is better than the other. We work on a simple assumption that SNEPI is more realistic as it captures the impacts of changing precipitation characteristics and increasing global temperatures in drought monitoring. Thus, on observing a strong disagreement of SNEPI with SPEI in the SM season, we suggest an enhancement of the operating value of the index during heavy rainfall months. It is to be noted that though SNEPI is more robust in the changing climate scenario, it is not free from limitations. Since SNEPI evaluates drought through the daily climatic variables, its application worldwide is subjected to the availability of daily data. Further the possibility of exploring non-stationarity is left unexplored while computing SNEPI.
A PCA-based attribution analysis aided the identification of variables sharing behavioral likelihood with D 1 and D 2 , through formulation of clusters. It is revealed that R SW clusters around D 1 , unanimously across India. This is expected as R SW represents the quantum of solar energy received by the Earth's surface from the Sun, governing the exchange of PET from the land to the atmosphere. Further, D 2 displays clusters with WS int in general and with WS dur in the monsoon season/climate zones, thereby highlighting the role of wet spell characteristics in defining seasonal drought events.

Conclusions
The present study attempts to unravel the seasonal character of the inherent divergence persisting among different meteorological drought indices like SPI, SPEI and SNEPI, along with the relevant variables dictating their evolution. Staged over a region with strong seasonality and a tropical climatology, the results from the study can be extended to similar climatological patches across the globe. It is revealed that temperature led traditional drought indices have been persistently evaluating stronger dry extremes. Further, the trend in the variability of precipitation spiked during the summer monsoons, indicating a possible disparate behavior of the season. This peculiar behavior of the monsoon season was also highlighted by the disagreements reported between SNEPI and SPEI, the former being able to report excessive season specific dryness which the latter overlooked. Thus, SNEPI successfully unravels the dominant occurrence of severe monsoon droughts. Further, the attribution analysis revealed a spatiotemporal dominance of shortwave radiation in driving the disagreements between SPEI and SPI. However, the divergence between SNEPI and SPEI was primarily driven by the wet spell characteristics (duration and intensity), with the relationship strengthening in monsoon season and tropical climate zones. Conclusively, in the changing climate, the pertinence of drought indices should be critically questioned before deploying them for operational usage, especially for regions with a monsoon driven agrarian socio-economy.