1 Introduction

Drought is an environmental calamity that can lead to a lack of drinking water for humans and livestock, reduced crop production, and increased economic costs, with considerable damaging effects on humankind (Dracup et al. 1980; Wilhite 2000). Usually, drought can be classified into four different types: meteorological, agricultural, hydrological, and socio-economic (Van Loon 2015; Wilhite and Glantz 1985). The first three types of droughts are caused by a deficiency of precipitation, soil moisture, and water resources (surface or groundwater), respectively (Mishra and Singh 2010; Zhu et al. 2016). Instead, the socio-economic drought represents adverse conditions in society, economy, and environment when water supply cannot meet water demand (Zhao et al. 2019; Zseleczky and Yosef 2014). Meteorological and hydrological droughts need to be analyzed to evaluate the impacts of the drought phenomenon on water resources. When meteorological drought conditions persist, soil moisture and surface water levels can decrease, eventually leading to hydrological drought (Van Loon 2015). However, the time it takes for meteorological drought to turn into hydrological drought can vary significantly depending on certain conditions. For example, if an area has high evaporation and high water demand, it may only take a short time; on the other hand, if evaporation and water demand are low in an area, it may take longer (Bakke et al. 2020; Zhang et al. 2022). The Standardized Precipitation Index (SPI) and the Standardized Runoff Index (SRI) are the most extensively employed indicators in literature for monitoring meteorological and hydrological drought, respectively (Katipoğlu and Acar 2022; Katipoğlu et al. 2021; McKee et al. 1993; Shukla and Wood 2008). Based on drought indicators, analyzing the relationship between different types of drought is one of the main approaches to understanding how drought spreads. Many researchers have evaluated the relationship between drought types. For example, Huang et al. (2015) employed the heuristic segmentation method to analyze hydrological and meteorological droughts in the Wei River Basin in China. Moreover, the cross-wavelet analysis has been used to detect the relationship between agricultural drought based on the Palmer Drought Severity Index (PDSI) and meteorological drought based on SPI. Results evidenced that the variability of the lag time of the first in response to the latter was significant on a seasonal scale: rapid response in summer and relatively slow in autumn. It was also revealed that vegetation influenced the transition from meteorological to agricultural drought. Wang et al. (2019) modeled the relationship between these two types of drought in the Yellow River Basin.

Revealing the link between meteorological and hydrological drought is essential to provide early warnings of the latter. The run theory method, the correlation analysis method, the signal process, and the nonlinear response methods have all been employed in literature to reveal the lag time between the two different types of drought (Lorenzo-Lacruz et al. 2013; Wu et al. 2017; Zhao et al. 2014). For example, to detect this time interval Zhao et al. (2014) applied the run theory in the Jinghe River Basin (northern China). Wu et al. (2016) analyzed the response of hydrological drought to meteorological drought in China’s Jinjiang River Basin, evaluating an average delay time of 3 months in winter and 1 month in the other seasons. Tijdeman et al. (2018) assessed the relationship between the two types of drought, considering natural and anthropogenic factors in the contiguous United States. Wu et al. (2021) compared run theory, correlation analysis, and nonlinear response methods to the analysis of the propagation of meteorological and hydrological drought in the Dongjiang River Basin, southern China, evidencing that the run theory is an essential tool for preventing and mitigating drought. Um et al. (2022) examined the extent of all types of droughts in the Yangtze River basin (China). As a result, the spread of meteorological to hydrological droughts took a relatively short time (1–2 months), while the spread of meteorological to agricultural droughts took longer (1–7 months). Considering meteorological and hydrological drought, He et al. (2022) applied a correlation analysis for modeling the drought transition time from the first type to the second in the Baiyangdian Lake in China. They evidenced that, although meteorological droughts have a high frequency, their duration is short, thus proving that meteorological drought does not always generate hydrological drought. Sarwar et al. (2022) employed the streamflow drought index (SDI) and SPI values, and by means of the Mann–Kendall test, a regression analysis, and a moving average approach revealed that there is a significant relationship between hydrological and meteorological droughts. Li et al. (2022) applied the Pearson correlation, the Cross-wavelet power spectra and the wavelet coherence and showed that the average duration and size of both meteorological and hydrological droughts increase with the increase of the drought time-period.

