## Abstract

Hydrological drought forecasting is considered a key component in water resources risk management. As sustained meteorological drought may lead to hydrological drought over time, it is conceptually feasible to capitalize on the dependency between the meteorological and hydrological droughts while trying to forecast the latter. As such, copula functions are powerful tools to study the propagation of meteorological droughts into hydrological droughts. In this research, monthly precipitation and discharge time series were used to determine Standardized Precipitation Index (SPI) and Standardized Hydrological Drought Index (SHDI) at different time scales which quantify the state of meteorological and hydrological droughts, respectively. Five Archimedean copula functions were adopted to model the dependence structure between meteorological/hydrological drought indices. The Clayton copula was identified for further investigation based on the *p*-value. Next, the conditional probability and the matrix of forecasted class transitions were calculated. Results indicated that the next month's SHDI class forecasting is promising with less than 10% error. Moreover, extreme and severe meteorological drought classes lead to hydrological drought condition with a more than 70% probability. Other classes of meteorological drought/wet conditions lead to normal hydrological (drought) condition with less than 50% probability and to wet hydrological condition with over 20% probability.

## INTRODUCTION

Drought is a recurring phenomenon that potentially leads to water shortage in comparison with the average condition in a given period of time. It temporarily affects economy, agriculture, environment, and even domestic water supply in more extreme cases. Thus, a drought early warning system based on a robust drought forecasting scheme is in demand.

There are various drought classifications. Wilhite & Glantz (1985) proposed drought types as meteorological, hydrological, agricultural, and socio-economic. Meteorological drought is defined as insufficient precipitation relative to normal condition (Mishra & Singh 2010). Also, hydrological drought is defined as a period of time in which the amount of available water, river discharge or reservoir level is less than the normal condition. Hydrological drought often leads to water scarcity when the water supply fails to meet demands (Dehghani *et al.* 2014). As streamflow has been a major source to supply water for the environment and human needs, hydrological drought forecasting is considered a key component in water resources risk management.

The majority of drought impacts are due to hydrological drought (Van Loon 2015). A number of techniques have been applied during the past decades for drought forecasting, including regression analysis (e.g., Stagge *et al.* 2015; Meng *et al.* 2016; Sohn & Tam 2016), time series modeling (e.g., Chen *et al.* 2016; Mahmud *et al.* 2016; Tian *et al.* 2016; Karthika *et al.* 2017), probability models (e.g., Khadr 2016; Rahmat *et al.* 2016; Sun *et al.* 2016; Zhang *et al.* 2017), machine learning (e.g., Borji *et al.* 2016; Deo & Sahin 2016; Kousari *et al.* 2017; Seibert *et al.* 2017), hybrid models (e.g., Deo & Sahin 2016; Memarian *et al.* 2016; Prasad *et al.* 2017), and dynamic modeling (e.g., Bowden *et al.* 2016; Dehghani *et al.* 2017). All forecast techniques have limitations and advantages, as reviewed by Mishra & Singh (2011) and Fung *et al.* (2019). Since occurrence, duration, severity, start and end of droughts are all probabilistic phenomena, the methods used to capture drought characteristics have a statistical and probabilistic nature (Dracup *et al.* 1980).

For accurate drought forecasting, it is essential to understand drought mechanism including how the drought changes from one type to the other. It is widely known that droughts originate from precipitation deficit, i.e., meteorological drought. Prolonged meteorological drought may then result in hydrological drought. Thus, meteorological and hydrological droughts are closely related (Heim 2002). The dependency of hydrological drought on meteorological drought has been studied, among others, by Wong *et al.* (2013), Van Lanen *et al.* (2013), Niu *et al.* (2015), and Van Loon & Laaha (2015). From past literature as well as the physical understanding of droughts, it is inferred that meteorological drought can be a decisive predictor in hydrological drought forecasting.

Statistical models are suitable tools to study the dependency structure between meteorological and hydrological droughts. Given this and the probabilistic nature of drought, copula functions are often recommended to simulate the dependency structure. Copulas have been given more attention in hydrology in recent years (Singh & Strupczewski 2007). There are a number of reported cases that successfully adopted copulas in drought analysis (e.g., Kao & Govindaraju 2010; Golian *et al.* 2012; Kuchment & Demidov 2013; Maity *et al.* 2013; Huang *et al.* 2014; Saghafian & Mehdikhani 2014; Borgomeo *et al.* 2015; Chen *et al.* 2015; Jeong & Lee 2015; Hao *et al.* 2016; Amirataee *et al.* 2018; Ayantobo *et al.* 2018; Montaseri *et al.* 2018; Van de Vyver & Van den Bergh 2018; Yang *et al.* 2018).

The objective of this study is to propose a novel approach for hydrological drought forecasting using copula functions while taking antecedent meteorological drought as the main predictor. In the case of drought prediction, two approaches are considered: forecasting the index first then convert it into class and forecasting the drought class. Both approaches are suitable for drought class forecasting although both have limitations. Forecasting the drought class transition and determining the matrix of the probability of drought class transition can provide an overview of drought class transition. In the present research, for the first time, the probabilistic forecasting of hydrological drought class change on the basis of meteorological drought is addressed. For this purpose, the widely used SPI quantifies the meteorological drought, while, similar to the SPI, the SHDI characterizes the hydrological drought. Archimedean copula is used to forecast SHDI in various aggregation time scales and the matrix of forecasted transition classes is extracted. The remainder of the paper is organized as follows. The section below presents the case study and data sets. Then the methodology of model development is described. The following sections represent the case study and methodology and discuss the results of forecasting and extracting the matrix of the probability of drought class transition. Finally, conclusions of the study are presented.

## CASE STUDY AND DATA

Karun basin is located between 49°30′E to 52°00′W latitude and 30°00′ to 32°30′N longitude in southwest Iran (Figure 1). The study basin's outlet lies at Karun 3 Dam reservoir. The basin area at Pole-Shaloo hydrometric station, just upstream of Karun 3 Dam, is 24,202 km^{2}. Karun is a mountainous basin with a wide range of altitudes from 70 m in Pole-Shaloo to over 4,500 m in Dena and Zardkouh mountains. The data of 14 meteorological stations and one hydrometric station were provided by the Iranian Ministry of Energy. Figure 2 shows the meteorological stations in the basin.

Available monthly precipitation and discharge time series cover the 1973–2011 period. To check data quality, nonparametric tests were utilized. For this purpose, the Mann–Kendall test was used for trend detection. As the Mann–Kendall test is sensitive to serial correlation in the time series, at first, the Durbin–Watson test was adopted to test the serial correlation in the precipitation series. Also the normal standard homogeneity test (NSHT) and Buishand's test (BT) were used for homogeneity. The detailed description of these tests is available, for example, in Winingaard *et al.* (2003) and Vezzoli *et al.* (2012). Data gaps were filled using linear regression with nearby stations.

## METHODOLOGY

### SPI and SHDI calculations

There are several indices used for hydrological drought analysis. Surface Water Supply Index (SWSI) and Palmer Hydrological Drought Index (PHDI) are two of the most well known. However, these two indices are data intensive and difficult to apply. For hydrological drought monitoring and forecasting, a simple index based on readily available discharge data is desirable. Such index should be normal in time and space so that comparison of drought severity and duration at different temporal and spatial scales are possible.

McKee *et al.* (1993) proposed SPI for meteorological drought. Monitoring SPI requires long-term (at least 30 years) monthly precipitation data where the probability distribution function (PDF) of precipitation is assumed to follow the gamma distribution. The cumulative distribution is then transformed, using equal probability, to a normal distribution with a mean of zero and standard deviation of unity. A given precipitation for a specified time period corresponds to a particular SPI value consistent with the normal probability. Positive SPI values indicate greater than median precipitation, while negative values indicate less than median precipitation.

*et al.*(2014) proposed SHDI as a hydrological drought index by replacing precipitation with discharge, singling out the assumption of gamma distribution. First, the best PDF is fitted to the discharge time series at different time scales. Contrary to the SPI, however, the best PDF does not have to be limited to gamma function. Then, the PDF is transformed to an equal probability normal distribution time series with a mean of zero and standard deviation of unity. The SHDI reports hydrological drought/wet status and is relatively easy to compute. The following steps outline the procedure for SPI calculation. The SHDI procedure is similar. Suppose that

