Observed trends in timing and severity of streamflow droughts across global tropics

We considered the daily streamflow observations throughout the tropics (23.5° N–23.5° S), with a buffer area extending from 30° N to 30° S, to accommodate even those river basins for which the headwater lies outside the tropics. This allows the inclusion of over 800 catchments that are archived across the Global Runoff Data Centre (GRDC 2015) and the India-Water Resource Information System (WRIS; https://indiawris.gov.in/wris/#/) repositories. The gauging sites are selected based on the following criteria: (A) the catchment size of at least 1000 km² or more. (B) A minimum of 70% of daily discharge time series must be available. (C) The catchment outlet location should not fall within a 10 km radius of a 'Mega Dam' (i.e. Dam height >15 m, Best 2019) (see supplementary information (SI); SI 1.1 for details). Following earlier assessments (Konapala and Mishra 2020, Kam 2021), we have identified catchments with near-natural flow and the discharge time series, which is least impacted by human influence. To this end, we have eliminated stations within a 10 km radius of a mega-dam (Kam 2021). Based on these criteria, 472 stream gauges were selected initially. Further, considering quality of streamflow records and availability of streamflow drought events, we scrutinize 135 sites over 19 IPCC reference climate regions across the tropic (figure S1, tables S1 and S2). Consistent with the literature (Fritsch and Carlson 1980, Gnauck 2004), we estimated the missing records in the streamflow series by time series interpolation using a shape-preserving piecewise cubic polynomial function at a daily time step (see SI 1.2 for details). On comparing observed versus infilled records, we find a good agreement between the two time series (see table S3 and figure S2).

We compared streamflow drought responses over two non-overlapping climate normal periods (CNPs): the baseline climate period is referred to as CNP1 (1961–1990), whereas the recent climate period is denoted CNP2 (1991–2020). These CNPs represent a 30 year average of global climate as described by the World Meteorological Organization (WMO 2020). We use the daily variable threshold approach for the identification of streamflow droughts. Considering a relatively low threshold, for instance, a threshold value of 70% would result in several drought events, which may not qualify for extremes. In contrast, considering a relatively higher threshold, for instance, a 90% threshold value would yield a very few drought events. Therefore, similar to earlier assessments and being the most common choice (Heudorfer and Stahl 2017, Ahmadi and Moradkhani 2019, Ganguli et al 2022a, Brunner et al 2023), we have adopted a threshold of an 80% exceedance probability value for streamflow drought identification. Following previous studies (Brunner et al 2021, 2023) and being the most common choice, we implemented a centered moving average of 31 days (d) to smooth the variable threshold time series. Further, a smoothing window excludes several short-term deficits resulting from variable thresholds, which cannot be considered droughts (Ahmadi and Moradkhani 2019). We further investigate the sensitivity of streamflow drought characteristics at different smoothing windows, such as 15 (±7 d), 21 (±10 d), 31 (±15 d), and 51 (±25 d) across a few selected representative catchments (figure S3). We find that the onset day and the drought frequency (table S4) remained relatively constant across catchments irrespective of the window lengths. We identify streamflow drought, if the discharge value is below the daily variable threshold persistently for 30 d or more (Brunner et al 2023). Next, we identify primary drought attributes, such as time to onset, duration and deficit volume (figure 1). Since streamflow timing follows a cyclical continuum, i.e., the ordinal day of an event at the end of a year can be close to an event occurring at the start of the following year, circular statistics are applied to investigate streamflow drought timing. We performed at-site trend analysis for individual drought properties using Sen slope (see SI 1.3 for details), while trend significance is determined using standard bootstrap-based procedure (Yue and Pilon 2004).

Zoom In Zoom Out Reset image size

We further explored, the nonlinear changes in joint occurrences of deficit volume and onset pattern over the two CNPs considering nonparametric kernel density estimators (Van Loon et al 2014). For modeling drought deficit volume, we evaluate the applicability of several probability distribution families, including Gamma, GEV, log-normal and Log-logistic distributions. However, for modeling drought onset, we use the von Mises distribution (see SI 1.4). Following earlier assessments (Dhakal et al 2015, Rutkowska et al 2018, Dhakal and Palmer 2020), von Mises distribution is used for fitting drought onset due to its (i) analytical tractability over other complex circular distributions, and (ii) closer resemblance to a simple linear Gaussian distribution, but with the unique added value that its coordinates exist on a circular plane, which is the case for random variables describing the timing of an event. The circular distribution is focused on the perimeter of a circle with unit radius (Jammalamadaka and SenGupta 2001). To assess the uniformity of the drought onset data within the circular framework, we performed Rayleigh's test, which is a powerful statistical tool for measuring the degree of concentration or dispersion in circular data (Villarini 2016, Wasko et al 2020a). Table S5 shows for almost all gauges onset records qualify for circularity, which further confirms the applicability of von Mises distribution for modeling drought onset. For drought deficit volume, we select the best probability distribution across CNPs (table S6) following the minimum value of the Akaike Information Criteria with the small sample corrections (Chakrabarti and Ghosh 2011). Table S7 shows the fidelity of different distribution in fitting drought deficit volume for a few selected representative catchments. Finally, we evaluate the goodness-of-fit of selected distribution using (i) Kolmogorov–Smirnov test for drought deficit volume (table S8) (ii) Watson test for drought onset time (table S9). The goodness-of-fit was assessed at 10% significance level with p-value > 0.10.

Next, we compared the probabilistic shift of individual drought properties, i.e., onset time and deficit volume during recent versus the retrospective CNPs across each catchment using the entropy-based Kullback–Leibler (KL) divergence metric (Santhosh and Srinivas 2013, Ji et al 2022). Moreover, we assess the nonlinear changes in the dependence strengths using the non-parametric form of total linear-circular dependency (see SI 1.5; Mardia and Jupp 1999) and an extremal dependence metric, Capéraá–Fougères–Genest (CFG; Frahm et al 2005). The total dependency metric (TDM; linear-circular dependency) can better represent the pairwise clustered behavior of two variables towards the center of distribution. In contrast, the Upper tail dependency metric (UTDM) focus on the likelihood of two random variables simultaneously reaching extreme values (Hao and Singh 2016).

The comparative assessment of trend in spatial pattern of onset, deficit volume and durations for two non-overlapping time windows (1961–1990 versus 1991–2020) are presented in figure 2. In the recent time window, 43% (58 out of 135) of catchments show increasing trends in deficit volume, and 68% (92 out of 135) of sites show a delayed onset. Further, we found only 27% (37 out of 135) sites depict an increase in deficit volume accompanied by delayed onset. The mean onset time of 60% of sites are clustered between January and March, with the interquartile range (75th–25th percentile) of median deficit volume remains 11.5 mm in CNP1 (figures S4(a) and (c)). In contrast, in CNP2, the mean onset time for around 65% of sites spread between January and April, while the number of sites showing onset during April in CNP2 is three times that of CNP1 (figure S4(b)). Further, in CNP2, we find a decrease in the median deficit volume across most regions, with an interquartile range of 8.5 mm, except in the South American Monsoon zone, where there has been an increase in deficit volume in recent period (figure S4(d)). In South America region, >50% of catchments showed a delayed onset (figure 2(a)). Interestingly, overall >60% of catchments illustrate significant (p-value < 0.1) decreasing trends in deficit volume in CNP1 (figure 2(b)). On the other hand, the percentage share of sites showing increasing trends in deficit volume with decreasing trends in durations are nearly equal (figure 2(c)). Compared to the CNP1, we find in CNP2, more sites show a delayed arrival of streamflow droughts accompanied by an increase in deficit volume with a longer drought duration (figure 2; right panel). A persistently prolong drought with an increase in deficit volume is apparent in 14% (10 out of 73) of catchments in South America and in a few catchments in northern Australia, such as Gregory Downs. In contrast, droughts over Southern Africa tend to become lesser severe with decrease in duration (figures 2(b) and (c)). Over the eastern coast of South Asia, an earlier onset of streamflow droughts is generally accompanied by a shift towards longer drought duration during CNP2 (figures 2(a) and (b): right). The South American continent shows the maximum changes in drought attributes (figure S5). Furthermore, the increase in relative frequency (number of events/year) of droughts are more prominent across the Northeast and South American monsoon region and the extratropical region of eastern coast of Australia (figure 3). However, there is a decrease in relative frequency of droughts in Western Africa, eastern coast of South Asia and a few catchments in Northern Australia (figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 portrays joint occurrences of shifts in onset time versus deficit volume of streamflow droughts. While the marginal distribution (shown using density plots) of changes in deficit volume and onset in recent time window (CNP2) shows a decrease in variance, only the shift in variance of onset time is statistically significant (p-value < 0.05 with F-test). This is further accompanied by a shift towards delayed onset as reflected by a rightward shift in onset PDF. In addition, the deficit volume change shows a leftward elongated tail, accompanied by a shift towards decrease in deficit volume. The trends in onset time versus changes in deficit volume mostly follow bivariate probability density fields (figure 4). A few exceptions include those contributing to outliers, such as a few catchments across the South American monsoon region in CNP1 (marked using a filled square) and North-East South America (marked using a filled circle) in CNP2, which show an increasing trend (significant) in deficit volume. Further, a few gauges show a non-significant decrease in deficit volume that is accompanied by a shift towards earlier onset. The slope in bilinear regression, incorporating the joint anomalies of onset time versus deficit volume, shows an approximate 19% increase from CNP1 to CNP2 in the second quadrant, confirming a further shift towards a decrease in deficit volume with a delayed arrival of drought.

Zoom In Zoom Out Reset image size

The probabilistic shift in distribution of drought properties is assessed using the KL divergence metric. Figures 5(a) and (b) shows changes in distribution of onset time and deficit volume. The higher values of KL divergence indicate a substantial shift in the shape of the distribution. Nearly one-third of catchments show a large shift in the deficit volume, while only 8% show a substantial change in the onset timing. Interestingly, sites in the Northern Hemisphere show a smaller shift as depicted by the mean KL divergence values for both deficit volume and onset time to Southern Hemisphere. However, the fraction of sites showing a large shift in deficit volume is nearly the same across both hemispheres. The interquartile range of the KL divergence value for deficit volume in the Southern Hemisphere is approximately double that of the Northern Hemisphere, indicating an increase in variability in drought deficit volume in CNP2. Approximately, 13% of sites show minimal shift (KL divergence ∼ 0) for both onset and deficit volume. These sites are also equally spread across both hemispheres. However, ∼50% sites show minimal shift in deficit volume over two CNPs in the Northern Hemisphere. Furthermore, the sites across North-East South America and the South American Monsoon region show significant shifts in the distribution of onset timing and deficit volume (figures 4 and 5).

Zoom In Zoom Out Reset image size

Figure 6(a) compares the TDM using linear (deficit volume)-circular (onset timing) dependence metrics of CNP2 versus CNP1. The difference in the frequency of TDM (figure 6(a); right) shows a weak asymmetric pattern with a negative skewness of close to zero. This suggests a little shift in TDM metrics over the two non-overlapping CNPs. However, at at-site level, >50% of sites across Asia and Africa, and >45% of sites in South America show notable positive changes in TDM dependence strengths, indicating stronger correlations in the recent period. In contrast, around 66% of sites in North America and 57% of sites across Australia display negative changes, indicating a decrease in dependence strengths in recent periods. Overall, 23% more sites show the dependence strengths of >0.5 in CNP2 compared to CNP1. In particular, Southern Asia, North-East South America, and North Australia regions showed stronger TDM in the CNP2 compared to CNP1 (figure S6(a)).

Zoom In Zoom Out Reset image size

The spatial variation in the UTDM shows a bimodal nature of distribution for CNP2 (figure 6(b)). The leftward tail of CNP2 shows a slight more increase in frequency (i.e., counts) compared to CNP1, suggesting a leftward shift of overall distribution. The decrease in strength of dependence in CNP2 shows a decrease in bivariate hazard over South-East Asia, North-West South America and South East South America regions. Taken together, considering the UTDM >50% sites show a decrease in dependence strengths in CNP2. Although, this overall shift in CNP2 relative to CNP1 is not statistically significant. Nevertheless, a few sites show an increase in UTDM, explaining an increase in bivariate hazard. Such sites are clustered around Eastern-South Africa, Western-South Africa, East-North America and North-Central America regions (figures 7 and S6(b)). More sites show higher uncertainty bounds in UTDM than that of TDM, as depicted by the interquartile range. This could be due to the limited sample size at the upper tail quadrants of bivariate pair onset time-deficit volume. In contrast, in TDM all samples are considered for determining dependence strength across the central part of the bivariate distribution. Considering UTDM, notable uncertainty is observed in Western Africa, followed by the South-West South America regions (figure 7). On the other hand, considering TDM, North-Central America shows the highest uncertainty, followed by West-Southern Africa.

Zoom In Zoom Out Reset image size