The statistical correlation analysis between hydrological and meteorological drought is performed by many researchers (Barker et al. 2016; Kazemzadeh and Malekian 2016; Katipoğlu et al. 2020). Analyzing the transition between different types of drought is important for water resources planning and management (Maity et al. 2016). As hydrological drought preparedness is essential to make policy decisions related to water and agricultural sectors at the river basin scale, it is highly essential to study the transition from meteorological drought to hydrological drought. In order to capture such a transition, running (dynamic) correlation analysis may be helpful. However, as many of the hydro–meteorological datasets including derived indices, such as the drought index, possess multiscale features, the correlation analysis performed by means of the Pearson correlation alone is insufficient. Therefore, associations between the drought indices need to be investigated, for which a multiscale time–frequency transformation technique is considered to be the most appropriate tool. Continuous wavelet transforms, and their wavelet coherence features are popular in hydro-climatic teleconnections, which were also applied earlier to investigate the association between predictor-predictant datasets in geophysical and hydro-meteorological studies (Su et al. 2019; Das et al. 2020; Song et al. 2020; Sreedevi et al. 2022). About 2 decades back, the difficulties in the implementation of wavelet were transformed into a powerful mathematical transform, namely Hilbert–Huang Transform (HHT) (Huang et al. 1998), which involves a data-adaptive decomposition process called Empirical Mode Decomposition (EMD) followed by Hilbert Transform (HT). The decomposition of signals using the traditional EMD proposed was found to be inefficient because of its inability to separate the components based on distinct frequency characteristics. The presence of multiple frequencies in signal components often called as ‘mode mixing’ was well addressed in the past using noise-assisted ensemble variants like ensemble EMD (EEMD) and complete ensemble EEMD with adaptive noise (CEEMDAN) (Wu and Huang 2009; Torres et al. 2011). However, researchers later on found that even the CEEMDAN variant is contaminated by residual noise and posses challenges to the time–frequency analysis. As a solution to this issue, Colominas et al. (2014) provided an improved version of CEEMDAN (ICEEMDAN) which is reported to be efficient in minimizing the noise residual content (Plocoste et al. 2022, 2023). Moreover, the advancements of HHT have led to a powerful multiscale dynamic (running) correlation method, namely Time Dependent Intrinsic Correlation (TDIC). TDIC, propounded by Chen et al. (2010), offers cleaner and separable correlation features in the time–frequency domain by preserving the stationarity property during the correlation estimation (Adarsh and Janga Reddy 2016). Therefore, integrating ICEEMDAN with TDIC can efficiently run correlation analysis in a multiscale standpoint. The TDIC method was applied for multiscale hydro-climatic teleconnection studies (Adarsh and Janga Reddy 2016, 2018, 2019). Johny et al. (2020) performed the teleconnection analysis of drought, but the study focused on the link between drought index with large-scale climatic oscillations. The application potential of TDIC for analyzing the multiscale correlation between meteorological and hydrological drought characteristics has not been explored yet. According to our knowledge, no dynamic (running) correlation studies between different types of drought have been reported in the multiscale domain. The outputs of such a study are thought to help water managers and policymakers to develop effective strategies to manage and mitigate the effects of drought and assess the impact of climate change on global warming.

This study aims to evaluate the characterization and change of two types of drought, i.e., meteorological and hydrological drought, and how they relate over different time scales. The specific objectives of this paper include (i) decomposition of SPI and SRI drought series of different scales pertaining to five different hydrometric and rain-gauge stations of the Wadi Mina basin using the robust ICEEMDAN method; (ii) performing the TDIC analysis to investigate the dynamic association between the selected types of droughts. The next section of the paper presents the study area and data statistics. Then, the details of different algorithms are presented. The Results and Discussion section describes the outcomes of ICEEMDAN and TDIC analysis of SPI and SRI of the basin. Important conclusions drawn from the study are presented in the last section.

2 Study area and data

The Wadi Mina basin (00° 22′–01° 09′ E and 34° 41′–35° 35′ N) is located in the northwest of Algeria and has a total area of 4900 km2 (Fig. 1) and an altitude ranging from 164 to 1327 m (Achite et al. 2022). The Wadi Mina involves four major tributaries: Wadi Mina, Wadi Haddad, Wadi Abd and Wadi Taht. Following the Köppen–Geiger classification (Köppen 1936), the climate of the Wadi Mina basin can be considered a Mediterranean climate (classes Csa and Csb), with hot summers and cold winters when the highest percentage of the mean annual rainfall (approximately 260 mm) occurs. As regards temperatures, they range between 16 and 19.5 °C. Nearly half of the basin is covered by vegetation of varying densities, with scrubs (32%) and forest and cereal crops (35.8%) representing the highest percentages.

Fig. 1
figure 1

Study area, pluviometric and hydrometric network

Full size image

The National Agency of Water Resources has collected hydro-meteorological data in the Wadi Mina basin. In particular, monthly rainfall and runoff series collected between 1973 and 2012 in five rain gauges and five hydrometric stations have been considered in this work (Fig. 1 and Table 1).

Table 1 Rainfall and hydrometric stations characteristics
Full size table

3 Methodology

3.1 Standardized Precipitation Index (SPI)

According to McKee et al. (1993), given the precipitation amount recorded at a station, the SPI allows us to evaluate wet and dry spells and their characteristics (e.g., duration, severity, intensity). The SPI is one of the most used and well-known indexes and has been described in depth in several studies (see, e.g. Achite et al. 2023). The method of drought estimation involves (i) the computation of aggregated monthly precipitation series for a specific aggregation scale, k; (ii) fitting of a cumulative distribution function (CDF) using appropriate probability distribution (usually gamma distribution function) on aggregated monthly precipitation series (say k = 3, 6 and 12 months); (iii) equipropability transformation of fitted CDF to CDF of standard normal distribution with zero mean and unit variance, to get the SPI. The details of distribution fittings are given below.

The probability density function (PDF) of Gamma function is evaluated as:

$$g\left( x \right) = \frac{1}{{\beta^{\alpha } {\Gamma }\left( \alpha \right)}}x^{\alpha - 1} e^{ - x/\beta } \;{\text{for}}\;x > 0\;{\text{and}}\;{\text{with}}\;\Gamma \left( \alpha \right) = \mathop \smallint \limits_{0}^{\infty } y^{\alpha - 1} e^{ - y} {\text{dy}}$$
(1)

in which x is the precipitation amount and α and β are the shape and the scale parameters which have been estimated by applying the of Thom (1958). Γ(α) is the gamma function that is undefined for a null rainfall amount, and thus a modified CDF must be considered.

$$H\left( x \right) = q + \left( {1 - q} \right) G\left( x \right),$$
(2)

where q is the probability of null precipitation and G(x) is the CDF:

$$G\left( x \right) = \mathop \smallint \limits_{0}^{x} g\left( x \right){\text{dx}} = \frac{1}{{\beta^{\alpha } \Gamma \left( \alpha \right)}}\mathop \smallint \limits_{0}^{x} x^{\alpha - 1} e^{ - x/\beta } {\text{dx}},$$
(3)

H(x) is then transformed into the standard normal random variable Z, with a mean of 0 and a variance of 1, which is the value of the SPI following Abramowitz and Stegun (1965)