*x*represents precipitation time series of a given time scale. Assuming that gamma probability density function is the best fit for monthly precipitation series, we have: where is the gamma function and and are the shape and scale parameters, respectively. Shape and scale parameters are estimated using maximum likelihood estimation (MLE). Also the Akaike information criterion (AIC) and the Kolomogrov–Smirnov and Anderson–Darling test statistics are used to test the goodness of fit. Finally, an equal probability transformation is used to transform the gamma cumulative probability into the standard normal function, thus resulting in the SPI values.

Table 1 shows the classification of SPI and SHDI. For more information on the SPI and SHDI, one may refer to Edwards & McKee (1997), Hayes *et al.* (1999), Lana *et al.* (2001), Hughes & Saunders (2002), Wu *et al.* (2005, 2007), and Dehghani *et al.* (2014).

Class | Value |
---|---|

Extremely wet | 2 and above |

Severely wet | 1.5 to 1.99 |

Moderately wet | 1 to 1.49 |

Slightly wet | 0.5 to 0.99 |

Near normal | −0.49 to 0.49 |

Slightly dry | −0.5 to −1 |

Moderately dry | −1 to −1.49 |

Severely dry | −1.5 to −1.99 |

Extremely dry | −2 and less |

Class | Value |
---|---|

Extremely wet | 2 and above |

Severely wet | 1.5 to 1.99 |

Moderately wet | 1 to 1.49 |

Slightly wet | 0.5 to 0.99 |

Near normal | −0.49 to 0.49 |

Slightly dry | −0.5 to −1 |

Moderately dry | −1 to −1.49 |

Severely dry | −1.5 to −1.99 |

Extremely dry | −2 and less |

### Copula function

A two-dimensional copula is a function from to which satisfies the following two properties:

For modeling the bivariate and multivariate distributions, it is necessary to consider both marginals' effect and dependency between the marginals. The copula is able to analyze the dependency structure of bivariate and multivariate distributions without considering the marginals' effects. As the copula is defined on , every marginal is desired which can be useful in the modeling and simulation.

*C*such that:

*C*is unique. From Equation (1), the copula

*C*is given by: where and are quasi-inverses of and , respectively.

*et al.*2013; Zhang

*et al.*2015; Chang

*et al.*2016). An Archimedean copula is defined as follows: where the function is called a generator of the copula and is a mapping which satisfies the following conditions (Genest

*et al.*2006):

It is necessary to choose the best copulas which are fitted to the data. For this purpose, a goodness of fit (GOF) test is required that examines whether the unique underlying copula *C* belongs to the parameter class of Archimedean copulas (Genest *et al.* 2006). Thus, the null hypothesis is expressed by .

Several GOF tests are available; we adopted the one proposed by Genest *et al.* (2006). The procedure is as follows:

In this study, five Archimedean copulas were used for meteorological/hydrological drought modeling as presented in Table 2.

- Genest
*et al.*(2006) proposed two statistics for testing the GOF as follows: where known as Kendall's process. Equation (12) can be expanded by a simple mathematical calculation into: and

Based on these statistics, the null hypothesis may be rejected when the observed value of or is greater than percentile of its distribution (critical value). The critical value and *p*-value for any test statistic may be calculated using the parametric bootstrap method, whose procedure is outlined as follows.

Since depends on the unknown association parameter (), using , the copula parameter is estimated by .

Generate

*N*random samples of size*n*from copula and determine the value of the test statistics and for each generated sample.Suppose and are the ordered values of the test statistics and calculated from

*N*samples generated previously. Then, the critical values of the test at level based on and are and , where denotes the integer part of*a*. Also, estimates of the*p*-value based on the values of and are given by and , where and are the observed values of and .

### Conditional probability

*X*and

*Y*, respectively. Then, the conditional cumulative distribution function of

*X*given is given by: where and are the probability density functions of and and is the density of the copula . Furthermore, the conditional probability of SHDI class given the SPI in previous time steps may be determined as follows: where, and

Now, it is required to assign a marginal distribution to the SPI and SHDI. In this study, six commonly used univariate probability distributions were selected as the candidate marginal distribution (Table 3). They are normal, Pearson type III, log-normal (3P), generalized extreme value, logistic, and generalized logistic. The parameters of these marginal distributions were estimated via MLE method and the goodness of fit was conducted based on Kolomogrov–Smirnov and Anderson–Darling test statistics. Akaike information criterion (AIC), log-likelihood (LL), and *p*-value were calculated to choose the best fitted marginal distribution.

Copula | |||||
---|---|---|---|---|---|

Clayton | |||||

Frank | |||||

Gumbel–Hougaard | |||||

Joe | |||||

Ali–Mikhail–Hag |

Copula | |||||
---|---|---|---|---|---|

Clayton | |||||

Frank | |||||

Gumbel–Hougaard | |||||

Joe | |||||

Ali–Mikhail–Hag |

Distribution | CDF | Parameters |
---|---|---|

Normal | _{;} is the CDF of standard normal distribution | , location; , scale |

Pearson type III | , shape; , scale; , location | |

Log-normal (3P) | _{;} is the CDF of standard normal distribution | , location; , scale; , shape |

Generalized extreme value | , location; , scale; , shape | |

Logistic | , location; , scale | |

Generalized logistic | , location; , scale; , shape |

Distribution | CDF | Parameters |
---|---|---|

Normal | _{;} is the CDF of standard normal distribution | , location; , scale |

Pearson type III | , shape; , scale; , location | |

Log-normal (3P) | _{;} is the CDF of standard normal distribution | , location; , scale; , shape |

Generalized extreme value | , location; , scale; , shape | |

Logistic | , location; , scale | |

Generalized logistic | , location; , scale; , shape |

### Forecast model development

In order to address the challenge of determining the best meteorological/hydrological drought forecast combination, the correlation coefficients and cross correlation between the SPI and SHDI at different time scales were calculated. Combinations with the highest correlation value were selected for further analysis. In the next step, for each copula function in Table 2, the dependence parameter under the null hypothesis was estimated by inversion of Kendall's tau; that is, by setting . Then, using parametric bootstrap, 2,000 samples of size *n* for each copula were generated and the *p*-value and critical value were computed in order to select the best copula and determine the marginal distributions.

The SHDI forecasting framework is shown in Figure 3.

## RESULTS

In order to extract reliable SPI and SHDI time series, good quality precipitation and discharge time series are needed. Rigorous quality tests were carried out to guarantee the quality of raw data. Results of SHNT, BT, Durbin–Watson, and Mann–Kendal tests are presented in Table 4, according to which, one could claim that the precipitation and discharge time series are homogeneous, serially independent with no trend.

Stations | Homogeneity tests | Serial correlation test | Trend test | ||
---|---|---|---|---|---|

SHNT | BT | Durbin–Watson | Mann–Kendal | ||

Precipitation | Armand | 0.77 | 0.86 | 0.87 | 0.77 |

Barez | 0.7 | 0.63 | 0.76 | 0.68 | |

Borojen | 0.39 | 0.13 | 0.12 | 0.1 | |

Chehelgerd | 0.95 | 0.98 | 0.46 | 0.68 | |

Hana | 0.8 | 0.91 | 0.81 | 0.66 | |

Kata | 0.86 | 0.93 | 0.8 | 0.81 | |

Lordegan | 0.74 | 0.98 | 0.26 | 0.75 | |

Pataveh | 0.76 | 0.47 | 0.29 | 0.32 | |

Pole-Shaloo | 0.85 | 0.97 | 0.66 | 0.83 | |

Sepidar | 0.77 | 0.48 | 0.31 | 0.8 | |

Shahid | 0.66 | 0.71 | 0.32 | 0.71 | |

Shahrekord | 0.89 | 0.61 | 0.83 | 0.54 | |

Solakan | 0.49 | 0.16 | 0.38 | 0.4 | |

Yasooj | 0.77 | 0.9 | 0.34 | 0.45 | |

Discharge | Pole-Shaloo | 0.2 | 0.18 | 0.6 | 0.06 |

Stations | Homogeneity tests | Serial correlation test | Trend test | ||
---|---|---|---|---|---|

SHNT | BT | Durbin–Watson | Mann–Kendal | ||

Precipitation | Armand | 0.77 | 0.86 | 0.87 | 0.77 |

Barez | 0.7 | 0.63 | 0.76 | 0.68 | |

