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 Standardized Precipitation Index (SPI) quantifies the meteorological drought, while, similar to the SPI, the Standardized Hydrological Drought Index (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 km2. 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.

Figure 1

Location of the study area.

Figure 1

Location of the study area.

Figure 2

Location of meteorological stations within Karun Basin.

Figure 2

Location of meteorological stations within Karun 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.

Dehghani 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:  
formula
 
formula
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).

Table 1

SPI and SHDI classification

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:

  • 1.

  • 2.
    For every such that and ,  
    formula

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.

Let be a joint distribution function with marginal cumulative distribution functions and . Based on Sklar's theorem (1959), there exists a copula C such that:  
formula
(1)
If and are continuous, then the copula C is unique. From Equation (1), the copula C is given by:  
formula
(2)
where and are quasi-inverses of and , respectively.
One of the most popular copula types is the Archimedean copula which is utilized for describing the dependency between the random variables and widely used in drought analysis (Maity et al. 2013; Zhang et al. 2015; Chang et al. 2016). An Archimedean copula is defined as follows:  
formula
(3)
where the function is called a generator of the copula and is a mapping which satisfies the following conditions (Genest et al. 2006):  
formula
(4)

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:

  1. Calculate the Kendall's tau (Kendall correlation coefficient) as:  
    formula
    (5)
    In an Archimedean copula with generator , this parameter is given by:  
    formula
    (6)
  2. Estimate the Kendall's tau as:  
    formula
    (7)
    where n is the number of observations and if and then the indicator function else .
  3. Suppose C be a copula function and let . Define as:  
    formula
    (8)
    For an Archimedean copula with generator , can be calculated as:  
    formula
    (9)

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

  4. as the empirical representation of is defined as follows:  
    formula
    (10)
    where are pseudo-observations defined by:  
    formula
    (11)
  5. Genest et al. (2006) proposed two statistics for testing the GOF as follows:  
    formula
    (12)
    where known as Kendall's process. Equation (12) can be expanded by a simple mathematical calculation into:  
    formula
    (13)
    and  
    formula
    (14)

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.

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

  2. Generate N random samples of size n from copula and determine the value of the test statistics and for each generated sample.

  3. 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

Due to the lag between meteorological and hydrological droughts, it is possible to forecast the hydrological drought based on history of the meteorological drought. Owing to the probabilistic characteristics of drought and the dependence between hydrological and meteorological droughts, conditional probability may be adopted for hydrological forecasting. For this purpose, SPI and SHDI time series were set up in different time scales. Suppose that SHDI and SPI are denoted by X and Y, respectively. Then, the conditional cumulative distribution function of X given is given by:  
formula
(15)
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:  
formula
(16)
where  
formula
(17)
and  
formula
(18)

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.

Table 2

Families of bivariate Archimedean copulas

Copula      
Clayton      
Frank      
Gumbel–Hougaard      
Joe      
Ali–Mikhail–Hag           
Copula      
Clayton      
Frank      
Gumbel–Hougaard      
Joe      
Ali–Mikhail–Hag           
Table 3

Characteristics of univariate cumulative distribution functions

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.

Figure 3

Framework of hydrological drought forecasting.

Figure 3

Framework of hydrological drought forecasting.

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.

Table 4

The p-value of nonparametric tests of precipitation and discharge time series

  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.

Table 5

The fitted gamma parameters to precipitation data and the results of nonparametric tests for SPI6 calculation

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 0.143 0.999 
May 12.836 20.574 −235.767 479.535 0.062 0.998 0.123 
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 0.143 0.999 
May 12.836 20.574 −235.767 479.535 0.062 0.998 0.123 
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 
Table 6

The fitted log-normal parameters to discharge data and the results of nonparametric tests for SHDI1 calculation

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.

Table 7

Matrix of correlation coefficients between SPI and SHDI in different time scales

 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 
Figure 4

Cross correlation between SPI and SHDI in different time scales.

Figure 4

Cross correlation between SPI and SHDI in different time scales.

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.

Table 8