$${\text{Z = SPI}} = \left\{ {\begin{array}{*{20}l} { - \left( {t - \frac{{c_{0} + c_{1} t + c_{2} t^{2} }}{{1 + d_{1} t + d_{2} t^{2} + d_{3} t^{3} }}} \right)} \hfill & {{\text{0 < H}}\left( x \right) \le 0.5\;with\;t = \sqrt {{\text{ln}}\left( {\frac{1}{{\left( {H\left( x \right)} \right)^{2} }}} \right)} } \hfill \\ { + \left( {t - \frac{{c_{0} + c_{1} t + c_{2} t^{2} }}{{1 + d_{1} t + d_{2} t^{2} + d_{3} t^{3} }}} \right)} \hfill & {0.{\text{5 < H}}\left( x \right) \le 1.0\;with\;t = \sqrt {{\text{ln}}\left( {\frac{1}{{\left( {1 - H\left( x \right)} \right)^{2} }}} \right)} } \hfill \\ \end{array} ,} \right.$$
(4)

with ci and di coefficients. Table 2 shows the SPI classification proposed by McKee et al. (1993).

Table 2 Drought classification based on SPI/SRI value and corresponding event probabilities (Mckee et al. 1993)
Full size table

The aggregation scale (k) is significant when evaluating the effects of drought. For example, an SPI with a monthly scale (i.e., SPI1) indicates the monthly fluctuation of wet and dry period relative to its mean. Longer time scales, such as SPI3, SPI6, and SPI12, display precipitation variations over the last 3, 6, and 12 months, which may be used as indicators of prolonged wet or dry periods (Thomas et al. 2015).

3.2 The Standardized Runoff Index (SRI)

SRI was developed by Shukla and Wood (2008) to assess hydrological drought considering streamflow data. The concept employed by McKee et al. (1993) for the SPI has been applied in defining the SRI. It involves fitting of suitable distribution to flow records of a particular location. After this, PDF is calculated, and it is transformed into a standardized Gaussian distribution with mean zero and unit variation that gives SRI.

The classification of SRI values is done similarly to the SPI in Table 2.

3.3 Multiscale running correlation analysis

In order to capture the association between SPI and SRI, a novel multiscale running correlation analysis is proposed. This proposed framework is developed based on Hilbert Huang Transform (Huang et al. 1998). In general, the procedure involves, (i) the decomposition of the two correlated series (here SPI and SRI) to a multitude of components of distinct periodic properties called Intrinsic Mode Functions (IMFs) and a final residue using EMD or its variants, (ii) determining the instantaneous frequencies of the IMFs; (iii) execute a dynamic correlation analysis between IMFs; (iv) compare the analysis of the low-frequency components of the two signals. The general framework followed in the research is expressed in Fig. 2. Steps (ii) and (iii) together form the steps of running correlation analysis, namely TDIC.

Fig. 2
figure 2

Framework for multiscale correlation analysis

Full size image

In short, the proposed involves a data decomposition phase and running correlation phase. EMD is a popular data adaptive signal decomposition method, which decomposes a given time series x(t) into a certain number of mean zero components and a final residue in the following form

$$x(t) = \sum\limits_{i = 1}^{N} {IMF_{i} (t) + R(t)}$$
(5)

where R is the final residue and N is the number of IMFs.

The IMFs are zero mean series, extracted through an iterative process by identifying the instantaneous peak or trough points. The process involves developing curves connecting the extremes via interpolation function, computation of the mean, differencing the signal from the mean etc., that are repeated until the differenced signal becomes a zero-mean series to get the first mode. Then, the repetition of the reminder signal (R1(t) = x(t) − IMF1(t)) is used as a new signal until it becomes a truly monotonic or single peaked signal. The signals and its component ideally should maintain the orthogonality property, which implies the quality of decomposition. It is computed using the Index of Orthogonality (IO) proposed by Huang et al. (1998) as

$$IO = \sum\limits_{t = 0}^{T} {\left( {\sum\limits_{j = 1}^{N + 1} {\sum\limits_{k = 1}^{N + 1} {\frac{{C_{j} (t)C_{k} (t)}}{{x^{2} (t)}}} } } \right)}$$
(6)

in which the square of the signal (denominator) is computed as

$$x^{2} (t) = \sum\limits_{j = 1}^{N + 1} {C_{j}^{2} (t)} + 2\sum\limits_{j = 1}^{N + 1} {\sum\limits_{k = 1}^{N + 1} {C_{j} (t)C_{k} (t)} }$$
(7)

C represent component including residue, N is the number of IMFs.

The complete process of EMD is well described in the literature (Huang et al. 1998). However, during the decomposition of a signal into IMFs, it can sometimes happen that signals with similar frequencies are not assigned to separate IMFs. This problem is described as ‘mode mixing’ in literature (Wu and Huang 2009; Olvera-Guerrero et al. 2017).

After all, each IMF consists of several oscillation modes, which causes them to lose their physical meaning. Noise-assisted data analysis methods have been proposed to eliminate this problem (Wu and Huang 2009; Torres et al. 2011; Colominas et al. 2014). To avoid the problem of mode mixing, this study followed the Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN) propounded by Colominas et al. (2014) for the decomposition of SPI and SRI time series, and multiscale correlation analysis of IMFs are performed subsequently, using the Time Dependent Intrinsic Correlation (TDIC). Hence, the ICEEMDAN–TDIC approach is used to analyze the dynamic relationships between SPI and SRI across different time scales.

3.3.1 ICEEMDAN

The Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) is a noise-assisted multiscale method based on EMD. In the case of CEEMDAN, some residual noise will be preserved and the signal may provide some spurious modes in the early stages. Colominas et al. (2014) suggested an Improved CEEMDAN (ICEEMDAN), which is claimed to give the modes lesser noise contamination andmore physical meaning.