Borojen | 0.39 | 0.13 | 0.12 | 0.1 | |

Chehelgerd | 0.95 | 0.98 | 0.46 | 0.68 | |

Hana | 0.8 | 0.91 | 0.81 | 0.66 | |

Kata | 0.86 | 0.93 | 0.8 | 0.81 | |

Lordegan | 0.74 | 0.98 | 0.26 | 0.75 | |

Pataveh | 0.76 | 0.47 | 0.29 | 0.32 | |

Pole-Shaloo | 0.85 | 0.97 | 0.66 | 0.83 | |

Sepidar | 0.77 | 0.48 | 0.31 | 0.8 | |

Shahid | 0.66 | 0.71 | 0.32 | 0.71 | |

Shahrekord | 0.89 | 0.61 | 0.83 | 0.54 | |

Solakan | 0.49 | 0.16 | 0.38 | 0.4 | |

Yasooj | 0.77 | 0.9 | 0.34 | 0.45 | |

Discharge | Pole-Shaloo | 0.2 | 0.18 | 0.6 | 0.06 |

After performing the nonparametric tests and filling data gaps, different probability distributions were fitted to the monthly precipitation and discharge time series. Gamma and log-normal functions were determined to be the best fit for SPI and SHDI, respectively (samples shown in Tables 5 and 6). Next, the SPI and SHDI time series were determined for time scales of 3, 6, 9, 12 and 1, 2, 3 months, respectively. Hereafter, for example, SHDI1 stands for SHDI corresponding to one month time scale. Basin average of SPI was determined using kriging method. It is worth noting that estimated values from equations of SPI and SHDI were used directly and then the predicted values categorized based on the drought intensity for conditional drought forecasting.

Month | Log-likelihood | AIC | K-S test | A-D test | ||||
---|---|---|---|---|---|---|---|---|

Statistic | p-value | Statistic | p-value | |||||

January | 6.273 | 49.291 | −230.456 | 437.265 | 0.12 | 0.824 | 0.225 | 0.962 |

February | 8.173 | 51.892 | −240.326 | 460.332 | 0.074 | 0.932 | 0.156 | 0.944 |

March | 7.043 | 178.436 | −227.711 | 463.423 | 0.09 | 0.918 | 0.222 | 0.983 |

April | 12.983 | 47.406 | −228.736 | 465.473 | 0.057 | 1 | 0.143 | 0.999 |

May | 12.836 | 20.574 | −235.767 | 479.535 | 0.062 | 0.998 | 0.123 | 1 |

June | 9.164 | 59.229 | −231.174 | 470.348 | 0.084 | 0.943 | 0.253 | 0.969 |

July | 6.325 | 115.036 | −223.543 | 455.086 | 0.088 | 0.917 | 0.213 | 0.987 |

August | 5.015 | 242.911 | −213.741 | 435.482 | 0.1 | 0.827 | 0.219 | 0.984 |

September | 1.236 | 86.662 | −197.78 | 403.559 | 0.129 | 0.542 | 0.629 | 0.62 |

October | 2.736 | 18.034 | −165.25 | 376.223 | 0.14 | 0.432 | 0.724 | 0.425 |

November | 1.546 | 44.723 | −180.42 | 400.432 | 0.125 | 0.526 | 0.313 | 0.652 |

December | 2.309 | 77.194 | −209.373 | 426.746 | 0.084 | 0.948 | 0.274 | 0.956 |

Month | Log-likelihood | AIC | K-S test | A-D test | ||||
---|---|---|---|---|---|---|---|---|

Statistic | p-value | Statistic | p-value | |||||

January | 6.273 | 49.291 | −230.456 | 437.265 | 0.12 | 0.824 | 0.225 | 0.962 |

February | 8.173 | 51.892 | −240.326 | 460.332 | 0.074 | 0.932 | 0.156 | 0.944 |

March | 7.043 | 178.436 | −227.711 | 463.423 | 0.09 | 0.918 | 0.222 | 0.983 |

April | 12.983 | 47.406 | −228.736 | 465.473 | 0.057 | 1 | 0.143 | 0.999 |

May | 12.836 | 20.574 | −235.767 | 479.535 | 0.062 | 0.998 | 0.123 | 1 |

June | 9.164 | 59.229 | −231.174 | 470.348 | 0.084 | 0.943 | 0.253 | 0.969 |

July | 6.325 | 115.036 | −223.543 | 455.086 | 0.088 | 0.917 | 0.213 | 0.987 |

August | 5.015 | 242.911 | −213.741 | 435.482 | 0.1 | 0.827 | 0.219 | 0.984 |

September | 1.236 | 86.662 | −197.78 | 403.559 | 0.129 | 0.542 | 0.629 | 0.62 |

October | 2.736 | 18.034 | −165.25 | 376.223 | 0.14 | 0.432 | 0.724 | 0.425 |

November | 1.546 | 44.723 | −180.42 | 400.432 | 0.125 | 0.526 | 0.313 | 0.652 |

December | 2.309 | 77.194 | −209.373 | 426.746 | 0.084 | 0.948 | 0.274 | 0.956 |

Month | Log-likelihood | AIC | K-S test | A-D test | ||||
---|---|---|---|---|---|---|---|---|

Statistic | p-value | Statistic | p-value | |||||

January | 0.423 | 6.497 | −254.024 | 512.049 | 0.142 | 0.423 | 0.664 | 0.589 |

February | 0.475 | 6.283 | −250.425 | 504.85 | 0.116 | 0.676 | 0.469 | 0.777 |

March | 0.45 | 5.806 | −231.334 | 466.668 | 0.163 | 0.264 | 0.812 | 0.471 |

April | 0.423 | 5.361 | −213.137 | 430.275 | 0.128 | 0.555 | 0.899 | 0.414 |

May | 0.393 | 5.038 | −198.898 | 401.796 | 0.119 | 0.641 | 0.544 | 0.701 |

June | 0.346 | 4.813 | −186.159 | 376.317 | 0.118 | 0.653 | 0.596 | 0.651 |

July | 0.275 | 4.779 | −176.735 | 357.471 | 0.106 | 0.771 | 0.349 | 0.897 |

August | 0.378 | 4.999 | −195.99 | 395.98 | 0.158 | 0.299 | 0.928 | 0.397 |

September | 0.554 | 5.358 | −222.678 | 449.356 | 0.194 | 0.115 | 1.707 | 0.134 |

October | 0.515 | 5.52 | −225.951 | 455.903 | 0.101 | 0.825 | 0.329 | 0.915 |

November | 0.525 | 5.749 | −234.854 | 473.708 | 0.077 | 0.971 | 0.268 | 0.96 |

December | 0.458 | 6.202 | −246.203 | 496.406 | 0.084 | 0.943 | 0.268 | 0.96 |

Month | Log-likelihood | AIC | K-S test | A-D test | ||||
---|---|---|---|---|---|---|---|---|

Statistic | p-value | Statistic | p-value | |||||

January | 0.423 | 6.497 | −254.024 | 512.049 | 0.142 | 0.423 | 0.664 | 0.589 |

February | 0.475 | 6.283 | −250.425 | 504.85 | 0.116 | 0.676 | 0.469 | 0.777 |

March | 0.45 | 5.806 | −231.334 | 466.668 | 0.163 | 0.264 | 0.812 | 0.471 |

April | 0.423 | 5.361 | −213.137 | 430.275 | 0.128 | 0.555 | 0.899 | 0.414 |

May | 0.393 | 5.038 | −198.898 | 401.796 | 0.119 | 0.641 | 0.544 | 0.701 |

June | 0.346 | 4.813 | −186.159 | 376.317 | 0.118 | 0.653 | 0.596 | 0.651 |

July | 0.275 | 4.779 | −176.735 | 357.471 | 0.106 | 0.771 | 0.349 | 0.897 |

August | 0.378 | 4.999 | −195.99 | 395.98 | 0.158 | 0.299 | 0.928 | 0.397 |

September | 0.554 | 5.358 | −222.678 | 449.356 | 0.194 | 0.115 | 1.707 | 0.134 |

October | 0.515 | 5.52 | −225.951 | 455.903 | 0.101 | 0.825 | 0.329 | 0.915 |

November | 0.525 | 5.749 | −234.854 | 473.708 | 0.077 | 0.971 | 0.268 | 0.96 |