Goodness-of-fit test based on Sn and Tn statistics

Copula Scenario    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    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.

Figure 5

Contour plot of bivariate Clayton copula and scatter plots of the joint SPI-SHDI.

Figure 5

Contour plot of bivariate Clayton copula and scatter plots of the joint SPI-SHDI.

Figure 6

The 3-D perspective plot of theoretical Clayton copula and the scatter plot of joint SPI-SHDI combinations.

Figure 6

The 3-D perspective plot of theoretical Clayton copula and the scatter plot of joint SPI-SHDI combinations.

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.

Table 9

SPI marginal distribution goodness of fit and estimated parameters

  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 
Table 10

SHDI marginal distribution goodness of fit and estimated parameters

  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 79 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 79, 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.

Figure 7

Conditional cumulative distribution function for SPI6-SHDI1(t + 1) combination.

Figure 7

Conditional cumulative distribution function for SPI6-SHDI1(t + 1) combination.

Figure 8

Conditional cumulative distribution function for SPI9-SHDI2(t + 2) combination.

Figure 8

Conditional cumulative distribution function for SPI9-SHDI2(t + 2) combination.

Figure 9

Conditional cumulative distribution function for SPI12-SHDI3(t + 3) combination.

Figure 9

Conditional cumulative distribution function for SPI12-SHDI3(t + 3) combination.

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 79. 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.

Table 11

Forecasted transition probability matrix in SPI-SHDI 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.

Table 12

Extended forecast transition probability matrix in different SPI-SHDI combinations

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.

Table 13

Class transition probability between SPI and SHDI

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