The main steps involved in ICEEMDAN are (Colominas et al. 2014):

  • Computation of local means of \(I\) realizations \(s^{\left( i \right)} = s + \beta_{0} E_{1} \left( {w^{\left( i \right)} } \right)\) by EMD to get the first residue:

    $$R_{1} = M\left( {s^{\left( i \right)} } \right)$$
    (8)

    Here s is the signal (SPI or SRI in this study), w is zero mean unite variance white Gaussian noise and \(\beta_{0} = \frac{{_{0} std \left( s \right)}}{{std\left[ {E_{1} (w^{\left( i \right)} } \right)]}}\), where \(\varepsilon_{0}\) is the noise amplitude and std stands for standard deviation of the first mode of white noise series obtained by EMD and M stands for local mean.

  • Computation of the first mode (with k = 1):

    $$IMF_{1} \widetilde{ = s} - R_{1}$$
    (9)

  • Estimation of the second residue as the average of local means of the realizations \(R_{1} + \beta_{1} E_{2} \left( {w^{\left( i \right)} } \right)\)

    $$\widetilde{{IMF_{2} }} = R_{1} - R_{2} = R_{1} - M\left( {R_{1} + \beta_{1} E_{2} \left( {w^{\left( i \right)} } \right)} \right)$$
    (10)

  • Computation of the k-th residue for \(k = 3, . . . , K\),

    $$R_{k} = M\left( {R_{{\left( {k - 1} \right)}} + \beta_{{\left( {k - 1} \right)}} E_{k} \left( {w^{\left( i \right)} } \right)} \right)$$
    (11)

  • Computation of k-th mode iteratively

    $$\widetilde{IMF}_{k} = R_{{\left( {k - 1} \right)}} - R_{k}$$
    (12)

    Here, E1(.), E2 (.) and Ek(.) stands for the operators for EMD.

3.4 Time-dependent intrinsic correlation (TDIC)

TDIC is a dynamic correlation analysis method that first determines the instantaneous frequencies (IFs) of the IMF components obtained by ICEEMDAN (or other variants of EMD). To determine the instantaneous frequencies (IFs), the methods of Hilbert transform or its variants like Normalized Hilbert Transform-Direct Quadrature (NHT-DQ), Zero crossing (ZC), Generalized Zero Crossing (GZC), Teager Energy Operator (TEO) etc. can be used (Huang et al. 2009). The superiority of NHT-DQ method over traditional HT was well documented by one of the co-authors in Adarsh and Janga Reddy (2019). The method is a combination of a normalization scheme of the resulting IMFs, and the direct quadrature (DQ) method, preserving the mathematical correctness and the physical meaning of the transformation. In this study, the NHT-DQ scheme is used to find the IFs required to perform the TDIC analysis. The detailed documentation on NHT-DQ scheme is available in Adarsh and Janga Reddy (2019). Then, the correlation estimation is carried out in a sliding window framework, in which the sliding window size is fixed adaptively based on the maximum number of instantaneous periods' IPs. This ensures stationarity property within each local sliding window, and enables the use of traditional correlation between the windows from two signals. As the sliding windows are fixed adaptively based on IPs, the computed correlation coefficients will capture the non-stationarity of the original physical processes (Chen et al. 2010).The steps of ICEEMDAN–TDIC method involves the following steps:

  1. (1)

    Decomposition of the two associated time series using ICEEMDAN

  2. (2)

    Application of HT on the IMFs using the NHT-DQ methodto get the IF and (and hence instantaneous periods, IPs)

  3. (3)

    Minimum sliding window size (\(t_{d} )\) for the local correlation estimation is computed as \(t_{d}\) = max (T1i (\(t_{k}\)), T2i (\(t_{k}\))) where \(T{}_{1,i}\) and \(T{}_{2,i}\) are the IPs of the pairs of IMFs at specific scale. In this study IMFs of SPI and SRI are the pairs;i is the number of IMFs varying from 1 to N; 1 and 2 are representing SPI and SRI series; k is indicating any instant for which the sliding window is computed.

  4. (4)

    Finding the size of the sliding window (SSW) as \(t_{w}^{n}\) = [\(t_{k}\) − n \(t_{d}\)/2:\(t_{k}\) + n \(t_{d}\)/2] for a specific IMF of the first signal (say SPI3), and \(t_{w}^{n}\),τ = [\(t_{k}\) − n \(t_{d}\)/2: \(t_{k}\) + n \(t_{d}\)/2] for the corresponding IMF of the second signal (say SRI3), where n is any positive number, normally is selected as 1 (Huang and Schmitt 2014).

  5. (5)

    The running correlation between the IMFpairs is computed as

\(R_{i} (t_{k}^{n} ) = Corr(IMF_{1i} (t_{w}^{n} ),IMF_{2i} (t_{w}^{n} ))\) at any time tk.

Corr is the general correlation coefficient between the two time series(here the windows from two IMF pairs at any time tk)

The statistical significance of correlations was determined by Student’s t-test. This can be repeated until the sliding window's boundary exceeds the time series' endpoints.

4 Results and discussion

4.1 Meteorological and hydrological drought analysis

The temporal variability analysis of SPI and SRI values at five stations in the study area are presented in Figure S1 and Figure S2, respectively. As regards the 1-month SPI, Extreme Drought (ED) values were noted at 3 out of 5 stations in October 1983. In the same month, all the stations showed SPI values lower than − 2 for the 3-month SPI. For the same time aggregation, ED values have also been evaluated at 4 out of 5 stations in May 1987 and in March and April 2000, while 3 out of 5 stations showed these values in November and December 1981 and December 1989. The 6-month SPI showed ED values, especially in December 1981 (3 out 5 stations), June, September, and October 1983 (5 out 5 stations), March 1997 (5 out 5 stations), and June and July 2000 (4 out 5 stations). As regards the 9-month time scale, SPI values lower than − 2 have been mainly detected between September and December 1983 for all the stations for 4 consecutive months and in September 2000 in 4 out 5 stations. A similar behavior was identified for the 12-month-SPI, between November 1983 and January 1984, with all the stations showing ED values, while the 24-month-SPI didn’t show any particular behavior (Figure S1).