December | 0.458 | 6.202 | −246.203 | 496.406 | 0.084 | 0.943 | 0.268 | 0.96 |

The forecasting technique relies on the correlation between the SPI and SHDI in different time scales. High correlation indicates that the antecedent meteorological drought has a significant impact on hydrological drought. For this purpose, the correlation coefficient and cross correlation between the SPI and SHDI in different time scales were calculated and presented in Table 7. One may observe that the correlation coefficients between SPI6, SPI9, and SPI12 with SHDI1(t + 1), SHDI2(t + 2), and SHDI3(t + 3) are very close, while dropping to lower values for SPI3. The results indicate that a lag of more than three months between meteorological and hydrological droughts exists. Among the remaining SPI-SHDI combinations, the correlation coefficient varies negligibly. Thus, all three SPI time scales (i.e., 6, 9, and 12) were considered in conjunction with one SHDI time scale. As a result, SPI6-SHDI1(t + 1), SPI9-SHDI2(t + 2), and SPI12-SHDI3(t + 3) predictor–predictant combinations were selected as the candidates for further investigation. This implies that, for instance, in the case of SPI6-SHDI1(t + 1), SPI6 (embodying precipitation of six months up to a given month t) was used to forecast the SHDI1 of the next month (t + 1). The cross correlation results in Figure 4 also verify the results in Table 7.

SHDI1(t + 1) | SHDI2(t + 2) | SHDI3(t + 3) | |
---|---|---|---|

SPI3 | 0.46 | 0.44 | 0.42 |

SPI6 | 0.58 | 0.55 | 0.51 |

SPI9 | 0.59 | 0.54 | 0.47 |

SPI12 | 0.53 | 0.49 | 0.50 |

SHDI1(t + 1) | SHDI2(t + 2) | SHDI3(t + 3) | |
---|---|---|---|

SPI3 | 0.46 | 0.44 | 0.42 |

SPI6 | 0.58 | 0.55 | 0.51 |

SPI9 | 0.59 | 0.54 | 0.47 |

SPI12 | 0.53 | 0.49 | 0.50 |

For all copula functions in Table 2, the observed value of Kendall's tau and the associated dependence parameter were estimated as presented in Table 8. For and combinations, the Kendall's tau values were 0.385 and 0.373, respectively. Thus, there is no for Ali–Mikhail–Hag copula function which satisfies . As a result, only four copulas remained for further analysis. The critical values and *p*-values reported in Table 8 were computed using 2,000 iterations within the parametric bootstrap procedure. Results indicated that at 0.05 level, only the Clayton copula was acceptable for the and combinations based on and statistics while other copulas were rejected. For the third combination, i.e., , Clayton and Ali–Mikhail–Haq copulas were acceptable while other copulas were rejected. However, we selected the Clayton over the Ali–Mikhail–Haq copula due to higher *p*-value.

Copula | Scenario | N | Critical value of | p-value of | Critical value of | p-value of | ||||
---|---|---|---|---|---|---|---|---|---|---|

Clayton | SPI6-SHDI1(t + 1) | 0.385 | 1.251 | 426 | 0.151 | 0.215 | 0.150 | 0.901 | 1.103 | 0.191 |

Frank | 0.385 | 3.956 | 426 | 0.372 | 0.173 | 0.001 | 1.374 | 1.001 | 0.002 | |

Gumbel–Hougaard | 0.385 | 1.626 | 426 | 0.515 | 0.194 | 0.000 | 1.624 | 1.028 | 0.000 | |

Joe | 0.385 | 2.141 | 426 | 1.681 | 0.722 | 0.000 | 2.203 | 1.032 | 0.000 | |

Ali–Mikhail–Hag | 0.385 | – | 426 | – | – | – | – | – | – | |

Clayton | SPI9-SHDI2(t + 2) | 0.373 | 1.191 | 422 | 0.104 | 0.214 | 0.314 | 0.618 | 1.074 | 0.717 |

Frank | 0.373 | 3.804 | 422 | 0.429 | 0.175 | 0.000 | 1.595 | 0.997 | 0.000 | |

Gumbel–Hougaard | 0.373 | 1.596 | 422 | 0.626 | 0.199 | 0.000 | 1.829 | 1.042 | 0.000 | |

Joe | 0.373 | 2.085 | 422 | 1.837 | 0.718 | 0.000 | 2.397 | 1.039 | 0.000 | |

Ali–Mikhail–Hag | 0.373 | – | 422 | – | – | – | – | – | – | |

Clayton | SPI12-SHDI3(t + 3) | 0.321 | 0.947 | 418 | 0.054 | 0.236 | 0.747 | 0.853 | 1.111 | 0.280 |

Frank | 0.321 | 3.163 | 418 | 0.571 | 0.198 | 0.000 | 1.990 | 1.038 | 0.000 | |

Gumbel–Hougaard | 0.321 | 1.473 | 418 | 0.812 | 0.207 | 0.000 | 2.099 | 1.048 | 0.000 | |

Joe | 0.321 | 1.856 | 418 | 1.969 | 0.668 | 0.000 | 2.594 | 1.042 | 0.000 | |

Ali–Mikhail–Hag | 0.321 | 0.981 | 418 | 0.064 | 0.238 | 0.643 | 0.940 | 1.103 | 0.155 |

Copula | Scenario | N | Critical value of | p-value of | Critical value of | p-value of | ||||
---|---|---|---|---|---|---|---|---|---|---|

Clayton | SPI6-SHDI1(t + 1) | 0.385 | 1.251 | 426 | 0.151 | 0.215 | 0.150 | 0.901 | 1.103 | 0.191 |

Frank | 0.385 | 3.956 | 426 | 0.372 | 0.173 | 0.001 | 1.374 | 1.001 | 0.002 | |

Gumbel–Hougaard | 0.385 | 1.626 | 426 | 0.515 | 0.194 | 0.000 | 1.624 | 1.028 | 0.000 | |

Joe | 0.385 | 2.141 | 426 | 1.681 | 0.722 | 0.000 | 2.203 | 1.032 | 0.000 | |

Ali–Mikhail–Hag | 0.385 | – | 426 | – | – | – | – | – | – | |

Clayton | SPI9-SHDI2(t + 2) | 0.373 | 1.191 | 422 | 0.104 | 0.214 | 0.314 | 0.618 | 1.074 | 0.717 |

Frank | 0.373 | 3.804 | 422 | 0.429 | 0.175 | 0.000 | 1.595 | 0.997 | 0.000 | |

Gumbel–Hougaard | 0.373 | 1.596 | 422 | 0.626 | 0.199 | 0.000 | 1.829 | 1.042 | 0.000 | |

Joe | 0.373 | 2.085 | 422 | 1.837 | 0.718 | 0.000 | 2.397 | 1.039 | 0.000 | |

Ali–Mikhail–Hag | 0.373 | – | 422 | – | – | – | – | – | – | |

Clayton | SPI12-SHDI3(t + 3) | 0.321 | 0.947 | 418 | 0.054 | 0.236 | 0.747 | 0.853 | 1.111 | 0.280 |

Frank | 0.321 | 3.163 | 418 | 0.571 | 0.198 | 0.000 | 1.990 | 1.038 | 0.000 | |

Gumbel–Hougaard | 0.321 | 1.473 | 418 | 0.812 | 0.207 | 0.000 | 2.099 | 1.048 | 0.000 | |

Joe | 0.321 | 1.856 | 418 | 1.969 | 0.668 | 0.000 | 2.594 | 1.042 | 0.000 | |

Ali–Mikhail–Hag | 0.321 | 0.981 | 418 | 0.064 | 0.238 | 0.643 | 0.940 | 1.103 | 0.155 |

The contour plot of bivariate Clayton copula and scatter plots of the (SPI, SHDI) joint variables for different combinations are shown in Figure 5. One can observe that the dependence structure of (SPI, SHDI) is similar to that of the Clayton copula. Also, the 3-D perspective plot of theoretical Clayton copula and the scatter plot of the joint SPI-SHDI combinations are drawn in Figure 6, indicating an excellent fit. Therefore, based on Table 8 statistics, it is fair to state that the Clayton copula is appropriate to model SPI-SHDI.

