Introduction

Droughts are highly complex and less understood phenomena as compared to other weather and climate extremes, that propagate slowly but with long-lasting and devastating effects1. It has now become a recognized fact that rising temperature and climate change cause the frequency and severity of drought events to increase in many regions around the globe2,3. While a simple definition of the onset of a drought is reduction in the atmospheric/climatic water balance (also referred as a meteorological drought)4 i.e. unusually long time period without precipitation, it can manifest into a more complex natural hazard. For instance, a prolonged deficit in the climatic water balance can lead to reduced streamflows and induce a hydrological drought; and/or it may lead to significant reduction in soil moisture and thus cause an agricultural drought5. While meteorological, hydrological and agricultural droughts are direct manifestations of reduction of water in the hydrological cycle, they in-turn lead to numerous adverse socio-economic, socio-political and environmental impacts, including but not limited to water insecurity, food insecurity, economic loss and water quality deterioration etc6., which will hit more severly in a future warmer world. The extent of these impacts depends on the onset, duration and severity of different types of droughts (e.g., meteorological hydrological etc.)7. Hence, it was never as important to efficiently measure, monitor, understand and predict these drought characteristics as it is now in these climate emergency times.

The use of standard indices is a highly popular and effective mechanism in monitoring and predicting the onset, duration and severity of droughts1,8,9,10. The Standard Precipitation Index (SPI)11,12, which is based on the deviation of observed precipitation with respect to the climatological average over a given time period, is the most prominent index in this regard, that is recommended by the World Meteorological Organization (WMO) as an effective indicator for inclusion in national and regional Drought Early Warning Systems (DEWS)4. SPI is a part of numerous national and regional DEWS, for instance, the North America Drought Monitor (https://www.ncdc.noaa.gov/temp-and-precip/drought/nadm/), United States Drought Monitor (https://droughtmonitor.unl.edu), European Drought Observatory and Pakistan’s National Drought Monitoring Centre (NDMC) (http://ndmc.pmd.gov.pk/ndmc.php). Drought monitoring agencies and DEWS across the globe typically use multiple such indicators to monitor different types of droughts. NDMC Pakistan, for instance, uses SPI for meteorological drought monitoring and Normalized Difference Vegetation Index (NDVI)13 for agricultural drought monitoring6. provide a comprehensive review of different drought indicators used by national and regional state agencies in operational DEWS. A key outcome of this review is that development of effective DEWS and the selection of appropriate drought indicators within, requires a deep understanding of the inter-linkages between water deficits in the hydrological cycle, human activities and associated impacts on society and the environment. This notion is echoed by WMO4,14 who argue that there is no universal definition of droughts, and hence, development of DEWS and effective use of drought indicators should incorporate an understanding of the regional and local contexts that are impacted by water deficits, e.g., water use, water management practices and water politics etc.

In Pakistan drought events are strongly linked to water security, food security and agro-economics of a majority of the country’s population. More than 60% of Pakistan’s population is rural15 and mostly resides in the semi-arid to hyper-arid plains of the lower Indus basin16. Agricultural activities of these rural dwellers are heavily dependent on streamflows originating from the Upper Indus catchments, that are diverted towards vast crop fields via the massive Indus Basin Irrigation System17. These streamflows generally exhibit high seasonal variability and prone to extreme events and climate change18,19. Historical droughts in the Upper Indus catchments have had extreme adverse effects on Pakistan’s economy and especially the rural sector in the past. For instance, the 1998–2002 droughts severely affected the rural population and resulted in negative growth in the agricultural sector, which had seen significant growth in the prior decade20. Remnants of this worst drought of the country’s history17 can still be seen in the parts of the Baluchistan Province in the form of several abandoned fruit orchards (Eye-witnessed). However, data, research and analyses on retrospective drought assessment, and development of drought warning systems and management plans for Pakistan are very limited15,21. The importance of further research in this regard is also observed by the recently coined and Pakistan’s first National Water Policy.

The use of multiple standard indices for development of DEWS for Pakistan is also relatively unexplored. We found only a few relevant studies in recent literature, that used SPI/SPEI to investigate characteristics of historical droughts in Pakistan22,23,24,25,2622. used a new agricultural Standardized Precipitation Index (aSPI), along with SPI and SPEI to compare them for meteorological drought assessment23. studied the propagation of meteorological drought to hydrological drought across various sub-basins of the Upper Indus Basin (UIB) using the combined application of principal component analysis (PCA) and wavelet analysis24. performed a historical spatial analysis of droughts in the Upper Indus Basin (UIB) for the period of 1991–201725. computed SPI time-series at different accumulation periods using observed precipitation data from 61 meteorological stations across Pakistan, to analyze duration, severity and frequency of historical droughts.

As reported by22, the SPI seems more popular than other indices e.g. SPEI or PDSI as it is more appropriate to use in all climate regimes and situations where there are challenges related to quality and availability of long-term data. Also, in Pakistan most of the droughts are precipitation driven, hence SPI, which is simpler than other indices as it only requires precipitation data, seems more robust choice26. computed SPI using the gridded APHRODITE precipitation dataset27 to discover spatio-temporal trends and cyclical behaviours of drought recurrence. There is no prior literature on the investigation of multiple indices for assistance with drought monitoring and early prediction using a distributed approach in terms of multiple important sub-basins of the UIB that are mainstay for downstream water supplies. A key research gap here is in the investigation of the coherence between meteorological and hydrological droughts, especially in the upper catchments of the Indus basin. As mentioned earlier, a key identifier of drought onset in Pakistan is the streamflow originating from Upper catchments of the Indus basin. Hence, identification of relevant indicators that can assist in analyzing hydrological drought characteristics at key streamflow locations in the Upper catchments of the Indus basin, can be extremely beneficial towards improving the effectiveness of existing DEWS in the country. Moreover, if coherences between hydrological and meteorological droughts can be unearthed through these indicators, they could also be integrated to provide early drought warning mechanisms.

The primary purpose of this study is to investigate the coherence between meteorological and hydrological droughts in four key upper catchments of the most important geographical segment of Pakistan i.e. the famous UIB. In order to investigate this coherence, we will utilize the Standardized Precipitation and Evaporation Index (SPEI)28, (computed using gridded precipitation and evaporation datasets) as an indicator of meteorological droughts, and the Standardized Streamflow Index (SSI)29 as a measure of hydrological droughts (computed on observed streamflows at selected catchments outlets). We assess coherence between meteorological droughts and hydrological droughts by analyzing cross-correlations and lagged cross-correlations between SPEI and SSI for different durations. Hydrological flows in the Upper catchments of the Indus basin are highly seasonal. Hence, monthly cross-correlations are investigated to incorporate seasonality. The study supplements the existing knowledge in our assessment of lagged cross-correlations particularly in context of snow dominated regions. In addition to this our study also captured the impact of seasonality. We believe that the findings of this study may be very useful in operational drought monitoring and early warning within Pakistan.

Study area and data

Study area selection

Pakistan’s freshwater resources primarily come from the upper catchments of the Indus basin16, a trans-boundary basin that spans four countries, Pakistan, India, Afghanistan and China (see Fig. 1). In terms of climatology and land-use the Indus basin can be divided into two key zones, i.e., the Upper Indus Basin and the Lower Indus Basin.

Fig. 1
figure 1

Overview of spatial variation in climatology across the Indus Basin (Upper and Lower zones): (a) Average annual gridded precipitation for the period 1901 to 2018, and (b) Average annual gridded Potential Evapotranspiration (PET) for the period 1901 to 2018. Gridded precipitation and PET are derived from the CRU TS4.03 data set39,40.

The Upper Indus Basin is primarily a high elevation zone that includes the Himalaya-Karakoram-Hindukush mountain ranges and is the dominant source of renewable freshwater within the Indus Basin. The Lower Indus Basin includes the highly irrigated central Indus Plain where irrigation is made possible through diversion of surface flows via the Indus Basin Irrigation System (IBIS), and also includes deserts of the Indus plain and the Indus delta15.

Climatology of the Indus basin (Fig. 1) clearly illustrates the higher concentration of precipitation in the Upper Indus. This coupled with the low PET in the mountain ranges of the Upper Indus, makes the entire basin heavily reliant on the high-altitude upper catchments for renewable water resources, especially in Pakistan and the Punjab province in India16.

According to numerous water availability assessments15,16,30,31 more than 70% of renewable freshwater in Pakistan comes from the Upper Indus Basin catchments, and predominantly from upper catchments of the ’Western’ rivers of the Indus, namely, Chenab, Jhelum and Indus (see Fig. 2), and the Kabul river basin. ‘Western rivers’ is a term that originates from the Indus Waters Treaty (IWT), a transboundary treaty signed between India and Pakistan in year 1960 to resolve post-independence water conflicts between the two countries by amicably sharing the surface water resources of the Indus Basin. The IWT gave consumptive use rights of surface flows of the Ravi, the Sutlej and the Beas rivers of the Indus Basin (also called Eastern Rivers) to India, and consumptive rights of the Chenab, the Jhelum and the Indus rivers (also called Western Rivers) to Pakistan. The Kabul is a transboundary right bank sub-basin of the Indus River that spans Afghanistan and Pakistan (see Fig. 2), and is thus not a part of the IWT. Kabul river, however, is a critical source of freshwater in the lower Indus Basin of Pakistan15,32 and joins the Indus River in Pakistan downstream of Nowshera (see Fig. 2).

Fig. 2
figure 2

Overview of study area for meteorological and hydrological analysis, that includes the four western river catchments of the Upper Indus Basin, i.e., Upper Indus, Kabul Upper Chenab and Upper Jhelum, and the studied streamflow stations at respective catchment outlets, i.e., Nowshera, Tarbela, Mangla and Marala.

Figure 2 delineates the four Indus basin catchments selected for Pakistan’s drought analysis in this study, i.e., the Upper Chenab, the Upper Jhelum, the Upper Indus and Kabul. As mentioned earlier, these catchments account for more than 70% of Pakistan’s renewable freshwater availability15,31, and via the outlets shown in Fig. 2, i.e. Chenab at Marala Barrage, Jhelum at Mangla Reservoir, Indus at the Tarbela reservoir and the Kabul at Nowshera. The streamflow gauging stations at these outlets (excluding Nowshera) are also called rim stations33, since they essentially signal the beginning of the Indus Basin Irrigation System (IBIS) (also called the Indus River System) of Pakistan. IBIS is the world’s most dense irrigation system, and is essential to the water and food security of Pakistan, and accounts for more than 90% of Pakistan’s food production34.

Since water supply for IBIS is predominantly derived from the four catchments delineated in Fig. 2, drought conditions in these catchments have significant impacts on water security, food security and economic well-being of Pakistan25,31,35. Hence, these catchments are the focus of the meteorological and hydrological drought analysis conducted for Pakistan in this study.

Data

Observed climate data for Upper Indus catchments is limited, both spatially and temporally36. In such data scarce scenario, use of a gridded climate dataset is a feasible alternative. Prior studies on data scarce Asia Pacific catchments observe that observation based gridded data sets are better than reanalysis-based data sets, especially for hydrological modelling and historical drought analysis37,38. Hence, we use the observation based gridded CRU TS4.0339 climate dataset for our drought analysis. CRU is a monthly gridded dataset that shows good correlations with station observations in the South Asian region, both for precipitation and Potential Evapotranspiration40. Despite the multiple advantages it offers in a situation where observed station data is not sufficient, it has some limitations due to its 0.5° spatial resolution and dependence on sparse station networks. Therefore, it cannot fully replace a dense network of the long-term, good quality observed data.

We also use streamflow data of outlets of the four catchments discussed in the previous section and delineated in Fig. 2 for the hydrological drought analysis of this study. Monthly streamflow records from 1961 to 2018 for the relevant flow stations, i.e., Chenab at Marala Barrage, Jhelum at Mangla Reservoir (upstream), Indus at Tarbela Reservoir (upstream) and Kabul at Nowshera gage (see Fig. 2) are used in our hydrological drought analysis. Figure 3 shows monthly distribution of flows at these stations and depicts the high seasonality and variation in outflows from the Upper catchments of the Indus basin.

Fig. 3
figure 3

Monthly flow distribution (using streamflow data from 1961–2018) at the outlets of Upper Indus catchments of Pakistan: (a) Chenab at Marala, (b) Jhelum at Mangla, (c) Indus at Tarbela and (d) Kabul at Nowshera.

Methodology

Drought indicators

The drought analysis of this study makes use of standard indicators for monitoring, assessment and early prediction of droughts in Pakistan. While many indicators are available for drought analysis and early warning11,28,29,41, the choice of an appropriate indicator may depend on the type of drought (e.g., meteorological, hydrological agricultural drought etc.) under investigation and/or the purpose of drought analysis4,6,7. Since, the focus of this study is on analyzing meteorological and hydrological droughts within the context of Pakistan, we will use indicators that are designed for quantifying hydrological and meteorological droughts.

Standard precipitation and evapotranspiration index (SPEI)

The Standard Precipitation and Evapotranspiration Index (SPEI)28 is the indicator that is used for quantification and monitoring of meteorological droughts in the Upper Indus catchments analyzed in this study. SPEI is theoretically similar to the widely used Standard Precipitation Index (SPI), which is computed by fitting a parametric distribution on precipitation accumulation over a user-defined period, and subsequent transformation of resultant probabilities to standard normal variates11,12. SPI is strongly recommended by the World Meteorological Organization (WMO)4, and is the standard index used for meteorological drought analysis by many drought monitoring authorities across the world5 including Pakistan Meteorological Department’s (PMD) National Drought Monitoring Centre (NDMC).

However, a distinct advantage of using SPEI over SPI for meteorological drought monitoring is that it uses both precipitation and temperature data for calculation and thus, may better represent meteorological droughts and how they may manifest into hydrological and agricultural droughts425. argue that temperature abnormalities play a critical role in drought onset and propagation, especially in cold-climate catchments (such as the Upper Indus catchments of this study). For instance, higher than normal temperatures in winters can result in less snow accumulation and low snowmelt runoff in spring.

SPEI is calculated in the same way as SPI, however the underlying climate variable is the difference between precipitation and evapotranspiration (also called the climatic water balance)28.

Another advantage of using SPEI instead of SPI is that the underlying climate variable will avoid having too many zero values. Precipitation can be zero for many time periods, especially in semi-arid and arid climate zones, and this can cause difficulties in adequate fitting of distributions during SPI computation43. In a comprehensive distribution fitting assessment of both SPI and SPEI for Europe8, observed that fit of distributions was indeed better for SPEI, and recommended its use as an adequate drought index. Hence, we have chosen SPEI as the index for analyzing meteorological droughts for our selected arid zone catchments.

We calculate SPEI for the four catchments analyzed in this study (see Fig. 2) at both grid-level and catchment level using gridded climate data from the CRU data set39 (as discussed earlier). For catchment level SPEI calculations, gridded precipitation and Potential Evapotranspiration values are averaged over the catchment at each time interval before use. SPEI is typically calculated over a user-defined accumulation period between 1 and 24 months. Hence, the monthly resolution of the CRU data set is sufficient for our analysis. Moreover, we fit the log-logistic distribution for computing SPEI at different accumulation periods (as recommended by42, and use the R package ‘SPEI’ for our computations44 due to a number of advantages it offers e.g. time scales can be adjusted, efficiently and seamlessly processes gridded datasets, tailored for climate studies, validated in global drought studies (e.g45. and cited in IPCC reports.

Since SPEI results from different accumulation periods are discussed in this study, we use the term ’SPEI-x’ in subsequent discussions, where x denotes the accumulation period.

Standard streamflow index (SSI)

The Standard Streamflow Index (SSI)29,46 (also called the Streamflow Drought Index and the Streamflow Runoff Index) is used to quantify and analyze hydrological droughts at catchment outlets in this study. The SSI is a widely used hydrological drought index, particularly effective for monitoring drought conditions in regions with conditions similar to Pakistan. Other similar indices include the Streamflow Drought Index (SDI), Surface Water Supply Index (SWSI), and Standardized Runoff Index (SRI). While these indices are valuable, SSI offers distinct advantages over these, e.g. it is derived directly from streamflow data, providing a standardized measure that facilitates comparison across different temporal and spatial scales47,48,49. This standardization allows for consistent monitoring and assessment of hydrological droughts. Given these strengths, SSI is often preferred for hydrological drought monitoring, offering a reliable and efficient means to assess and compare drought conditions across various regions.

The SSI computation mechanism is similar to SPI and SPEI. However, streamflow is the underlying time-series data used and the recommended interval may vary between 1 and 12 months.

Within the context of Pakistan and especially the catchments being assessed in this study, analysis of hydrological droughts over accumulation periods of 1–6 months might be the most suitable50,51. This is because hydrological outflows from the four catchments of our study are central to seasonal planning and operations (conducted in cycles of six months) of the IBIS in lower Indus50,51. A major proportion of the outflows of Upper Chenab, Upper Jhelum, Upper Indus and Kabul catchments are diverted for agricultural consumption in the lower Indus via IBIS15. Since, Pakistan’s agricultural system operates on rotation of crops in two six-month seasons, called ‘Rabi’ (October-March) and ‘Kharif’ (April-September), seasonal water supply planning and operation of IBIS is conducted in cycles of six months33.

In hydrological drought analysis, selecting an appropriate probability distribution is crucial for accurately modeling streamflow data and computing indices such as the Standardized Streamflow Index (SSI)52. Unlike SPI and SPEI, there is no probability distribution that is universally recommended for computation of SSI53. However54, observe that distributions which are commonly used for hydrological frequency analysis may be suitable for computing SSI, especially the more flexible three parameter distributions. Moreover54, recommend that different distributions should be tested for separate months/accumulation periods before selection of an appropriate distribution for SSI computation. We compare five distributions, namely, log-logistic, Gamma, Pearson Type III, Lognormal and GEV, at 1-Month accumulation periods (i.e., for each month) for all four streamflow locations, to select an appropriate distribution for the SSI computations of this study.

Log-logistic distribution

This distribution is particularly effective in modeling data with heavy tails, making it suitable for representing extreme drought events.

Gamma distribution

Widely used due to its flexibility in modeling skewed data, the Gamma distribution is adept at handling non-negative datasets like precipitation and streamflow.

Pearson type III distribution

Known for its ability to model skewed hydrological data, this distribution is often employed in flood frequency analysis and has been recommended in various hydrological studies.

Lognormal distribution

This distribution is suitable for modeling positively skewed data and is commonly applied in hydrology for variables like rainfall and streamflow.

Generalized extreme value (GEV) distribution

GEV is particularly useful for modeling the extremes of hydrological variables, capturing the behavior of maximum or minimum values over specified periods.

Moreover, we use historical streamflow data from January 1961 to December 2018 for SSI time-series computations (at 1- to 12-month accumulation periods), and employ the R package ’SPEI’ for the computations44. In subsequent discussions, the term ’SSI-x’ is used to denote SSI with ‘x’ accumulation period. Therefore, SSI-3 denotes SSI values for a 3-month streamflow accumulation period. We have used the lag of 2 months for the SSI-3 in order to avoid the double accumulation. Once SSI is computed, cross-correlations between SSI and SPEI (at different accumulation periods) will be calculated to analyze how meteorological droughts may manifest into hydrological droughts.

Drought characteristics

This study also includes an indicator-based analysis of the characteristics of historical drought events (hydrological and meteorological) observed in the Upper Indus catchments under consideration (see Fig. 2). The definition proposed by12 is used to identify drought events, i.e., any period where an indicator’s (SSI or SPEI) value is consistently below zero and also goes below a pre-defined threshold is identified as a drought event. We use a drought threshold of − 1.0 in our analysis, and record the number of drought events for each catchment (see Table 1 and Sect. 4.2 for discussion). We also analyze statistics of two key drought characteristic variables, i.e., drought duration and drought severity for the identified drought events. Both of these characteristics are widely used by drought monitoring agencies and water resources managers12,55,56 to analyze historical drought events. Drought severity has been defined in different ways in prior literature12,55,57. We use the definition proposed by McKee et al. (1993), which is as follows:

Table 1 Summary characteristics (number of events and average drought duration) of historical meteorological and hydrological drought events observed in the study catchments. Characteristics are calculated using − 1 as the threshold for defining a drought event.
$$\:{S}_{e}=\sum\:_{i={j}_{e}}^{{x}_{e}}\left|{I}_{i}\right|$$

In the above equation, \(\:{S}_{e}\) is the severity of drought event \(\:e\), \(\:{I}_{i}\) is the drought index value for month \(\:i\), \(\:{j}_{e}\) is the starting month of drought event \(\:e\) (i.e., starting month \(\:i\) when \(\:{I}_{i}<-1\) and \(\:{I}_{i-1}>=-1\) (if exists)) and \(\:{x}_{e}\) is the ending month of drought event \(\:e\). Please note that the predefined threshold used for defining the onset of drought is set to −1 for our analysis (also used by Halwatura et al., (2015)). The purpose of the drought characteristic analysis is to assess the duration and magnitude (and the correlation between duration and magnitude) of historical drought events observed in the Upper Indus catchments, in order to consequently understand vulnerability of these catchments to drought events. In order to achieve this purpose, we analyze the correlations between drought duration and severity for all catchments, and also the empirical distribution of drought duration. Section 4.2 discusses the characteristics of drought events observed in the Upper Indus catchments, by also using the variables defined in this section.

Integrating meteorologic and hydrologic drought indices

In order to develop effective early hydrological drought forecasting and monitoring systems, it is important to understand if climatic indicators can adequately predict the onset of hydrological droughts5. Hence, a key question that this study aims to address is that can the onset of hydrological droughts in the Upper Indus catchments be predicted or characterized through meteorological drought indices? Numerous methods have been used in the past to investigate the coherence between hydrological and meteorological droughts58,59,60,61. A popular method in this regard is the use of cross-correlations and lagged cross-correlations between meteorological and hydrological drought indices. For instance58, observed strong cross-correlations and lagged cross-correlations between SPI and SSI for catchments in the UK. However, they did not consider the impact of seasonality in their analysis.

Climate variables and streamflows in the Upper Indus catchments of this study are strongly seasonal. For instance, average catchment outflows in the wet season of April-September (also called Kharif season) are 4–5 times higher than dry season flows (October-March; also called Rabi season)15. Even within the wet season, a significantly high proportion of flows occur in a 3-month window. This trend is illustrated in Fig. 3 via monthly flow distribution plots at the outlets of the catchments. It is important to incorporate the effect of seasonality in the analysis of coherence between meteorological and hydrological droughts15,55. Hence, we develop monthly cross-correlations and lagged cross-correlations between SSI and SPEI (at different accumulation periods) to (i) analyze how meteorological droughts may manifest into hydrological droughts and (ii) to identify appropriate lead times (if any; via lagged correlations), for early warning/prediction of hydrological droughts using SPEI.

Results and discussion

SSI distributions

The calculation of SSI for the streamflow sites of this study was preceded by selection of an appropriate distribution for SSI computations (for each streamflow site). As discussed in Sect. 3.1.2 and observed by53, there is no single distribution that is universally suitable for computing SSI for any streamflow site. Hence, we compared five different distributions, i.e., log-logistic, Gamma, GEV, Pearson-III and lognormal, separately, for the 4 streamflow sites of this study. Moreover, quality of distribution fit was analyzed for each month per site. Figure 4 provides a visual comparison of the fit of each distribution for the four sites (for selected months), by plotting density functions of fitted distributions against histograms for respective months and catchment outlets. While the difference between distributions is hard to ascertain from Fig. 4, overall performance of log-logistic is marginally better. In a more comprehensive assessment of different distributions for SSI computations54, also recommended log-logistic distribution as an appropriate option, if a single distribution is being used. Hence, we use the log-logistic distribution for computing SSI for all four streamflow gauges of this study.

Fig. 4
figure 4

Comparison of five theoretical distributions, i.e., log-logistic, Gamma, Pearson Type III, lognormal and GEV, that were tested for computing Standardized Streamflow Index (SSI), for streamflows at the outlets of Upper Chenab, Upper Jhelum, Upper Indus and Kabul catchments (see Fig. 2).

Drought characteristics

We analyze drought characteristics of the selected Upper catchments of the Indus basin by first analyzing the time-series plots of SSI and SPEI. Figure 5 shows the time-series plot of SSI-3 for streamflow gauges at the outlets of the four catchments. SSI-3 provides a good representation of historical droughts observed in the Upper Indus catchments. According to prior literature prominent drought events were observed in Pakistan in early 1970 s, mid 1980 s and in the years 1998–200215,25,26. These drought events are prominent in Fig. 5 as well. Another key observation from Fig. 5 is that many hydrological drought events are cross-correlated across the four basins. This observation implies a higher impact on downstream water users due to such cross-correlated events. Moreover, in such drought periods, water demands of the lower Indus basin may be highly dependent on surface reservoirs and may put further pressure on the already stressed groundwater resources62.

Fig. 5
figure 5

Time-series plots of 3-Month SSI (also called SSI-3) for streamflows at outlet of Chenab, Jhelum, Indus and Kabul rivers (see Fig. 2 for gauge locations). The droughts of early 1998–2002 are well represented by SSI-3.

The SPEI-3 time-series plots of Fig. 6 also prominently represent the historical droughts reported in history, including the 1998–2002 droughts, that are considered the worst in Pakistan’s history20. Moreover, Fig. 6 is also indicative of the coherence between SPEI and SSI for respective catchments (discussed in further detail in Sect. 4.3). The temporal oscillations (i.e., oscillations between negative and positive SPEI-3 values) in the basin-wide SPEI-3 values, however, are higher than in corresponding SSI-3 values. A plausible reason here could be that snowmelt and glacial melt dominate the hydrology of these Upper Indus catchments18,50. Hence, SPEI-x with a larger accumulation period could correspond better with SSI-3, in terms of oscillatory patterns and drought duration.

Fig. 6
figure 6

Time-series plots of 3-Month SPEI (also called SPEI-3) for the Upper catchments of Chenab, Jhelum and Indus rivers and Kabul river basin.

Table 1 provides a summary of characteristics of meteorological and hydrological droughts observed in the Upper Indus catchments, by reporting the number and average duration of drought events (using − 1.0 as threshold for defining a drought event; also see Sect. 3.2). Unsurprisingly, the number of drought events reduce with an increase in accumulation period, for both SPEI and SSI, whereas average duration of droughts increases with increase in accumulation period12,58. Table 1 also indicates that frequency of meteorological droughts (represented by SPEI) is generally higher than hydrological droughts (represented by SSI), whereas duration of hydrological droughts is relatively longer. The observation regarding drought durations is also evident in Fig. 7 where distributions of drought durations are visualized (via violin plots, that combine kernel density plots and boxplots), for both SSI and SPEI with different accumulation periods. Figure 7 also depicts that SPEI based drought durations are shorter than SSI based drought durations, for the same accumulation periods. Table 1; Fig. 7 imply that multiple meteorological droughts events can collectively result in a hydrological drought event. Hence, SPEI values with relatively higher accumulation periods may correlate better with corresponding SSI values for streamflow-based drought characteristic analysis. We explore these correlations in more detail in Sect. 4.3.

Fig. 7
figure 7

Empirical distributions of drought durations computed using (a) SSI with 1-Month, 3-Month and 6-Month accumulation periods (i.e., SSI-1, SSI-3 and SSI-6), (b) SPEI with 1-Month, 3-Month and 6-Month accumulation periods (i.e., SPEI-1, SPEI-3 and SPEI-6).

We also analyze the relationship between duration and severity (as defined in Sect. 4.2) of historical droughts observed in the Upper Indus catchments. To this effect, Figures S1 and S2 (in supplementary document) provide a visual overview of the strong mutual correlation between hydrological drought duration and severity. This relationship is logical, since, longer droughts tend to be more severe, and is also corroborated by prior probabilistic duration-severity analyses conducted using copulas57,63 that report high mutual correlation between drought duration and severity for other basins.

SPEI and SSI coherence

The SPEI and SSI-based drought analysis performed so far shows some coherence between (i) SPEI and SSI drought characteristics deduced from both of these metrics (see Figs. 5 and 6), and (ii) distributions of SPEI and SSI based drought durations. In this section we further explore the coherence between these indices through analysis of cross-correlations between them.

Figure 8 represents the overall cross-correlation between SPEI and SSI for all four catchments, and with different respective accumulation periods (1–12 months for SPEI and SSI-1 and SSI-3). The cross-correlations between SPEI and SSI are better for higher SSI accumulation periods (e.g., SSI-3). Moreover, cross-correlations are better when SPEI accumulation periods are greater than respective SSI accumulation periods. This is expected for the large catchments of this study (see Fig. 2) because these catchments have seasonal climates with sub-zero temperatures, highly variant topography, and snow accumulation in winter. Hence, climatic water deficits tend to accumulate and propagate slowly into water deficits in surface flow.

Fig. 8
figure 8

Heatmap of cross-correlations between SSI-1 and SPEI (for different accumulation periods) for the four Upper Indus catchments, i.e., Chenab, Jhelum, Indus and Kabul.

There are no SPEI and SSI accumulation periods where cross-correlation values in Fig. 8 are greater than 0.7, for any of the four catchments. Moreover, the cross-correlations are particularly poor for Upper Indus. A key reason for these overall low cross-correlation values is the lack of consideration of seasonality in their computations5. explains that in cold-region catchments with highly seasonal climates, the hydrological processes behind drought propagation are highly time-dependent. For instance, winter season droughts in cold regions may be triggered by abnormal temperature drops leading to early snow accumulation and low winter discharge. Or precipitation deficit in winter could lead to reduced snowmelt dominated discharge in spring5. Hence, it is important to analyze cross-correlations between SPEI and SSI for the cold-region catchments of this study, with inclusion of time-dependence and seasonality.

Seasonal trends

Figure 9 shows heatmaps of monthly correlation between SPEI (for accumulation periods between 1 and 12 months) and SSI (for accumulation periods 1 and 3) for the four Upper Indus catchments of this study. It is immediately apparent that monthly cross-correlations values in Fig. 9 are better than the overall cross-correlation values depicted in Fig. 8. For example, if we compare the correlation for Jhelum between SSI and SPEI, overall correlation for SPEI-6 is around 0.5, whereas monthly correlation for the months of February, March and April are around 0.7. There are multiple SPEI accumulation periods where SPEI cross-correlates strongly (i.e., where cross-correlation \(\:r>0.7\)) with SSI (at different accumulation periods), specifically in the late winter and spring months (Jan-May), for Chenab and Jhelum, and spring to summer months (Mar-Aug) for Kabul. For instance, in May, cross-correlations between SSI-6 and SPEI-8 are greater than 0.8 for Chenab, Jhelum and Kabul basins. Also, the cross-correlations are generally better for higher SSI accumulation periods.

Fig. 9
figure 9

Heatmaps of monthly cross-correlations between SSI and SPEI for the four Upper Indus catchments, i.e., Chenab, Jhelum, Indus and Kabul (depicted in rows). Each column corresponds to a different SSI accumulation period (i.e., SSI-1, SSI-3 and SSI-6). SSI-3 and SPEI cross-correlations are computed with a lag of two months to avoid double accumulation of precipitation.

The cross-correlations between SPEI and SSI are not strong for (i) shorter SSI accumulation periods in summer and fall months (June-October) for all catchments and (ii) the Upper Indus. Plausible reasons for the poor cross-correlations observed in the Indus catchment are explored in Sect. 4.4. The relatively weak cross-correlations for summer and fall months could be attributed to errors in the gridded CRU precipitation dataset in the summer months. Summer precipitation and streamflows in the Upper Indus catchments are strongly influenced by monsoon rains, which are not adequately represented in most historical climate data sets64.

Lagged cross-correlations

The strong cross-correlation observed between SPEI and SSI, for Chenab, Jhelum and Kabul (especially in late winter and spring months) is a valuable outcome of our analysis, that can be used effectively in Drought Early Warning Systems (DEWS) in Pakistan. For instance, if lagged cross-correlations exist between SPEI and SSI (where SPEI values are lagged), SPEI values could be used in forecasting and early detection of droughts. Figure 10 (for SSI-3 and SSI-6) report lagged cross-correlations between SPEI (with different accumulation periods) and SSI, for Chenab, Jhelum, Indus and Kabul basins.

Fig. 10
figure 10

Heatmaps of lagged (1–5 months) monthly cross-correlations between (a) SSI-3 and SPEI and (b) SSI-6 and SPEI for the four Upper Indus catchments, i.e., Chenab, Jhelum, Indus and Kabul (depicted in rows). Each column corresponds to a different lag (in months; increasing from left to right).

When cross-correlations are lagged, strong correlation values (\(\:r>0.7\)) are found between lagged SPEI and SSI-3 from March to June for Chenab, from March to August for Jhelum, and from April to September for Kabul. Hence, SPEI values can be used to provide early warning of surface outflow deficit from these basins in the above-mentioned. Moreover, the information embedded in the lagged cross-correlations can be valuable in (i) seasonal drought planning pertaining to surface water use in the lower Indus basin21 and (ii) seasonal planning of operations of the Mangla reservoir (in Jhelum) and operations of the Indus Basin Irrigation System (IBIS)15,65. For instance, early prediction/detection of surface flow deficit in the Upper Indus catchments in April to June, also referred to as the early Kharif season, is extremely valuable for water management and operations in the lower Indus basin, including operations of storage reservoirs50,51,65. SPEI values linked to Chenab, Jhelum and Kabul basins could thus be used in both seasonal drought management planning and water systems’ operations.

Spatial validation

Figure 11 provides a spatial validation of SPEI’s capability of adequately representing drought events in the Upper catchments of the Indus basin, especially in the Upper Chenab, Upper Jhelum, Upper Indus and Kabul sub-basins. Figure 11 reports gridded SPEI-12 values (obtained from66 of the Chenab, Jhelum, Indus and Kabul sub-basins of the Upper Indus basin. It is evident from Fig. 11 that the drought condition in June 2001 (accumulated over 12 months) is spatially prevalent in all sub-basins except Indus. Within the Indus sub-basin, a significant part of the catchment depicts a high 12-month water surplus in June 2001, according to SPEI. This anomaly, reported in the upper portion of the catchment, is primarily due to inaccurate representation of precipitation (by the CRU data set) in this part of the catchment. Numerous prior studies report that, primarily due to the lack of meteorological stations in the upper portion of the Upper Indus basin, gridded precipitation data (of all datasets tested) is highly inaccurate in this region of the catchment36,64. Figure S3 in the supplement is also consistent with this observation, where the extreme flood event of August 2010 is not represented in the upper region of the Indus sub-basin.

Fig. 11
figure 11

Spatial representation of the Droughts of 1998-02 in the Upper Indus sub-basins of Pakistan, i.e., (a) Chenab, (b) Jhelum, (c) Indus and (d) Kabul, using gridded SPEI-12 values of June 2001. Gridded SPEI values are obtained from the SPEIbase data set45,66.

We believe that the poor cross-correlations observed between SSI and SPEI for the Indus catchment in the Upper Indus basin are partially due to inaccuracies in gridded precipitation data in the upper regions of the catchment36,64 suggest mechanisms for correcting/improving precipitation records in these catchments. These improvements could be used in future studies to improve cross-correlations between SSI and SPEI in this catchment. In a hydrological modelling study51, split the Upper Indus catchment into further two sub-catchments to significantly improve modelling results for Indus at the Tarbela. A similar split may also be experimented for the SPEI-SSI cross-correlation analysis, in order to develop an effective DEWS for the Upper Indus catchment in future.

Conclusions

We investigated the combined use of the Standard Precipitation Evaporation Index (SPEI) and the Standard Streamflow Index (SSI), to analyze their capability in monitoring and early detection of droughts in the four key upper catchments of the Indus Basin of Pakistan, i.e., Chenab, Jhelum, Indus and Kabul.

Our combined indicator-based drought analysis shows that both SPEI and SSI are able to identify historical hydrological droughts and streamflow deficits at the outlets of these catchments. Moreover, a brief analysis of two key drought characteristics, duration and severity, shows that there is a high correlation between drought duration and severity. Moreover, empirical distribution plots of drought duration were also analyzed and visual coherence between SPEI-based and SSI-based drought duration statistics (and respective empirical distributions) was also observed (see Fig. S1).

Since a key purpose of this study was to analyze coherences between SPEI and SSI, in order to unearth new insights for improvement of existing DEWS, we also analyzed cross-correlations and lagged cross-relations between SPEI and SSI. When seasonality was not considered, weak cross-correlations were observed between SPEI and SSI for all catchments (see Fig. 8). In order to incorporate seasonality into our analysis, we also computed monthly cross-correlations between SPEI and SSI. Strong cross-correlations (i.e., \(\:r>0.7\)) were observed for Chenab, Kabul and Jhelum catchments, especially in late winter and early spring months (see Fig. 9).

Our lagged monthly cross-correlation analysis showed strong correlations between SPEI and SSI in (at-least) the early Kharif season (April to June)50, for the Chenab, Jhelum and Kabul catchments (with lags of up to two months). This is a very important outcome, since monitoring and prediction of streamflow deficits in these basins in early Kharif months is extremely critical for reservoir operations and irrigation planning in the lower Indus, and for seasonal drought forecasting and planning in the country15,21,50. Given these strong lagged cross-correlations, we believe that SPEI could be used in operational drought forecasting and warning systems within Pakistan in future, for early streamflow deficit prediction in Kabul, Chenab and Jhelum basins in the early Kharif season.

Our cross-correlation analysis also showed poor correlations between SPEI and SSI for the Indus sub-catchment of the Upper Indus. This is primarily due to inaccuracies in gridded precipitation data in upper regions of this catchment36. In future, datasets other than CRU39, including corrected precipitation datasets developed for this catchment36,64 could be tested to improve SPEI and SSI correlations in the Indus sub-catchment. Another future research direction is a more comprehensive multivariate distribution-based analysis of drought characteristics (e.g., severity, duration and frequency) and associated risks for the Upper Indus catchments using both SSI and SPEI63.