Concerning the SRI, ED values have rarely been reached considering the 1-month timescale, while in 2 out of 5 stations, these values were identified for the 3-month SRI in March 2004 and March 2008. Regarding the 6-month SRI, values lower than − 2 have been detected in 3 out of 5 stations in July 2000 and 2 out of 5 stations in December 2006 and July 2008, considering the 9-month SRI. Finally, ED values in more than one station were identified in January and February 2007 for the 12-month SRI and in almost 16 consecutive months, between December 2005 and March 2007 for the 24-month SRI (Figure S2).

Table 3 notes the positive and negative peak values, the critical period duration, and the starting and ending dates for all rain gauge stations. The longest period in which SPI and SRI values remain negative without interruption is called the critical period. In rain gauge station S1, the index values were less than − 2 and greater than + 2 in all time scales in the most severe wet and drought periods, respectively. The drought duration increased with the increase in the time scale for station S1. Moreover, the critical periods at this station usually occurred just before and after 2000, and only for SPI3 was the critical period detected in 1983. Looking at the results of this station, in general, it is clear that the effect of drought has been felt more recently in the region than in past times, and has been experienced longer. Generally, in station S2, the drought lasted shorter than in station S1, and the most critical period mainly occurred between 1980 and 1985. The critical period duration, start and end dates in station S3 are similar to the ones evaluated for station S2.In addition, in this station, the duration values increased with the increase of the time scale. The results of station S4 are similar to those obtained in station S1, with the critical period mainly identified in the 2000s, except for SPI1 and SPI3. At station S5, the most critical period was between 1982 and 1984. In this station, the duration value changed between 7 and 61 months as the time scale increased from 1 to 48 months. As a general comment, ED and EW values occurred in all the regions, with positive peak values higher than + 2 detected in all the stations and for all the time scales. In contrast, stations S2 (SPI1) and S3 (SPI1 and SPI24) did not show negative peak values lower than − 2.

Table 3 The positive and negative peak values, duration, date of start, and end of critical period details for the different stations
Full size table

Moreover, Table 3 shows the negative and positive peak values obtained at different time scales at the hydrometric stations, the start and end dates of the critical period, and the duration values. The driest period at station HS1 mainly started in 2003 and 2004, except for the 3-month time scale, and ended in 2007–2008. Absolute negative peak values obtained at this station are larger than positive peak values except for SRI1. The critical period at the station HS2 started in the late 1990s and continued until the mid-2000s. It can be clearly stated that the critical period occurred in the recent past at both of these stations. Critical periods were experienced at stations HS3, HS4, and HS5 in the 1980s. According to the average duration values obtained for all the timescales, the average hydrological drought lasted 39 months for station HS1, 54 months for station HS2, 35 months for station HS3, 38 months for station HS4, and 89 months for station HS5. The average duration value for all regional stations was calculated as 15.6, 32.2, 48, 62, and 69 months and 77.8 months for SRI1, SRI3, SRI6, SRI9, SRI12, and SRI24, respectively. From the evaluation of the duration values obtained at the different stations, it is understood that the duration value generally increases with the increase of the time scale. Among the stations examined, the most critical station in terms of hydrological drought is the HS5 station, according to the duration values. It is suggested that the results of the drought analysis at HS1 and HS2 stations, which have been determined to be quite dry in recent years, are remarkable, and the decrease in the planning of the activities carried out depending on the streamflow values in this region should be taken into account.

4.2 Multiscale correlation between SPI and SRI

In this study, three kinds of correlation analyses between SPI and SRI are performed: (i) correlation between SRI of a specific hydrometric station, with SRI of all other stations, in order to assess the most correlated pairs, considering the complete length of time series; (ii) for the most correlated SPI–SRI pairs, the correlation between components (IMFs and Residue), to assess the strength of correlation, considering the complete length of time series, which justifies the necessity of TDIC analysis very well; (iii) running correlation between IMF pairs, for in-depth investigation between SPI and SRI. Thus, before proceeding with the multiscale correlation between SPI and SRI, it is worth investigating the stations in which the SPI and SRI shows the highest correlation. The correlation coefficients between the SPI and SRI values obtained at different rainfall and hydrological stations at 3-, 6- and 12-month time scales are shown in Fig. 3. The analysis showed that the SRIs of station 1 are correlated with the SPIof station 4; the SRI of station 2 is correlated with the SPI of station 1; SRIs of stations 3, 4, and 5 are correlated most with the SPIs of station 5, irrespective of the scale and hence the type of drought (like medium- or long-term represented by 3-month, 6-month and 12-month).

Fig. 3
figure 3

Correlation coefficients between the SPI and SRI indexes for different periods

Full size image