For conditional probability calculation, it is required to determine the marginal distribution of SPI and SHDI in different time scales. For this purpose, six univariate probability distributions, as shown in Table 3, were fitted to the SPI and SHDI in different time scales. Results are presented in Tables 9 and 10. Accordingly, for SPI of different time scales, normal, log-normal (3P), generalized extreme value, logistic, and generalized logistic were fitted successfully. Similarly, log-normal (3P), generalized extreme value, logistic, generalized logistic, and Pearson type III were fitted for SHDI of different time scales. The parameters of marginal distributions estimated via MLE method and the test statistics, AIC and LL are presented in Tables 9 and 10. Based on the AIC and LL, the generalized extreme value and generalized logistic were the best fitted marginal distributions for SPI and SHDI in all time scales, respectively.

Margin | Parameters | Log-likelihood | AIC | K-S test | A-D test | |||
---|---|---|---|---|---|---|---|---|

Statistics | p-value | Statistics | p-value | |||||

SPI6 | Normal | (,) = (0.068,0.999) | −603.89 | 1,211.79 | 0.035 | 0.663 | 0.57 | 0.677 |

Log-normal (3P) | (,,) = (3.461,0.031, −31.799) | −605.56 | 1,217.12 | 0.041 | 0.438 | 0.798 | 0.445 | |

Generalized extreme value | (,,) = (−0.252,1.035, − 0.35) | −600.91 | 1,207.8 | 0.035 | 0.662 | 0.612 | 0.637 | |

Logistic | (,) = (0.088,0.57) | −607.93 | 1,219.86 | 0.034 | 0.675 | 0.625 | 0.623 | |

Generalized logistic | (,,) = (0.106,0.569, − 0.046) | −606.89 | 1,219.79 | 0.028 | 0.873 | 0.479 | 0.767 | |

SPI9 | Normal | (,) = (0.075,0.995) | −596.55 | 1,197.11 | 0.048 | 0.282 | 0.538 | 0.71 |

Log-normal (3P) | (,,) = (3.461,0.031, − 31.8) | −597.99 | 1,201.99 | 0.054 | 0.109 | 0.715 | 0.441 | |

Generalized extreme value | (,,) = (−0.242,1.029, − 0.353) | −592.88 | 1,191.75 | 0.034 | 0.708 | 0.541 | 0.705 | |

Logistic | (,) = (0.091,0.572) | −602.05 | 1,208.11 | 0.04 | 0.478 | 0.796 | 0.487 | |

Generalized logistic | (,,) = (0.108,0.57, − 0.042) | –601.27 | 1,208.55 | 0.035 | 0.686 | 0.095 | 0.547 | |

SPI12 | Normal | (,) = (0.075,0.997) | −591.72 | 1,187.45 | 0.066 | 0.046 | 0.987 | 0.364 |

Log-normal (3P) | (,,) = (3.46,0.031, − 31.79) | −594.09 | 1,194.19 | 0.073 | 0.004 | 1.277 | 0.103 | |

Generalized extreme value | (,,) = (−0.228,1.046, − 0.384) | −585.97 | 1,177.95 | 0.044 | 0.365 | 0.764 | 0.509 | |

Logistic | (,) = (0.102,0.569) | −595.39 | 1,194.78 | 0.057 | 0.123 | 1.067 | 0.324 | |

Generalized logistic | (,,) = (0.13,0.567, − 0.070) | −593.14 | 1,192.29 | 0.045 | 0.353 | 0.868 | 0.435 |

Margin | Parameters | Log-likelihood | AIC | K-S test | A-D test | |||
---|---|---|---|---|---|---|---|---|

Statistics | p-value | Statistics | p-value | |||||

SPI6 | Normal | (,) = (0.068,0.999) | −603.89 | 1,211.79 | 0.035 | 0.663 | 0.57 | 0.677 |

Log-normal (3P) | (,,) = (3.461,0.031, −31.799) | −605.56 | 1,217.12 | 0.041 | 0.438 | 0.798 | 0.445 | |

Generalized extreme value | (,,) = (−0.252,1.035, − 0.35) | −600.91 | 1,207.8 | 0.035 | 0.662 | 0.612 | 0.637 | |

Logistic | (,) = (0.088,0.57) | −607.93 | 1,219.86 | 0.034 | 0.675 | 0.625 | 0.623 | |

Generalized logistic | (,,) = (0.106,0.569, − 0.046) | −606.89 | 1,219.79 | 0.028 | 0.873 | 0.479 | 0.767 | |

SPI9 | Normal | (,) = (0.075,0.995) | −596.55 | 1,197.11 | 0.048 | 0.282 | 0.538 | 0.71 |

Log-normal (3P) | (,,) = (3.461,0.031, − 31.8) | −597.99 | 1,201.99 | 0.054 | 0.109 | 0.715 | 0.441 | |

Generalized extreme value | (,,) = (−0.242,1.029, − 0.353) | −592.88 | 1,191.75 | 0.034 | 0.708 | 0.541 | 0.705 | |

Logistic | (,) = (0.091,0.572) | −602.05 | 1,208.11 | 0.04 | 0.478 | 0.796 | 0.487 | |

Generalized logistic | (,,) = (0.108,0.57, − 0.042) | –601.27 | 1,208.55 | 0.035 | 0.686 | 0.095 | 0.547 | |

SPI12 | Normal | (,) = (0.075,0.997) | −591.72 | 1,187.45 | 0.066 | 0.046 | 0.987 | 0.364 |

Log-normal (3P) | (,,) = (3.46,0.031, − 31.79) | −594.09 | 1,194.19 | 0.073 | 0.004 | 1.277 | 0.103 | |

Generalized extreme value | (,,) = (−0.228,1.046, − 0.384) | −585.97 | 1,177.95 | 0.044 | 0.365 | 0.764 | 0.509 | |

Logistic | (,) = (0.102,0.569) | −595.39 | 1,194.78 | 0.057 | 0.123 | 1.067 | 0.324 | |

Generalized logistic | (,,) = (0.13,0.567, − 0.070) | −593.14 | 1,192.29 | 0.045 | 0.353 | 0.868 | 0.435 |

Margin | Parameters | Log-likelihood | AIC | K-S test | A-D test | |||
---|---|---|---|---|---|---|---|---|

Statistics | p-value | Statistics | p-value | |||||

SHDI1(t + 1) | Pearson type III | (,,) = (56.46,0.13, − 7.48) | −600.6 | 1,207.2 | 0.067 | 0.039 | 2.363 | 0.046 |

Log-normal (3P) | (,,) = (2.361,0.093, − 10.647) | −600.4 | 1,206.9 | 0.067 | 0.04 | 2.311 | 0.065 | |

Generalized extreme value | (,,) = (−0.377,0.952, − 0.194) | −601.4 | 1,208.8 | 0.072 | 0.039 | 2.697 | 0.039 | |

Logistic | (,) = (−0.028,0.549) | −599.1 | 1,202.3 | 0.051 | 0.054 | 1.627 | 0.149 | |

Generalized logistic | (,,) = (−0.049,0.54,0.063) | −596.2 | 1,198.5 | 0.047 | 0.277 | 1.028 | 0.345 | |

SHDI2(t + 2) | Gamma (3P) | (,,) = (55.786,0.134, − 7.481) | −597.7 | 1,201.4 | 0.063 | 0.051 | 2.31 | 0.058 |

Log-normal (3P) | (,,) = (2.58,0.075, − 13.30) | −597.5 | 1,201.1 | 0.065 | 0.052 | 2.295 | 0.016 | |

Generalized extreme value | (,,) = (−0.369,0.966, − 0.213) | −598 | 1,202.1 | 0.067 | 0.108 | 2.628 | 0.043 | |

Logistic | (,) = (−0.019,0.556) | −597.5 | 1,199.1 | 0.053 | 0.049 | 1.7 | 0.135 | |

Generalized logistic | (,,) = (−0.036,0.554,0.053) | −595.7 | 1,197.4 | 0.047 | 0.3 | 1.257 | 0.247 | |

SHDI3(t + 3) | Gamma (3P) | (,,) = (55.20,0.137, − 7.483) | −597.9 | 1,201.9 | 0.052 | 0.217 | 2.022 | 0.092 |

Log-normal (3P) | (,,) = (2.9,0.055, − 18.10) | −596.7 | 1,199.4 | 0.05 | 0.225 | 1.794 | 0.12 | |