REFERENCES
Amirataee
,
B.
Montaseri
,
M.
Rezaie
,
H.
2018
Regional analysis and derivation of copula-based drought Severity-Area-Frequency curve in Lake Urmia basin, Iran
.
Journal of Environmental Management
206
,
134
144
.
Ayantobo
,
O. O.
Li
,
Y.
Song
,
S.
Javed
,
T.
Yao
,
N.
2018
Probabilistic modelling of drought events in China via 2-dimensional joint copula
.
Journal of Hydrology
559
,
373
391
.
Bowden
,
J. H.
Talgo
,
K. D.
Spero
,
T. L.
Nolte
,
C. G.
2016
Assessing the added value of dynamical downscaling using the standardized precipitation index
.
Advances in Meteorology
.
doi: 10.1155/2016/8432064
.
Chen
,
L.
Singh
,
V. P.
Guo
,
S.
Zhou
,
J.
Zhang
,
J.
2015
Copula-based method for multisite monthly and daily streamflow simulation
.
Journal of Hydrology
528
,
369
384
.
Chen
,
S.
Shin
,
J. Y.
Kim
,
T. W.
2016
Probabilistic forecasting of drought: a hidden Markov model aggregated with the RCP 8.5 precipitation projection
.
Stochastic Environmental Research and Risk Assessment
31
,
1061
1076
.
Dehghani
,
M.
Saghafian
,
B.
Nasiri Saleh
,
F.
Farokhnia
,
A.
Noori
,
R.
2014
Uncertainty analysis of streamflow drought forecast using artificial neural networks and Monte-Carlo simulation
.
International Journal of Climatology
34
,
1169
1180
.
Dehghani
,
M.
Saghafian
,
B.
Rivaz
,
F.
Khodadadi
,
A.
2017
Evaluation of dynamic regression and artificial neural networks models for real-time hydrological drought forecasting
.
Arabian Journal of Geosciences
10
,
266
.
Dracup
,
J. A.
Lee
,
K. S.
Paulson
,
E. G.
Jr.
1980
On the statistical characteristics of drought events
.
Water Resources Research
16
(
2
),
289
296
.
Edwards
,
D. C.
McKee
,
T. B.
1997
Characteristics of 20th Century Drought in the United States at Multiple Scales
.
Atmospheric Science Paper No. 634, May
, pp.
1
30
.
Fung
,
K. F.
Huang
,
Y. F.
Koo
,
C. H.
Soh
,
Y. W.
2019
Drought forecasting: a review of modelling approaches 2007–2017
.
Journal of Water and Climate Change
.
https://doi.org/10.2166/wcc.2019.236
.
Genest
,
C.
Quessy
,
J. F.
Rémillard
,
B.
2006
Goodness-of-fit procedures for copula models based on the probability integral transformation
.
Scandinavian Journal of Statistics
33
(
2
),
337
366
.
Golian
,
S.
Saghafian
,
B.
Farokhnia
,
A.
2012
Copula-based interpretation of continuous rainfall-runoff simulations of a watershed in northern Iran
.
Canadian Journal of Earth Sciences
49
,
681
691
.
Hao
,
Z.
Hao
,
F.
Singh
,
V. P.
Sun
,
A. Y.
Xia
,
Y.
2016
Probabilistic prediction of hydrologic drought using a conditional probability approach based on the meta-Gaussian model
.
Journal of Hydrology
542
,
772
780
.
http://dx.doi.org/10.1016/j.jhydrol.2016.09.048
.
Hayes
,
M. J.
Svoboda
,
M. D.
Wilhite
,
D. A.
Vanyarkho
,
O.
1999
Monitoring the 1996 drought using the standardized precipitation index
.
Bulletin of the American Meteorological Society
80
(
3
),
429
438
.
Heim
,
R. R.
2002
A review of twentieth-century drought indices used in the United States
.
Bulletin of the American Meteorological Society
83
(
8
),
1149
1166
.
Hughes
,
B.
Saunders
,
M.
2002
A drought climatology for Europe
.
International Journal of Climatology
22
,
1571
1592
.
Kao
,
S.
Govindaraju
,
R. S.
2010
A copula-based joint deficit index for droughts
.
Journal of Hydrology
380
,
121
134
.
Karthika
,
K.
Krishnaveni
,
M.
Thirunavukkarasu
,
V.
2017
Forecasting of meteorological drought using ARIMA model
.
Indian Journal of Agricultural Research
51
,
103
111
.
Kousari
,
M. R.
Hosseini
,
M. E.
Ahani
,
H.
Hakimelahi
,
H.
2017
Introducing an operational method to forecast long-term regional drought based on the application of artificial intelligence capabilities
.
Theoretical and Applied Climatology
127
,
361
380
.
Kuchment
,
L. S.
Demidov
,
V. N.
2013
Probabilistic characterization of hydrological droughts
.
Russian Meteorology and Hydrology
38
(
10
),
694
700
.
Mahmud
,
I.
Bari
,
S. H.
Ur Rahman
,
M. T.
2016
Monthly rainfall forecast of Bangladesh using autoregressive integrated moving average method
.
Environmental Engineering Research
22
,
162
168
.
Maity
,
R.
Ramadas
,
M.
Govindaraju
,
R. S.
2013
Identification of hydrologic drought triggers from hydroclimatic predictor variables
.
Water Resources Research
49
,
4476
4492
.
McKee
,
T. B.
Doesken
,
N. J.
Kleist
,
J.
1993
The relationship of drought frequency and duration to time scales
. In:
Eighth Conference on Applied Climatology
,
Anaheim, CA, USA
, pp.
179
184
.
Meng
,
L.
Ford
,
T.
Guo
,
Y.
2016
Logistic regression analysis of drought persistence in East China
.
International Journal of Climatology
37
,
1444
1455
.
Mishra
,
A. K.
Singh
,
V. P.
2010
A review of drought concepts
.
Journal of Hydrology
391
(
1–2
),
202
216
.
Mishra
,
A. K.
Singh
,
V. P.
2011
Drought modeling – a review
.
Journal of Hydrology
403
(
1–2
),
157
175
.
Montaseri
,
M.
Amirataee
,
B.
Rezaie
,
H.
2018
New approach in bivariate drought duration and severity analysis
.
Journal of Hydrology
559
,
166
181
.
Rahmat
,
S. N.
Jayasuriya
,
N.
Bhuiyan
,
M. A.
2016
Short-term droughts forecast using Markov chain model in Victoria, Australia
.
Theoretical and Applied Climatology
129
,
445
457
.
Saghafian
,
B.
Mehdikhani
,
H.
2014
Drought characterization using a new copula-based trivariate approach
.
Natural Hazards
72
(
3
),
1391
1407
.
Seibert
,
M.
Merz
,
B.
Apel
,
H.
2017
Seasonal forecasting of hydrological drought in the Limpopo Basin: a comparison of statistical methods
.
Hydrology and Earth System Sciences
21
,
1611
1629
.
Singh
,
V. P.
Strupczewski
,
W. G.
2007
Editorial (Special issues in copula)
.
Journal of Hydrological Engineering
12
,
345
.
Sklar
,
A.
1959
Fonctions de répartitions à n dimensions et leurs marges Paris (Distriburtion functions with n dimensions and their margins)
.
Publications de L'Institut Statistiques de L'Universite de Paris
8
,
229
231
.
Stagge
,
J. H.
Kohn
,
I.
Tallaksen
,
L. M.
Stahl
,
K.
2015
Modeling drought impact occurrence based on meteorological drought indices in Europe
.
Journal of Hydrology
530
,
37
50
.
Sun
,
P.
Zhang
,
Q.
Singh
,
V. P.
Xiao
,
M.
Zhang
,
X.
2016
Transitional variations and risk of hydro-meteorological droughts in the Tarim River basin, China
.
Stochastic Environmental Research and Risk Assessment
31
,
1515
1526
.
Van de Vyver
,
H.
Van den Bergh
,
J.
2018
The Gaussian copula model for the joint deficit index for droughts
.
Journal of Hydrology
561
,
987
999
.
Van Lanen
,
H. A. J.
Wanders
,
N.
Tallaksen
,
L. M.
Van Loon
,
A. F.
2013
Hydrological drought across the world: impact of climate and physical catchment structure
.
Hydrology and Earth System Sciences
17
,
1715
1732
.
Van Loon
,
A. F.
2015
Hydrological drought explained
.
Wiley Interdisciplinary Reviews: Water
2
,
359
392
.
Vezzoli
,
R.
Pecora
,
S.
Zenoni
,
E.
Tonelli
,
F.
2012
Data analysis to detect inhomogeneity, change points, trends in observations: an application to Po river discharge extremes
.
CMCC Research Paper 138
.
Wilhite
,
D. A.
Glantz
,
M. H.
1985
Understanding the drought phenomenon: definitions
.
Water International
10
,
111
120
.
Winingaard
,
J. B.
Kleink Tank
,
A. M. G.
Konnen
,
G. P.
2003
Homogeneity of 20th century European daily temperature and precipitation series
.
International Journal of Climatology
23
,
679
692
.
Wong
,
G.
van Lanen
,
H. A. J.
Torfs
,
P. J. J. F.
2013
Probabilistic analysis of hydrological drought characteristics using meteorological drought
.
Hydrological Sciences Journal
58
(
2
),
253
270
.
Wu
,
H.
Hayes
,
M. J.
Wilhite
,
D. A.
Svoboda
,
M. D.
2005
The effect of data length on the standardized precipitation index calculation
.
International Journal of Climatology
25
,
505
520
.
Wu
,
H.
Svoboda
,
M. D.
Hayes
,
M. J.
Wilhite
,
D. A.
Wen
,
F.
2007
Appropriate application of the Standardized Precipitation Index in arid locations and dry seasons
.
International Journal of Climatology
27
,
65
79
.
Yang
,
J.
Chang
,
J.
Wang
,
Y.
Li
,
Y.
Hu
,
H.
Chen
,
Y.
Huang
,
Q.
Yao
,
J.
2018
Comprehensive drought characteristics analysis based on a nonlinear multivariate drought index
.
Journal of Hydrology
557
,
651
667
.
Zhang
,
T.
Li
,
J.
Hu
,
R.
Wang
,
Y.
Feng
,
P.
2017
Drought class transition analysis through different models: a case study in North China
.
Water Science and Technology: Water Supply
17
,
138
150
.