We have considered the SPI3, SPI6, and SPI12 as the representatives of short- medium- and long-term droughts (Thomas et al. 2015), and ICEEMDANis invoked to find the components of SPI and SRI at different time scales. Both SRI and SPI time series are decomposed using ICEEMDAN with a noise amplitude of 0.2 and ensemble number of 200 (number of realizations) following experimentations and past recommendations (Colominas et al. 2014; Plocoste et al. 2022, 2023). For the ICEEMDAN-based decomposition, we ensured the index of orthogonality defined by Huang et al. (1998) was found to be less than 10−4 in each case, which ensured the credibility of the decomposition process. The decomposition resulted in different numbers of modes in different cases, at the most, they differ by one in the SPI and SRI time series of a station. The periodicity of each mode is computed by the zero-crossing procedure stated by Huang et al. (2009), which is found to be closely matching the period computed from IFs estimated by the NHT-DQ method. The third mode shows annual periodicity in the case of the 3-month drought index series. The periodicity of SPI and SRI times series and their variability (i.e., a ratio of a mode variance to the variance of a signal expressed in %) is shown in Table 4 which indicates that the dominant variance contribution is at different scales for different SPI series. Similar results of other stations are presented in Tables of Supplementary Information (SI) file. Furthermore, the mean periods of all the modes of all SPI and SRI series of different stations are presented in Fig. 4. Figure 4 and Table 4 confirm that the periodicity of different modes not necessarily follows dyadic property (Huang et al. 1998). These periods could be connected with the processes governing the mechanism behind drought and streamflow variability. The local scale meteorological processes and global climatic oscillationsare well proven to be modulating the rainfall and hence the streamflow (Massei et al. 2007; Adarsh and Janga Reddy 2016; Shamna et al. 2022; Johny et al. 2020; Mohan et al. 2023). The climatic oscillations like El Nino Southern Oscillations (ENSO), Pacific Decadal Oscillations (PDO), Mediterranian Oscillations (MO) and North Atlantic Oscillations (NAO) are deciphered at diverse periodicity ranging from interannual to interdecadal scales. In the Algerian and Tunisian river basins also these teleconnections with climatic oscillations are well proven and are found to be influencing the rainfall processes,streamflowand drought (Ouachani et al. 2013; Meddi et al. 2014; Hallouz et al. 2020).

Table 4 Periodicity (months) and variability (V in %) of different modes of SRI of station 1 and SPI of station 4. LT refers long term
Full size table
Fig. 4
figure 4

Mean time periods of IMFs of SPI and SRI series of different scales. Upper panels show the short term drought, middle panels shows the results of moderate droughts and lower panels show the results of long term drought

Full size image

The correlation between the modes of SRI of station 1 and the modes of SPI of station 4 is shown in Table 5, which points out that the correlation between SRI mode and SPI mode is very weak except for a few low-frequency modes and residue. The cross-correlation values for other stations are provided in Supplementary Information (SI) file. It is worth mentioning that such a correlation is computed by considering the complete length of the series. However, it must be noted that the use of simple statistical correlation like Pearson correlation can capture only the linear association between two time series. Hydro-meteorological processes are usually non-stationary and nonlinear processes, for which the linear correlation between the two data series on the whole domain cannot reveal the possible relationship between them. The drought indices like SPI and SRI are governed by the processes of different periodic features, indicating a multiscale and nonlinear problem. In such a context, Pearson correlation alone cannot capture the true associations, as the SPI–SRI associations may be of different nature and different strengths across the scales and along the domain. In short, the associations SPI–SRI at different time scales and also, along the time domain need to be performed.

Table 5 Correlation between modes of SRI of hydrometric station 1 and modes of SPI of raingaugestation 4. Significant correlations (at 5% level) are underlined
Full size table

Since streamflow drought and meteorological drought are temporal phenomena, the evolution of the drought along the time spell needs to be investigated. This brings forth the notion that for a specific time scale, in some of the spells, the correlation may be high and negative (positive), but in some other spells, the correlation may be high and positive (negative), which brings down the overall correlation to a very low value on considering the complete time span of the study. A running correlation between SPI and SRI is to be performed to track such localized relationships between SPI and SRI across the scales. TDIC analysis is one such tool which can capture the associations at different distinct scales, as a decomposition based on EMD (or its variants) is involved. Also, as it performs the correlations along the time domain in a running correlation mode, the strength of association along the domain can be captured.

Thus the TDIC analysis stemmed from HHT, which uses a periodicity-based sliding window selection and ensures stationarity property within each time window, which is a more credible approach in running correlation studies.

Before proceeding with the TDIC analysis, comparison plots of IMFs of SPI and SRI are developed to understand the pattern of evolution of modes in terms of phase and phase difference. The plots are provided in Figs. 5, 6 and 7. Figure 5 shows that phase relationships are not consistent along the time domain and scale. When considering the case of short-term drought, in-phase relationships are in-phase at the annual scale of IMF5. Moreover, with localized spells lasting up to a decade, such relationships persist at annual scales. For example, for HS4, during 1990–2000, in-phase relationships emerge between SPI3 and SRI3. SRI leads over SPI (or vice versa) is visible mostly at inter-annual time scales. Say, during 1990s in HS2 at IMF4 and IMF5, during 1980s in IMF4 of HS1, etc. For moderate drought conditions also in-phase relationships between SPI and SRI are evident mainly in low-frequency IMFs of inter-annual periodicity. Long-term drought (SPI12) for lower-order IMFs lasts for longer durations up to a decade or even up to 2 decades in some cases (say IMF5 of HS2).

Fig. 5
figure 5

IMF comparison plots of short term drought (3-month)

Full size image
Fig. 6
figure 6

IMF comparison plots of moderate drought (6-month)

Full size image
Fig. 7
figure 7

IMF comparison plots of long term drought (12-month)

Full size image

SPI and SRI can be estimated for different aggregation scales (k) like 3, 6 and 12 months. The CDF of aggregated series was fitted using Gamma distribution following the recommendations from past studies (Thomas et al. 2015). However, an experimentation using different probability distributions can also be made to choose the distribution. (Lee et al. 2023). In such cases, a normality test is essential to be performed. Finally, an equipropability transformation of fitted CDF to CDF of standard normal distribution with zero mean and unit variance is performed, to get the SPI. In our case, as the most commonly used Gamma distribution is fixed a priori, standard normal probability distribution is used during the transformation of precipitation to SPI and runoff to SRI, ensuring the normality. However, before proceeding with TDIC analysis, this verification is carried out to confirm the normality of all SPI and SRI series by performing the best fit and its confirmation using Akaike Information Criteria (AIC) (Akaike 1973) and considering the minimum value as the best fit. The relevant confirmation of normality of the SRI and SPI series of all the stations across different scales is well proven by Figure S3, Figure S4 and Table S9 of SI file. After confirming the normality, TDI correlation analysis is performed by finding the running correlation between the corresponding IMFs of SRI and SPI at different process scales.