Generalized extreme value | (,,) = (−0.296,0.999, − 0.241) | −597.8 | 1,201.5 | 0.057 | 0.176 | 2.247 | 0.067 | |

Logistic | (,) = (0.067,0.562) | −595.6 | 1,196.8 | 0.04 | 0.394 | 0.962 | 0.377 | |

Generalized logistic | (,,) = (0.061,0.562,0.018) | −595.4 | 1,195.2 | 0.038 | 0.593 | 0.9 | 0.413 |

Margin | Parameters | Log-likelihood | AIC | K-S test | A-D test | |||
---|---|---|---|---|---|---|---|---|

Statistics | p-value | Statistics | p-value | |||||

SHDI1(t + 1) | Pearson type III | (,,) = (56.46,0.13, − 7.48) | −600.6 | 1,207.2 | 0.067 | 0.039 | 2.363 | 0.046 |

Log-normal (3P) | (,,) = (2.361,0.093, − 10.647) | −600.4 | 1,206.9 | 0.067 | 0.04 | 2.311 | 0.065 | |

Generalized extreme value | (,,) = (−0.377,0.952, − 0.194) | −601.4 | 1,208.8 | 0.072 | 0.039 | 2.697 | 0.039 | |

Logistic | (,) = (−0.028,0.549) | −599.1 | 1,202.3 | 0.051 | 0.054 | 1.627 | 0.149 | |

Generalized logistic | (,,) = (−0.049,0.54,0.063) | −596.2 | 1,198.5 | 0.047 | 0.277 | 1.028 | 0.345 | |

SHDI2(t + 2) | Gamma (3P) | (,,) = (55.786,0.134, − 7.481) | −597.7 | 1,201.4 | 0.063 | 0.051 | 2.31 | 0.058 |

Log-normal (3P) | (,,) = (2.58,0.075, − 13.30) | −597.5 | 1,201.1 | 0.065 | 0.052 | 2.295 | 0.016 | |

Generalized extreme value | (,,) = (−0.369,0.966, − 0.213) | −598 | 1,202.1 | 0.067 | 0.108 | 2.628 | 0.043 | |

Logistic | (,) = (−0.019,0.556) | −597.5 | 1,199.1 | 0.053 | 0.049 | 1.7 | 0.135 | |

Generalized logistic | (,,) = (−0.036,0.554,0.053) | −595.7 | 1,197.4 | 0.047 | 0.3 | 1.257 | 0.247 | |

SHDI3(t + 3) | Gamma (3P) | (,,) = (55.20,0.137, − 7.483) | −597.9 | 1,201.9 | 0.052 | 0.217 | 2.022 | 0.092 |

Log-normal (3P) | (,,) = (2.9,0.055, − 18.10) | −596.7 | 1,199.4 | 0.05 | 0.225 | 1.794 | 0.12 | |

Generalized extreme value | (,,) = (−0.296,0.999, − 0.241) | −597.8 | 1,201.5 | 0.057 | 0.176 | 2.247 | 0.067 | |

Logistic | (,) = (0.067,0.562) | −595.6 | 1,196.8 | 0.04 | 0.394 | 0.962 | 0.377 | |

Generalized logistic | (,,) = (0.061,0.562,0.018) | −595.4 | 1,195.2 | 0.038 | 0.593 | 0.9 | 0.413 |

In the next step, conditional probability was determined for all combinations based on the Clayton copula and the selected marginal distributions. Results plotted in Figures 7–9 indicate that for , the SHDI will be in the drought or normal state while the probability of SHDI wet class is negligible. This may be interpreted such that with a sustained meteorological drought, the surface water resources moves into the drought state with high probability, or to normal state with a very low probability. This means that droughts do not start or end abruptly. The state leads to a normal SHDI with high probability. For , the probability of SHDI drought, wet, and normal states are negligible, 0.3, and 0.7, respectively. This implies that a sustained wet meteorological condition will not necessarily result in a wet hydrological condition in the next month. In Figures 7–9, the curves corresponding to the meteorological drought states cover a wide range while the curves of wet states are very close or overlaid. This is due to the fact that will result in SHDI drought conditions in a wide range of probabilities while covers a narrow range of probability for . In general, it may be concluded that the case leads to hydrological drought state with a high probability, the case leads to normal hydrological state with high probability, and the case leads to normal hydrological state with high probability and wet hydrological state with lower probability in all SPI-SHDI combinations.

In the next step, for better judgment and to prepare a management tool for water policy (decision) makers, the forecasted transition probability matrixes of SPI-SHDI combinations were extracted. First a matrix showing the transition probability of SPI-SHDI drought-normal-wet conditions was extracted (Table 11). Decisions based on forecasts are quite straightforward if SPI and SHDI are to be classified into only three major states (i.e., drought, normal, and wet). Based on Table 11, for combination, meteorological drought state will lead to drought, normal, and wet hydrological drought states in the next month with about 64%, 28%, and 7% probability, respectively, that are in overall agreement with Figures 7–9. However, SPI normal state leads to SHDI normal state with over 50% probability and to drought and wet states with 22% and 25% probability, respectively. For the SPI wet state, SHDI normal and wet states with 46% and 41% of probability are expected respectively while the probability of drought occurrence is about 9%. Nearly similar probability values are assigned for two other combinations.

SPI/SHDI class | SPI6-SHDI1(t + 1) | SPI9-SHDI2(t + 2) | SPI12-SHDI3(t + 3) | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Drought | Normal | Wet | Drought | Normal | Wet | Drought | Normal | Wet | ||

Forecasted | Drought | 64.01 | 28.39 | 7.18 | 63.21 | 28.33 | 7.79 | 54.72 | 32.51 | 11.51 |

Normal | 22.49 | 51.09 | 25.37 | 22.98 | 49.65 | 26.12 | 20.89 | 47.41 | 30.03 | |

Wet | 9.37 | 46.65 | 41.88 | 9.93 | 45.53 | 42.09 | 10.24 | 43.27 | 43.63 | |

Observed | Drought | 61.20 | 32.76 | 6.03 | 54.92 | 36.88 | 8.20 | 51.67 | 33.30 | 15.00 |

Normal | 23.98 | 53.80 | 21.64 | 22.98 | 54.04 | 22.98 | 18.75 | 55.63 | 25.63 | |

Wet | 7.14 | 45.71 | 47.14 | 6.47 | 48.20 | 45.32 | 8.70 | 46.38 | 44.93 |

SPI/SHDI class | SPI6-SHDI1(t + 1) | SPI9-SHDI2(t + 2) | SPI12-SHDI3(t + 3) | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Drought | Normal | Wet | Drought | Normal | Wet | Drought | Normal | Wet | ||

Forecasted | Drought | 64.01 | 28.39 | 7.18 | 63.21 | 28.33 | 7.79 | 54.72 | 32.51 | 11.51 |

Normal | 22.49 | 51.09 | 25.37 | 22.98 | 49.65 | 26.12 | 20.89 | 47.41 | 30.03 | |

Wet | 9.37 | 46.65 | 41.88 | 9.93 | 45.53 | 42.09 | 10.24 | 43.27 | 43.63 | |

Observed | Drought | 61.20 | 32.76 | 6.03 | 54.92 | 36.88 | 8.20 | 51.67 | 33.30 | 15.00 |

Normal | 23.98 | 53.80 | 21.64 | 22.98 | 54.04 | 22.98 | 18.75 | 55.63 | 25.63 | |

Wet | 7.14 | 45.71 | 47.14 | 6.47 | 48.20 | 45.32 | 8.70 | 46.38 | 44.93 |

To examine the accuracy of the developed model, the probability of SHDI in different classes given an observed SPI class is presented in Table 11. The observed and forecasted probability of class transition among all SPI-SHDI combinations matched well. For SPI6-SHDI1(t + 1), the maximum difference between observed and forecasted probability of class transition is about 6%. For SPI9-SHDI2(t + 2) and SPI12-SHDI3(t + 3), the maximum difference is 9% and 8%, respectively. In general, it may be concluded that the copula model is capable of forecasting the class transition of drought appropriately with very minimum error.

An extended transition probability matrix was extracted in which the SPI classes increased to nine while for the SHDI classes remained three (Table 12). The results in Table 12 indicate that for SPI in extreme, severe, and moderate drought states, the probability of SHDI in normal and wet states are negligible while the probability of SHDI drought state falls between 95% and 60%. This shows that a 6- to 12-month meteorological drought leads to a persistent hydrological drought, such that transition to other classes even in future months beyond the next month is not expected. As drought is a creeping phenomenon, termination of meteorological drought may not necessarily lead to immediate termination of hydrological drought. Therefore, if a prolonged moderate to extreme meteorological drought dominates the study basin, the hydrological drought could well continue after meteorological drought termination.

SPI | SPI6-SHDI1(t + 1) | SPI9-SHDI2(t + 2) | SPI12-SHDI3(t + 3) | ||||||
---|---|---|---|---|---|---|---|---|---|

Drought | Normal | Wet | Drought | Normal | Wet | Drought | Normal | Wet | |

Extremely dry | 95.32 | 2.26 | 0.36 | 92.71 | 2.54 | 0.44 | 84.96 | 6.52 | 1.44 |

Severely dry | 88.28 | 9.92 | 1.70 | 87.37 | 10.45 | 1.98 | 76.90 | 17.92 | 4.48 |

Moderately dry | 70.53 | 24.31 | 4.99 | 69.73 | 24.52 | 5.55 | 58.74 | 31.31 | 9.58 |

Slightly dry | 47.43 | 40.83 | 11.33 | 47.28 | 40.12 | 12.17 | 39.84 | 42.48 | 17.15 |

Normal | 22.49 | 51.09 | 25.37 | 22.98 | 49.65 | 26.28 | 20.89 | 47.41 | 30.67 |

Slightly wet | 11.14 | 48.84 | 38.20 | 11.72 | 47.55 | 38.88 | 11.80 | 44.89 | 41.74 |

Moderately wet | 8.35 | 45.64 | 43.76 | 8.88 | 44.57 | 44.28 | 9.29 | 42.44 | 46.41 |

Severely wet | 7.10 | 43.54 | 46.84 | 7.60 | 42.62 | 47.24 | 8.13 | 40.88 | 48.95 |

Extremely wet | 6.64 | 42.63 | 48.09 | 7.14 | 41.79 | 48.42 | 7.74 | 40.27 | 49.89 |

SPI | SPI6-SHDI1(t + 1) | SPI9-SHDI2(t + 2) | SPI12-SHDI3(t + 3) | ||||||
---|---|---|---|---|---|---|---|---|---|

Drought | Normal | Wet | Drought | Normal | Wet | Drought | Normal | Wet | |

Extremely dry | 95.32 | 2.26 | 0.36 | 92.71 | 2.54 | 0.44 | 84.96 | 6.52 | 1.44 |

Severely dry | 88.28 | 9.92 | 1.70 | 87.37 | 10.45 | 1.98 | 76.90 | 17.92 | 4.48 |

Moderately dry | 70.53 | 24.31 | 4.99 | 69.73 | 24.52 | 5.55 | 58.74 | 31.31 | 9.58 |

Slightly dry | 47.43 | 40.83 | 11.33 | 47.28 | 40.12 | 12.17 | 39.84 | 42.48 | 17.15 |

Normal | 22.49 | 51.09 | 25.37 | 22.98 | 49.65 | 26.28 | 20.89 | 47.41 | 30.67 |

Slightly wet | 11.14 | 48.84 | 38.20 | 11.72 | 47.55 | 38.88 | 11.80 | 44.89 | 41.74 |

Moderately wet | 8.35 | 45.64 | 43.76 | 8.88 | 44.57 | 44.28 | 9.29 | 42.44 | 46.41 |

Severely wet | 7.10 | 43.54 | 46.84 | 7.60 | 42.62 | 47.24 | 8.13 | 40.88 | 48.95 |

Extremely wet | 6.64 | 42.63 | 48.09 | 7.14 | 41.79 | 48.42 | 7.74 | 40.27 | 49.89 |

For other SPI classes, the probability of SHDI in normal state oscillates between 40% and 50% while the probability of SHDI in wet state is in the range of 11%–50%. In all combinations, the probability of SHDI normal state will not exceed 50% given different SPI states. It is worth noting that for all combinations, the SPI drought state will be followed by the SHDI drought state with 40%–95% probability. However, the SPI wet state will be followed by SHDI wet state with 38%–50% probability, which implies that the impact of meteorological drought state on hydrological drought state is higher than that of the wet state. This may be due to the fact that the majority of the precipitation in the basin falls as rain rather than snow. The probability of SHDI in different classes, given the observed SPI class, matches well (more than 90%) with the forecasted probabilities (not reported).

A fully extended transition probability matrix, shown in Table 13, was extracted based on Equation (17) with nine SPI and SHDI classes. According to Table 13, extreme SHDI drought state will occur when , while in other meteorological drought states, the probability of extreme hydrological drought is negligible. The probability of severe hydrological drought is greater than 10% when . The moderate and slight drought states are expected with a probability greater than 10%, even for SPI in normal condition, which indicates a lag between meteorological and hydrological droughts. The probability of normal SHDI is less than 50% in almost all combinations. The probability of extreme and severe hydrological wet state () in all SPI-SHDI combinations does not exceed 8%, which implies fast streamflow response to precipitation. Furthermore, meteorological extreme wet state cannot guarantee hydrological extreme wet state, although hydrological drought condition is not going to prevail either with about 95% of probability. The probability of SHDI slight and moderate wet states is less than 25%. It may be concluded that as the time scale of SPI increases, the probability of SHDI wet conditions increases slowly whereas the probability of drought states decreases, particularly in combination.

SPI | Extremely dry | Severely dry | Moderately dry | Slightly dry | Normal | Slightly wet | Moderately wet | Severely wet | Extremely wet |
---|---|---|---|---|---|---|---|---|---|

SHDI1(t + 1) | |||||||||

Extremely dry | 42.24 | 32.00 | 15.77 | 5.31 | 2.26 | 0.21 | 0.09 | 0.04 | 0.02 |

Severely dry | 8.07 | 27.32 | 34.05 | 18.84 | 9.92 | 1.01 | 0.41 | 0.18 | 0.10 |

Moderately dry | 1.46 | 9.94 | 28.44 | 30.68 | 24.31 | 2.94 | 1.22 | 0.53 | 0.29 |

Slightly dry | 0.31 | 2.90 | 14.58 | 29.64 | 40.83 | 6.53 | 2.85 | 1.26 | 0.70 |

Normal | 0.06 | 0.67 | 4.77 | 16.99 | 51.09 | 13.69 | 6.73 | 3.14 | 1.80 |

Slightly wet | 0.02 | 0.22 | 1.83 | 9.07 | 48.84 | 19.29 | 10.61 | 5.22 | 3.08 |

Moderately wet | 0.01 | 0.15 | 1.29 | 6.90 | 45.64 | 21.26 | 12.43 | 6.30 | 3.77 |

Severely wet | 0.01 | 0.12 | 1.06 | 5.90 | 43.54 | 22.21 | 13.46 | 6.96 | 4.21 |

Extremely wet | 0.01 | 0.11 | 0.98 | 5.54 | 42.63 | 22.56 | 13.89 | 7.24 | 4.40 |

SHDI1(t + 2) | |||||||||

Extremely dry | 44.98 | 28.09 | 14.30 | 5.33 | 2.54 | 0.26 | 0.11 | 0.05 | 0.02 |

Severely dry | 11.10 | 27.34 | 31.02 | 17.92 | 10.45 | 1.17 | 0.48 | 0.21 | 0.11 |

Moderately dry | 2.30 | 11.34 | 27.37 | 28.71 | 24.52 | 3.23 | 1.38 | 0.60 | 0.31 |

Slightly dry | 0.53 | 3.64 | 15.01 | 28.10 | 40.12 | 6.92 | 3.09 | 1.38 | 0.71 |

Normal | 0.11 | 0.90 | 5.25 | 16.72 | 49.65 | 13.98 | 7.05 | 3.33 | 1.75 |

Slightly wet | 0.04 | 0.30 | 2.12 | 9.26 | 47.55 | 19.35 | 10.89 | 5.43 | 2.92 |

Moderately wet | 0.02 | 0.21 | 1.51 | 7.14 | 44.57 | 21.21 | 12.67 | 6.50 | 3.56 |

Severely wet | 0.02 | 0.17 | 1.26 | 6.16 | 42.62 | 22.10 | 13.67 | 7.15 | 3.95 |