The TDI correlation analysis shows some interesting observations for each case. For the short-term drought series (Fig. 8), the association between SPI and SRI is mostly negative, except for the annual periodic scale, where there is a dominance of a negative correlation between the two. However, the strength of their association differs over time scales and the time domain, which is typical of multiscale processes. In IMF6 of the station 1, association and IMF5 of station 3, etc., the SPI–SRI association is strong, long-range, positive, and consistent. One can see many alterations in the streamflow-precipitation drought associations, implying that many other local meteorological processes might have controlled the associations. The changes in precipitation patterns, including climatic changes or anthropogenic interferences on streamflow may also influence the precipitation (meteorological)—hydrological drought interplays. One can recollect that SRI with the most correlated case of SPI is considered in this study from a multi-site perspective. It need not be the SPI of the same station/location that influence the SRI of a given hydrometric station.

Fig. 8
figure 8

TDIC plots between 3-month SRI and SPI at different process scales. The plots in rows 1 to row 5 show the results of hydrometric stations 1 to 5 respectively. The void space in TDIC plots shows that the correlation is not significant at that time scale (window size) and time spell. Y-axis label is Size of Sliding Window (SSW), which indicates the time scale, while IMF indicates the process scale. X-axis label is showing the Centre Position of Sliding Window (CPSW)

Full size image

Analysis of the 6-month time scale also shows the overall dominance of a positive correlation between SRI and SPI for different stations. Here, the negative dominancy in an association is noted only at intra-seasonal scales of 1.5 months for stations 1 and 2. When comparing the TDIC plots of different IMFs of a given station, the bottom contours show a progressively upwards shift, which implies that significant correlations are noted only on considering larger windows. Strong and consistent correlations are pointed out in IMF1 and IMF2 of all the stations, even though a transition-like correlation emerges for many localized time spells.

Also, in the case of IMF4 of station 3, as in short-term drought, a very strong association between the two indices emerges for a longer time scale (Fig. 9). A similar association was also noted for station 5, except for a localized spell in 2004. Unlike the previous cases, negative correlations are not visible in the long term for the intra-seasonal scales in the case of long-term droughts. A weak negative correlation is evident in IMF6 of station 3, but visible only on considering a longer time series span, i.e., a longer window for the analysis. Similar plots of long-term droughts are shown in Fig. 10.

Fig. 9
figure 9

TDIC plots between 6-month SRI and SPI at different process scales. The plots in row 1 to row 5 shows the results of hydrometric station 1 to 5 respectively. The void space in TDIC plots shows that the correlation is not significant at that time scale (window size) and time spell. Y-axis label is Size of Sliding Window (SSW), which indicates the time scale, while IMF indicates the process scale. X-axis label is showing the Centre Position of Sliding Window (CPSW)

Full size image
Fig. 10
figure 10

TDIC plots between 12-month SRI and SPI at different process scales. The plots in rows 1 to 5 show the results of hydrometric stations 1 to 5 respectively. The void space in TDIC plots shows that the correlation is not significant at that time scale (window size) and time spell. Y-axis label is Size of Sliding Window (SSW), which indicates the time scale, while IMF indicates the process scale. X-axis label is showing the Centre Position of Sliding Window (CPSW)

Full size image

In all three cases, the switchover or transition from positive to a very strong negative correlation was noted between 1984 and 1994, especially in stations 4 and 5. It is to be recollected that the period of significant drought events recognized by other techniques also shows that the maximum droughts occurred between 1981–1984 and 1988–1990. In general, this precipitation variability is influenced over the complete basin and is reflected in different TDIC plots. This instability is primarily due to the influence of large-scale climatic oscillations like El Nino events on the drought characteristics of Northern Algeria (Zeroual et al. 2017; Taibi et al. 2022).

In general, considering IMF1 and IMF2 with intra annual periodicity, the SPI–SRI correlations are positive regardless of the drought type. Here the alterations in the nature of correlations from positive to negative (or vice versa) are very localized in time domain. At annual scales, the transitions last longer. For example at stations 1, 4 and 5 during 1990s, the SPI6-SRI6 correlation is negative, whereas for HS2, the relation is negative 2002–2005 period. For the moderate and long term drought types, the SPI–SRI correlation is very strong, long range and positive for IMF5 of inter annual periodic scale. The SPI–SRI correlations are mostly negative or weak positive at the intra annual time scale of periodicity (IMF4). Overall, the alterations time from positive to negative is not same across the IMFs, stations and the drought type, which indicates the dynamic behavior of the SPI–SRI correlations.

The correlation between the trends of SRI of different stations with SPI of the most correlated station is given in Table 6. The comparison of long-term trends using simple Pearson correlations shows that in station 3, for all the types of droughts, the SPI–SRI association is negative. The associations are negative for medium and long-term droughts in stations 3–5 of the basin. For station 1 and station 2, the association between the two is positive for all the types of drought. Furthermore, the inherent trend components of the SRI series of five HS are compared with that of the corresponding SPI to enable comparison of long-term evolution droughts, as shown in Fig. 11. The results show a good similarity between SPI3 and SPI12 of station 1. Considering the case of HS3, the zero crossing happens more or less similarly instantly (year), but the behavior is the opposite. In most cases, there is a clear shift in the points of zero crossing of the trend of the two types of drought series. This indicates that the meteorological drought of the basin in the form of precipitation alone does not represent its hydrological drought. Further analysis considering droughts considering, soil moisture, evapotranspiration and vegetation characteristics are required to reinforce this finding and its reasoning.

Table 6 Correlation between the trends of SRI of different stations with SPI of most correlated station. Significant correlations (at 5% level) are underlined
Full size table
Fig. 11
figure 11

Comparison of the long-term evolution of SRI trends of five HS with different SPI droughts. Upper panels show the results of short-term droughts (SPI3), middle panels show the results of moderate droughts and lower panels show the results of long-term droughts

Full size image

The multiscale correlation analysis has the potential to allow a better understanding of drought propagation studies. The cross-wavelet power spectra and wavelet coherence analysis help get an overall association in drought teleconnections. Still, the associations at definite scales can be deciphered in a detailed way using TDIC analysis. Further hybrid decomposition-machine learning models for drought predictions can be developed based on multiscale correlation studies (Adarsh and Janga Reddy 2018, 2019; Johny et al. 2022a, 2022b, 2022c). Thus, the study has a lot of potential in Algeria’s water resources planning and management. Even though we have proposed a novel framework, the study is not free from limitations. We have analyzed the link between SPI and SRI in a multiscale perspective in a statistical manner, but the exact physical reasoning behind such dynamic relationships is not analyzed. Here, the hydrological response of the watershed to precipitation is an integrated effect of area, climatic, altitude and slope. Any alternations in such driving factors may influence the hydrological drought of the region. Therefore, in order to decipher the reasoning behind such dynamic behavior in SPI–SRI, the in-depth local scale analysis of watershed and climatic characteristics must also be integrated with the TDIC-based multiscale statistical analysis. The dynamical correlation analysis performed in this paper is a purely statistical work performed based on data. In order to get the physical reasoning, the information of local scale meteorological variables for long term is required, which is not available for all the station locations. Once such data is available, the meteorological modeling, which considers local scale processes and the physics of the problem (like governing equations) needs to be applied in order to decipher the transitions of drought completely. In general such huge, precise and local scale data availability and high computational requirement forces the modelers to use alternative statistical methods. Moreover, the study can be performed in a lagged framework employing the time-dependent intrinsic cross-correlation (TDICC) method to capture the in-depth association between SPI and SRI accounting for the antecedent information in each time scale (Johny et al. 2022a, 2022b), to model the relationship between drought types in depth in future studies.

The study's main limitation is that the data's normality test is not applied following a formal procedure when calculating the indices. This may reduce the calculation accuracy of the SPI and SRI indicators used. For this purpose, Bargaoui and Jemai (2022) stated that when calculating SPI indices, it should be shown that the data is close to normal by applying normality tests such as Shapiro–Wilk, and Kolmogorov–Smirnov within a 95% confidence interval. However, while calculating SPI values in the study of Katipoğlu and Acar (2021) applied Kolmogorov–Smirnov, Anderson–Darling and Chi-square (χ2) tests to fit Gamma, Lognormal, Generalized Extreme Value (GED), Weibull, Normal, Logistic and Log-logistic distributions to precipitation data, and SPI values were calculated according to the best distribution. As a result of the analysis, it was revealed that the GED distribution represents the precipitation in the Euphrates basin of Turkey better than the Gamma distribution, and there is an average deviation of 2% in the calculated SPI values. However, before proceeding with TDIC analysis, this verification is carried out to confirm the normality of all SPI and SRI series by performing the best fit and its confirmation.

5 Conclusion and recommendations

Meteorological drought can be a warning sign for hydrological drought, and modeling this relationship is vital to plan the utilization of water resources effectively. This study firstly analyzed the characteristics of the Standardized Precipitation Index (SPI) and Standardized Runoff Index (SRI) of Wadi Mina Basin, Algeria, at different time scales of 1–24 months. In order to analyze the relationship between the two types of droughts in multiple time scales, a novel framework based on Time-Dependent Intrinsic Correlation (TDIC) is proposed. The most correlated SPI and SRI are decomposed using an Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (ICEEMDAN) method. Subsequently, the TDIC method is applied to analyze the running correlation between the components of SPI and SRI with similar periodicity. The important conclusions drawn from the study are:

  • The temporal variability of SPI and SRI values of different stations displayed similar fluctuations and the start and end times of the drought characteristics are found to be coinciding well. It was also found that the maximum droughts occurred between 1981–1984, 1988–1990, and 2001 and 2004 and the drought durations increased with the increase in the time scales.

  • The comparison of modes of SPI and SRI with similar periodicity showed that there exists in-phase relationships between SPI and SRI in all types of droughts at inter-annual scales. The low-frequency modes and trend components of SPI and SRI are well correlated, despite the differences in the pattern of their evolution.

  • TDIC analysis has proven that the SPI–SRI relationships of different time scales bear a dynamic behavior associated with frequent alterations in nature and strength, over the different time scales and along the time domain. However, the time at which alterations in correlation behavior switches from positive to negative is not same across the modess, stations and the drought type.

  • At the intra-annual scales the SPI–SRI correlations are weak positive with localized alterations to negative along the time domain, whereas the relationship is strongly positive and long-range in the inter-annual scales up to 4 years.

  • This dynamic behavior in the SPI–SRI association along with the comparison of components at diverse scales indicates that the relationship between rainfall processes are not directly transferred to streamflow drought but it also gets controlled to many other local meteorological processes.

Increasing drought severity may negatively affect agricultural production, leading to food insecurity, irrigation and drinking water problems, disruptions in hydroelectric energy production, and deterioration of natural ecosystems. For this reason, policymakers and water managers need to implement practices to manage drought risks and water conservation effectively.