Extremely wet | 0.02 | 0.16 | 1.17 | 5.80 | 41.79 | 22.42 | 14.07 | 7.42 | 4.11 |

SHDI1(t + 3) | |||||||||

Extremely dry | 34.79 | 23.45 | 17.18 | 9.54 | 6.52 | 0.85 | 0.36 | 0.15 | 0.06 |

Severely dry | 11.30 | 19.56 | 25.42 | 20.61 | 17.92 | 2.62 | 1.13 | 0.48 | 0.18 |

Moderately dry | 3.39 | 9.59 | 20.37 | 25.38 | 31.31 | 5.50 | 2.45 | 1.06 | 0.40 |

Slightly dry | 1.06 | 3.97 | 12.17 | 22.64 | 42.48 | 9.57 | 4.51 | 2.01 | 0.75 |

Normal | 0.29 | 1.27 | 5.15 | 14.18 | 47.41 | 16.03 | 8.48 | 3.98 | 1.53 |

Slightly wet | 0.11 | 0.52 | 2.47 | 8.70 | 44.89 | 20.52 | 12.00 | 5.91 | 2.32 |

Moderately wet | 0.08 | 0.37 | 1.85 | 6.99 | 42.44 | 22.07 | 13.58 | 6.87 | 2.73 |

Severely wet | 0.06 | 0.31 | 1.58 | 6.18 | 40.88 | 22.81 | 14.47 | 7.43 | 2.97 |

Extremely wet | 0.06 | 0.29 | 1.49 | 5.90 | 40.27 | 23.07 | 14.80 | 7.65 | 3.07 |

SPI | Extremely dry | Severely dry | Moderately dry | Slightly dry | Normal | Slightly wet | Moderately wet | Severely wet | Extremely wet |
---|---|---|---|---|---|---|---|---|---|

SHDI1(t + 1) | |||||||||

Extremely dry | 42.24 | 32.00 | 15.77 | 5.31 | 2.26 | 0.21 | 0.09 | 0.04 | 0.02 |

Severely dry | 8.07 | 27.32 | 34.05 | 18.84 | 9.92 | 1.01 | 0.41 | 0.18 | 0.10 |

Moderately dry | 1.46 | 9.94 | 28.44 | 30.68 | 24.31 | 2.94 | 1.22 | 0.53 | 0.29 |

Slightly dry | 0.31 | 2.90 | 14.58 | 29.64 | 40.83 | 6.53 | 2.85 | 1.26 | 0.70 |

Normal | 0.06 | 0.67 | 4.77 | 16.99 | 51.09 | 13.69 | 6.73 | 3.14 | 1.80 |

Slightly wet | 0.02 | 0.22 | 1.83 | 9.07 | 48.84 | 19.29 | 10.61 | 5.22 | 3.08 |

Moderately wet | 0.01 | 0.15 | 1.29 | 6.90 | 45.64 | 21.26 | 12.43 | 6.30 | 3.77 |

Severely wet | 0.01 | 0.12 | 1.06 | 5.90 | 43.54 | 22.21 | 13.46 | 6.96 | 4.21 |

Extremely wet | 0.01 | 0.11 | 0.98 | 5.54 | 42.63 | 22.56 | 13.89 | 7.24 | 4.40 |

SHDI1(t + 2) | |||||||||

Extremely dry | 44.98 | 28.09 | 14.30 | 5.33 | 2.54 | 0.26 | 0.11 | 0.05 | 0.02 |

Severely dry | 11.10 | 27.34 | 31.02 | 17.92 | 10.45 | 1.17 | 0.48 | 0.21 | 0.11 |

Moderately dry | 2.30 | 11.34 | 27.37 | 28.71 | 24.52 | 3.23 | 1.38 | 0.60 | 0.31 |

Slightly dry | 0.53 | 3.64 | 15.01 | 28.10 | 40.12 | 6.92 | 3.09 | 1.38 | 0.71 |

Normal | 0.11 | 0.90 | 5.25 | 16.72 | 49.65 | 13.98 | 7.05 | 3.33 | 1.75 |

Slightly wet | 0.04 | 0.30 | 2.12 | 9.26 | 47.55 | 19.35 | 10.89 | 5.43 | 2.92 |

Moderately wet | 0.02 | 0.21 | 1.51 | 7.14 | 44.57 | 21.21 | 12.67 | 6.50 | 3.56 |

Severely wet | 0.02 | 0.17 | 1.26 | 6.16 | 42.62 | 22.10 | 13.67 | 7.15 | 3.95 |

Extremely wet | 0.02 | 0.16 | 1.17 | 5.80 | 41.79 | 22.42 | 14.07 | 7.42 | 4.11 |

SHDI1(t + 3) | |||||||||

Extremely dry | 34.79 | 23.45 | 17.18 | 9.54 | 6.52 | 0.85 | 0.36 | 0.15 | 0.06 |

Severely dry | 11.30 | 19.56 | 25.42 | 20.61 | 17.92 | 2.62 | 1.13 | 0.48 | 0.18 |

Moderately dry | 3.39 | 9.59 | 20.37 | 25.38 | 31.31 | 5.50 | 2.45 | 1.06 | 0.40 |

Slightly dry | 1.06 | 3.97 | 12.17 | 22.64 | 42.48 | 9.57 | 4.51 | 2.01 | 0.75 |

Normal | 0.29 | 1.27 | 5.15 | 14.18 | 47.41 | 16.03 | 8.48 | 3.98 | 1.53 |

Slightly wet | 0.11 | 0.52 | 2.47 | 8.70 | 44.89 | 20.52 | 12.00 | 5.91 | 2.32 |

Moderately wet | 0.08 | 0.37 | 1.85 | 6.99 | 42.44 | 22.07 | 13.58 | 6.87 | 2.73 |

Severely wet | 0.06 | 0.31 | 1.58 | 6.18 | 40.88 | 22.81 | 14.47 | 7.43 | 2.97 |

Extremely wet | 0.06 | 0.29 | 1.49 | 5.90 | 40.27 | 23.07 | 14.80 | 7.65 | 3.07 |

In order to compare this study with past studies there are a few reports aimed at forecasting drought class transition, only a couple of which used meteorological drought as the predictor for hydrological drought forecasting. Most of the past studies were carried out to forecast the time series of hydrological drought. For example, Li *et al.* (2016) used log-linear regression to forecast hydrological drought based on SPI and SRI. They reported that the model is capable of forecasting hydrological drought in one and two months ahead with 90% and 80–90% accuracy, respectively. Yet, our study results outperform the accuracy reported by Li *et al.* (2016) in general, while our approach is also capable of forecasting drought up to three months ahead. In other studies, to our knowledge, meteorological drought was not used as a predictor and drought classes were not forecasted. Instead, the direct value of hydrological drought index was forecasted.

## CONCLUSION

In this research, for the first time, the copula was adopted as a tool to forecast hydrological drought/wet states based on meteorological drought index as a predictor. First, SHDI and SPI at different time scales were calculated. SPI was calculated using precipitation of a number of meteorological stations. Then, using the kriging method, the time series of basin average SPI was determined. In the next step, the correlation coefficient and cross correlation between SPI and SHDI at different time scales were determined in order to select the best SPI-SHDI forecast combinations. Afterwards, the marginal distributions were determined for SPI and SHDI. Five Archimedean copulas were identified as candidates to fit the joint SPI-SHDI. The dependence parameter of these copulas was estimated using the inversion of Kendall's tau. Then, the conditional probability and the forecasted class transition matrix were calculated.

Findings show that copula functions are suitable tools for drought class transition forecasting. Also, based on the results, extreme and severe meteorological drought leads to hydrological drought with a high probability, while other classes of meteorological drought/wet conditions lead to normal hydrological state with less than 50% probability and to wet state with over 20% probability. This implies that the hydrological drought/wet state in the study basin is highly sensitive to the precipitation, such that even sustained meteorological wet condition may not guarantee hydrological wet condition occurring in the next months. It may depend on the type of precipitation in the basin. Most of the precipitation in the basin is rainfall, any variation of which in monthly rainfall can lead to immediate variation in hydrological dry/wet condition.

The findings are valuable for water resources management in the study basin due to the presence of several large hydropower dams downstream. The proposed methodology could be adopted in other basins with no specific limitation.

## REFERENCES

*138*