A novel generalized combinative procedure for Multi-Scalar standardized drought Indices-The long average weighted joint aggregative criterion

Abstract

Drought hazards have complex climatic and spatio-temporal features. Therefore, its accurate monitoring is a challenging task in hydrological research. In recent, the use of standardized drought indices for drought monitoring is common in practice. However, the existence of several drought indices creates chaotic problems in data mining and decision making. This article presents a new weighting scheme for combining multiple drought indices. We propagated steady-state probabilities of Markov chain as weights in the Probabilistic Weighted Joint Aggregative Index (PWJADI) criterion. Hence, to aggregate drought characterization two or more indices, averaged long term behavior of drought classification states observed in the individual drought index is considered as a weighting characteristic. The proposed algorithm is rather general and can be used for any standardized type of drought indices. The new procedure is named as Long Averaged Weighted Joint Aggregative Criterion (LAWJAC). In this research, we focused on the three multi-scalar drought indices namely, Standardized Precipitation Index (SPI), Standardized Precipitation Evaporation Index (SPEI), Standardized Precipitation Temperature Index (SPTI). The selection of these indices is due to their similar computational procedures. In the evolution of LAWJAC, three meteorological stations of the Northern Area of Pakistan are considered. A comparison of LAWJAC is made with PWJADI. Results show significant deviations between existing and proposed methods. By the rationale of the proposed algorithm, these deviations strongly advocate the use of LAWJAC for more accuracy in drought characterization.

Keywords:

• Drought monitoring
• Markov chain
• multi-scalar drought indices
• probabilistic Weighted Joint Aggregative Index

1. Introduction

Drought is a complex natural phenomenon, and it is recurrently occurring in several regions across the globe because of climate change and an increasing trend of global warming (Hagman, 1984). Drought hazards largely affect agriculture, human life, ecology and livestock (Mechler et al., 2018). Due to continued worst climate conditions, the need for climate control and water resource management is to uncover hidden challenges of future drought conditions. For example, several regions of Pakistan are bearing a drought condition (Memon et al., 2018). Only in 2010, several peoples in various regions have suffered due to the severe effects of drought (Guha-Sapir et al., 2011). Hence accurate drought monitoring and characterization became challenging in hydrological research (Wilhite, 2006; Wilhite et al., 2014).

Numerous drought monitoring tools and methods are developed. Precipitation and temperature are key time series data in most of these methods. Initially, the Palmer Drought Severity Index (PDSI) was developed to identify dry and wet conditions at a meteorological station (Palmer, 1965). Several studies have used PDSI as a drought monitoring indicator in various regions. For example, Eder et al., (1987) examined spatio-temporal variability in drought using PSDI at South-Eastern United States. Dai et al. (2004) explored the relationship between soil moisture and surface warming based on PDSI. Other applications mainly include Klugman (1978), Bhalme and Mooley (1979), Karl and Quayle (1981), Karl and Koscielny (1982) and Rind et al. (1990). Later studies have reshaped PDSI by including new variables. For example, Yu et al. (2019) proposed a new drought index- the Modified Palmer Drought Severity Index (MPDSI). In MPDSI, the PDSI is improved by adding irrigation quotas and soil water deficit. Shen et al. (2019) developed the Integrated Drought Condition Index (IDCI) by combining precipitation, soil moisture, potential evaporation, temperature, and vegetation conditions.

After PDSI, McKee et al. (1993) introduced a new probabilistic procedure for drought index: the Standardized Precipitation Index (SPI). In SPI procedure, precipitation data is used to obtain quantitative value of the Cumulative Distribution Function (CDF) of the gamma distribution. The CDF is than standardized by appropriate approximation. Guttman (1999) has recommended SPI index as a standard method for all users and for all regions. Many authors have used SPI drought index in various regions. Some of them are Tsakiris and Vangelis (2004), Moreira et al. (2008), Nalbantis and Tsakiris (2009), Zhang, et al. (2012), etc. In later research, several authors have introduced various standardized procedures for monitoring drought. For example, Vicente-Serrano et al. (2010) and Ali et al. (2017) have introduced two drought indices namely Standardized Precipitation Evapotranspiration Index (SPEI) and Standardized Precipitation Temperature Index (SPTI). All these standardized indices have the similar mathematical procedure and drought classification criterion. Erhardt and Czado (2018) acknowledged all these types of indices as a Standardized Drought Indices (SDI). The main advantages of standardized drought indices are to monitor and compare different regions with different time scales (Payab and Türker, 2019).

However, the main problem in SDI-type drought indices is their probabilistic estimation. That is, the subjective choice of probability distribution greatly affects the operational definition of drought. Quiring (2009b) has concluded that probabilistic estimation-based values of drought indices (i.e SPI) are more sensitive to the choice of distribution. An alternative to probabilistic estimation, Farahmand and AghaKouchak (2015) have introduced a generalized procedure of non-parametric estimation of drought indices. They used Gingorten (Gringorten, 1963) probability position formula as an alternative to CDF of the gamma distribution. One problem in non-parametric estimation is the compatibility of probability plotting position formulas with probability distributions. Further, the use of one probability plotting position formula is not enough for all the time scales of drought indices and geographical regions as well. For example, analogous to the law of large numbers, the distribution for the time series data of higher time scales have a tendency to converge into the normal distribution. Moreover, instead of using only gamma distribution, the behavior of data can be captured by more advance and appropriate probability distributions Stagge et al. (2015). Equivalent to this, the compatibility of various plotting position formula such as Hanzen (Allen, 1914) probability plotting position formula with various probability distributions has great importance in accurate propagation of uncertainty. For instance, several authors have suggested the deployment of various probabilities plotting formula for specified probability distribution using correlation (Vogel, 1986; Mehdi and Mehdi, 2011).

In summation, the accuracy of drought characterization is badly affected by the existence of a certain amount of errors (Quiring, 2009a). Hao and Singh, 2015; It is witnessed that various methods of SDI type drought indices produce varying and contrary results in different climatological regions. For example, Tan et al. (2015) reported that SPEI is more applicable than SPI in the country of Ningxia. Previous research proven that the inconsistencies and deviations among drought indices are due to inappropriate selection of probability distributions (Stagge et al. 2015), subjective choice of probability plotting position formula (Looney and Gulledge, 1985; El-Shanshoury and El-Hemamy, 2013), geographical characteristics, regional compatibility to the choice of drought index, and the complex feature of regional climate (Ali et al., 2017).

However, to make accurate and efficient drought monitoring, several authors have recommended the use of multiple or combined drought indices for a certain region. In previous research, various authors have provided various combinatorial techniques. For example, Balint et al. (2013) proposed Combined Drought Index (CDI). CDI consists of the weighted average of the Vegetation Drought Index (VDI), 2-decked lagged Precipitation Drought Index (PDI) and Temperature Drought Index (TDI). Accordingly, Luetkemeier et al. (2017) have proposed Blended Drought Index (BDI) by integrating SPI, SPEI, and Standardized Soil Moisture Index (SSMI) in copula equation (Hao, and AghaKouchak, 2013). Recently, Ali et al. (2019a) have provided a Probabilistic Weighted Joint Aggregative Drought Index (PWJADI) criterion for combining temporal classification of three drought indices namely SPI, SPEI, and SPTI.

This paper proposes a new criterion to combine temporal classification of drought determined by various standardized tools. The proposed criterion is more likely to PWJADI (Ali et al., 2019a). However, its weighting scheme differs from those which have been used in the PWJADI criterion. In this research, long term behavior of drought classification states is quantified using steady-state probabilities, which are then used as a weighting characteristic in aggregation. We named the new criterion-the Long Averaged Weighting Joint Aggregative Criterion (LAWJAC). The rationale behind the use of steady-state probabilities as weights and procedure of accumulation presented in Section 3. The application of the proposed criterion is presented for three meteorological stations located in the Northern part of Pakistan. We used three drought indices, namely SPI, SPEI and SPTI at 1, 3, 6, 9, 12, and 24-month time scales. The performance of LAWJAC is assessed by comparing its results with PWJADI. The distribution of this paper is as follows: Section 2 consists of methods and existing literature on the relevant subject, where a detailed overview and brief description of the Markov chain, steady-state probabilities, SDI-type and PWJADI criterion is presented. The proposed method is presented in Section 3. Section 4 consists of the application of the proposed method. Results are summarized in Section 5, while discussion and conclusion are described in Sections 6 and 7, respectively.

2. Methods

2.1. Stochastic process, discrete Markov chain, and steady-states probabilities

A stochastic process with state space $E\in (e_1,e_2,\dots ,e_n)$ is a collection $\{z_t;t\in T\}$ of a family random variables $z_t$ defined on the same probability space and taking values in state space $E$ (Lawler, 2018). When the state space $E$ of the stochastic process $z_t$ is continuous, it is called a Markov process regardless of whether the parameter (or time) is discrete or continuous, and when the Markov process is discrete-valued (i.e. discrete state space $E$) , it is called a Markov chain.

Markov chain helps to predict the state of the process of any uncertain phenomena by past transient behavior of the process. Markov chain has wide range of applications in hydrology and related research (Akyuz et al., 2012). Sharma and Panu (2012) have concluded that Markov chain models are efficient tools for predicting drought characteristics using a standardized hydrological index. Mallya et al., (2013) have utilized the hidden Markov chain model on the temporal classification of drought states. They reported that hidden Markov chain models are appropriate for analyzing spatio-temporal characteristics of drought. Avilés et al., (2016) assessed the suitability of the Markov chain-based model and Bayesian network-based model for the forecasting of drought classification states. Sedlmeier et al., (2016) have analyzed climate dynamics of extreme events using the Markov chain model. Recently, Chen et al. (2019) proposed a hidden Markov model framework for forecasting categorical sequence of drought states.

A short description of the Markov chain process, transition probability matrices and steady-states probabilities are as follows:

Let $S=\mathrm{}\{s_1,\mathrm{}{s}_2,\dots \dots \mathrm{}{s}_r\}\mathrm{}$ be the states of the process, where the process may begin in one of these states and successively moves from one state to another. Here each of the next movements of the process is called a step. If the current position of the chain is in state$\mathrm{}{s}_i,$ then it passes to state $s_j$ with probability $p_{ij}$ for the next step. In the Markov chain, the probability of the next state is independent of the probability of a previous state (Runger and Wasserman, 1979; Poggi et al., 2000). Mathematically, the probability of moving one state to another state is represented by the Transition Probability Matrix (TPM). TPM quantifies “$n$” step transient behavior of the process. It is a square matrix, where the elements within the matrix are non-negative and real. The matrix can be represented as follows: $$P_{ij}^{(n)}=\left[\begin{array}{ccc}{p}_{11}& \cdots & p_{1n}\\ ⋮& ⋮\cdots & ⋮\\ p_{n1}& \cdots & p_{nn}\end{array}\right]$$

The elements within the above matrix quantify the transient probabilities of the state space $E$ of the process. And the sum of the row element (probabilities) is equal to 1.

However, the long-term behavior of the process states is quantified by the stationary probabilities or limiting distribution of the process. These probabilities are often called steady-state probabilities.

In another way, the average probabilities that the system will be in a certain state after a large number of transition periods are called steady state probabilities. The system does not need to remain in one condition. The system will remain to move from state to state in future periods; however, the average probability of moving from state to state for all periods will remain constant in the long run. In a Markov process, the probabilities will approach steady-state after some number of steps.

Let ${\pi }_j$ denotes the limiting probability of $\mathit{ith}$ after “$n$” step, the formula of steady-states probability is defined as follows. (1) $$\pi (j)=\underset{n\to \mathrm{\infty }}{\mathit{lim}}P(Z_n-j|Z_o-i)$$(1)

In other words, (2) $${\pi }_j=\underset{n\to \infty }{\lim }{P}_{ij}^{(n)}$$(2)

Detailed theory and mathematical description of steady-state probabilities of Markov chain can be found in Stewart, (2009). In this research, we have focused on the use of long-term behavior of the hydrological process as weights in updating and correcting information of various sources. In our proposal, we configured steady-state probabilities as weights in the aggregation of various information.

2.2. Standardized drought indices

In literature, many authors have provided various drought indicators for a standardized procedure of drought indices. Some of them are available in Svoboda et al. (2016). However, Standardized Drought Indices (SDI) procedures are the most commonly used and acceptable around the world (Erhardt and Czado, 2018). The important characteristics of SDI methods are that they are comparable.

McKee et al. (1993) proposed the Standardized Precipitation Index (SPI) to quantify drought characteristics. SPI uses time series data of precipitation ($P_i$) at a given location. The first standardization of SPI method involves Gamma distribution. The disadvantage of SPI drought index is that it uses only a single climatic variable. However, drought is a complex hazard that should be characterized using multiple climatic variables such as low/high temperature, relative humidity, and wind speed, etc.

Later, Vicente-Serrano et al. (2010) introduced another standardized index- the SPEI. SPEI has both the characteristics of SDI- the standardization of values and time scales. On the same methodological structure of SPI, SPEI uses a water balance equation (see equation (3)(3) $$DE{F}_i=P_i-PE{T}_{i\mathrm{}}$$(3) ). (3) $$DE{F}_i=P_i-PE{T}_{i\mathrm{}}$$(3)

In the above equation, $P_i$ is the monthly total amount of precipitation, ${\mathit{PET}}_i$ is the estimated amount of Potential Evapotranspiration (PET) and DEF_i denotes the difference between $P_i$ and ${\mathit{PET}}_i$. One of the main problem in SPEI index is its operational suitability and compatibility. For example, the original proposal of SPEI suggests Thornthwaite (Th) (Thornthwaite, 1948) equation for the estimation of potential evaporation. However, Th equation gives unreliable and inconsistent values of potential evaporation at low temperatures and arid regions (Jensen et al., 1990).

Recently, Ali et al. (2017) proposed another drought index-the SPTI. Contrary to SPEI, SPTI accounts for the direct role of temperature. In SPTI index, De Martone Aridity Index (DMAI) is standardized on the same step of SPI and SPEI procedures. Here, De Martone Aridity Index is computed using monthly total precipitation and average monthly temperature for the selected meteorological station (see equation 4(4) $$\text{DM}{\text{AI}}_i=\frac{P_i}{10+T_i}$$(4) ). (4) $$\text{DM}{\text{AI}}_i=\frac{P_i}{10+T_i}$$(4)

In the above equation, $P_i$ is the monthly total precipitation and, $T_i$ is the mean monthly temperature.

Further, the detailed computational procedure of SPI, and SPEI, SPTI, are provided in Section 3.

3. The combinative procedure for drought Characterization-The proposed framework

This section describes four phases of the proposed criterion for the aggregation of multiple drought indices (see Figures 1 and 2).

Fig. 1. Flowchart of the proposed criterion.

Fig. 2. Flowchart of the LAWJADC equation.

A detailed description of each phase is as follows:

3.1. Phase: 1 selection of drought indices

This phase consists of the selection of drought indicators from the list of all available drought indicators of the SDI procedure and the estimation procedures.

The major concern of this phase is to select the climatic parameters and the time scale for the estimation of multi-scalar drought indices. Depending on the nature of climatic, soil type and tropical status, various drought indices required various climatic parameters such as temperature, precipitation, solar radiation, and humidity, etc. Therefore, an optimized selection of drought indices and their estimation procedure can significantly contribute to accurate and reliable drought monitoring. In particular, this step requires a deep knowledge of the following issues:

• The identification of the nature of the gauging station and the accessibility of the time series data on the climatic parameters.

• The appropriate selection of multi-scalar drought indicator (i.e. SPI, SPEI, SPTI) that can be accomplished with the available data.

• Type of drought with their corresponding time scale. In this step, the time scale of multi-scalar drought indices is selected. For example, short time scales are recommended for meteorological, whereas longer time scale is specified for the monitoring of agricultural and hydrological drought.

Keeping the above points, this research comprised of the three multi-scaler drought indices named, SPI, SPEI, and SPTI. The selection of these indices is due to their similar mathematical procedure and drought classification criterion and the regional compatibility of data that we have considered in the application section. Important details and applications of these indices have been described in Section 2.2.

3.2. Phase: 2. Estimation of standardized values and markov chain theory

After the selection of drought indicators, the next step is to estimate standardized values using appropriate estimation methods. In the very first proposal of SPI index, McKee et al. (1993) have used the CDF of gamma distribution for obtaining standardized values. In a later study, various authors recommended varying distribution for the various climatic regions. Recently, Stagge et al. (2015) have introduced a compact of using varying distribution for various drought indices at their individual time scales. Parallel to probability distribution-based standardization, Farahmand and AghaKouchak (2015) have introduced a non-parametric way of obtaining standardized ways. In this paper, we used Stagge et al. (2015) guidelines for standardization. The procedure of standardization is given.

Let $DA{I}_i\in ({\mathrm{P}}_i,{\mathrm{}\text{DEF}}_i{,\mathrm{}\text{DMT}}_i\mathrm{)}\mathrm{}$ be time series data of the selected Drought Assessment Indicators (DAI). To obtain the standardized time series, the stepwise procedure is as follows.

• First step is to reveal the appropriate probability distribution for each time series. This step could be done by using the Kolmogorov test (Massey, 1951) and Anderson darling test (Scholz and Stephens, 1987).

• In the second step, those probability distributions that have a minimum value of Akaike Information Criterion (AIC) (Yamaoka et al., 1978) or Bayesian Information Criterion (BIC) (Watanabe, 2013) are recommended for standardization.

• The third step consists of the standardization process. In this step, a little modification is made to adjust null/undefined values. For example, equation (5)(5) $$G\left(DA{I}_i\right)=q+\left(1-q\right)F\left(DA{I}_i\right)$$(5) suggests the little modification in the CDF of the gamma distribution. (5) $$G\left(DA{I}_i\right)=q+\left(1-q\right)F\left(DA{I}_i\right)$$(5)

  The above modification is suggested due to zero values in the numerical vectors of DAI. Here, $q$ indicates the proportion of undefined values in this vector.

• In the fourth step, SDI values are obtained by standardizing the modified CDF under the following approximation procedure. (6) $$\mathit{SDI}=-(k+\frac{a_0+\mathrm{}{a}_{1}k+\mathrm{}{a}_2{k}^2}{1+\mathrm{}{b}_0+\mathrm{}{b}_{1}k+\mathrm{}{b}_2{k}^2+\mathrm{}{b}_3{k}^3})$$(6)

For (7) $$k=\sqrt{\hspace{0.17em}\text{ln}\hspace{0.17em}\left[\frac{1}{{\left\{G\left(DA{I}_i\right)\right\}}^2}\right]}$$(7)

When (8) $$0<G\left(\mathit{DAI}\right)<0.5$$(8) (9) $$\mathit{SDI}=+(k+\frac{a_0+\mathrm{}{a}_{1}k+\mathrm{}{a}_2{k}^2}{1+\mathrm{}{b}_0+\mathrm{}{b}_{1}k+\mathrm{}{b}_2{k}^2+\mathrm{}{b}_3{k}^3})$$(9)

And for (10) $$k=\sqrt{\hspace{0.17em}\text{ln}\hspace{0.17em}\left[\frac{1}{{\left\{1-G\left(DA{I}_i\right)\right\}}^2}\right]}$$(10)

When (11) $$0.5\le G(DA{I}_i)\le 1$$(11)

After the estimation and classification of time series data of drought indices, we have incorporated the theory of the Markov chain and its applications in hydrology.

Recall, a Markov chain process $Zt$ is the random variable where the probability of each state for all time points of the process is described by the following equation: (12) $$P(X_{t+1}=s|X_t=s_t,\dots .,X_o=s_o)=P(X_{t+1}=s|X_t=s_t)$$(12)

For our proposed method, this phase reflects the classification of drought states as the Markov chain process.

Consider, in general, $C$₁, $C$₂, … $Cn$ be the drought classification states of SDI type processes (for our cases we have seven drought classes, see table 1). One can imagine the time series data of various drought classes provided in table 1 as the Markov chain process. That is, we have considered that $\mathit{SDIt}$ (i.e SPI, SPEI and SPTI) is qualitative time series data of drought classes as discrete Markov process. Then the transition probability matrix of the process $\mathit{SDIt}$ can be defined as follows: $$\begin{array}{c}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{C}1\mathrm{}\mathrm{}\mathrm{}\mathrm{C}2\mathrm{}\mathrm{}.\mathrm{}\mathrm{}.\mathrm{}\mathrm{}.\mathrm{}\mathrm{}.\mathrm{}\mathrm{}\mathrm{C}7\\ {(P_{ij})}^t=\begin{array}{c}C1\\ C2\\ .\\ .\\ .\\ .\\ C7\end{array}\left|\begin{array}{ccccccc}{p}_{11}& p_{12}& .& .& .& .& p_{17}\\ p_{21}& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ p_{71}& .& .& .& .& .& p_{77}\end{array}\right|\end{array}$$

Table 1. Drought classification criterion.

Range of Index Values	Class
$\mathit{SDI}\ge 2$	Extreme Wet (EW)
$1.5\le \mathit{SDI}<2$	Severe Wet (SW)
$1\le \mathit{SDI}<1.5$	Moderate Wet (MW)
$-1\le \mathit{SDI}<1$	Normal Dry (ND) (Near Normal)
$-1>\mathit{SDI}\ge -1.5$	Moderate Dry (MD)
$-1.5>\mathit{SDI}\ge -2$	Severe Dry (SD)
$\mathit{SDI}<-2$	Extreme Dry (ED)

The above (transition probability) matrix, either first order on n^th order, gives a probability of moving one drought states to another drought states. In the previous study, Ali et al. (2019a) have suggested these probabilities as switching weights. They assumed that the behavior of the process follows first-order Markov chain. Contrary to one step transient probability, this research suggests accounting the long-run behavior of the chain. Consequently, we have proposed steady-state probabilities as weights in the aggregation criterion.

We have presented the theory and application of steady-state probability in section 2.1. Accordingly, the transition probability matrix of the corresponding drought indices can be written as follows: $$\begin{array}{c}\begin{array}{c}\begin{array}{c}\text{ED}\mathrm{}\text{SD}\mathrm{}\text{MD}\mathrm{}\text{ND}\mathrm{}\text{MW}\text{SW}\mathrm{}\text{EW}\\ TPM\mathrm{}\text{of}\mathrm{}\mathit{SPI}=\begin{array}{c}ED\\ \begin{array}{c}SD\\ MD\\ ND\\ MW\end{array}\\ SW\\ EW\end{array}\left[\begin{array}{cccc}{x}_{11}& x_{12}& x_{13}& x_{14}\\ x_{21}& x_{22}& x_{23}& x_{24}\\ x_{31}& x_{32}& x_{33}& x_{34}\\ \begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_{41}\\ x_{51}\end{array}\\ x_{61}\end{array}\\ x_{71}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_{42}\\ x_{52}\end{array}\\ x_{62}\end{array}\\ x_{712}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_{43}\\ x_{53}\end{array}\\ x_{63}\end{array}\\ x_{73}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{x}_{44}\\ x_{54}\end{array}\\ x_{64}\end{array}\\ x_{74}\end{array}\end{array}\mathrm{}\begin{array}{c}{x}_{15}\\ x_{25}\\ x_{35}\\ \begin{array}{c}\begin{array}{c}{x}_{45}\\ x_{55}\end{array}\\ x_{65}\\ x_{75}\end{array}\end{array}\mathrm{}\begin{array}{c}{x}_{16}\\ x_{26}\\ x_{36}\\ \begin{array}{c}{x}_{46}\\ x_{56}\\ x_{66}\\ x_{76}\end{array}\end{array}\mathrm{}\begin{array}{c}{x}_{17}\\ x_{27}\\ x_{37}\\ \begin{array}{c}{x}_{47}\\ x_{57}\\ \begin{array}{c}{x}_{67}\\ x_{77}\end{array}\end{array}\end{array}\right]\end{array}\end{array}\\ \text{ED}\mathrm{}\text{SD}\mathrm{}\text{MD}\mathrm{}\text{ND}\mathrm{}\text{MW}\text{SW}\mathrm{}\text{EW}\\ TPM\mathrm{}\text{of}\mathrm{}\mathit{SPEI}=\begin{array}{c}ED\\ \begin{array}{c}SD\\ MD\\ ND\\ MW\end{array}\\ SW\\ EW\end{array}\left[\begin{array}{cccc}{y}_{11}& y_{12}& y_{13}& y_{14}\\ y_{21}& y_{22}& y_{23}& y_{24}\\ y_{31}& y_{32}& y_{33}& y_{34}\\ \begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_{41}\\ y_{51}\end{array}\\ y_{61}\end{array}\\ y_{71}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_{42}\\ y_{52}\end{array}\\ y_{62}\end{array}\\ y_{712}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_{43}\\ y_{53}\end{array}\\ y_{63}\end{array}\\ y_{73}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_{44}\\ y_{54}\end{array}\\ y_{64}\end{array}\\ y_{74}\end{array}\end{array}\mathrm{}\begin{array}{c}{y}_{15}\\ y_{25}\\ y_{35}\\ \begin{array}{c}\begin{array}{c}{y}_{45}\\ y_{55}\end{array}\\ y_{65}\\ y_{75}\end{array}\end{array}\mathrm{}\begin{array}{c}{y}_{16}\\ y_{26}\\ y_{36}\\ \begin{array}{c}{y}_{46}\\ y_{56}\\ y_{66}\\ y_{76}\end{array}\end{array}\mathrm{}\begin{array}{c}{y}_{17}\\ y_{27}\\ y_{37}\\ \begin{array}{c}{y}_{47}\\ y_{57}\\ \begin{array}{c}{y}_{67}\\ y_{77}\end{array}\end{array}\end{array}\right]\\ \text{ED}\mathrm{}\text{SD}\mathrm{}\text{MD}\mathrm{}\text{ND}\mathrm{}\text{MW}\text{SW}\mathrm{}\text{EW}\\ TPM\mathrm{}\text{of}\mathrm{}\mathit{SPTI}=\begin{array}{c}ED\\ \begin{array}{c}SD\\ MD\\ ND\\ MW\end{array}\\ SW\\ EW\end{array}\left[\begin{array}{cccc}{u}_{11}& u_{12}& u_{13}& u_{14}\\ u_{21}& u_{22}& u_{23}& u_{24}\\ u_{31}& u_{32}& u_{33}& u_{34}\\ \begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_{41}\\ u_{51}\end{array}\\ u_{61}\end{array}\\ u_{71}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_{42}\\ u_{52}\end{array}\\ u_{62}\end{array}\\ u_{712}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_{43}\\ u_{53}\end{array}\\ u_{63}\end{array}\\ u_{73}\end{array}& \begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_{44}\\ u_{54}\end{array}\\ u_{64}\end{array}\\ u_{74}\end{array}\end{array}\mathrm{}\begin{array}{c}{u}_{15}\\ u_{25}\\ u_{35}\\ \begin{array}{c}\begin{array}{c}{u}_{45}\\ u_{55}\end{array}\\ u_{65}\\ u_{75}\end{array}\end{array}\mathrm{}\begin{array}{c}{u}_{16}\\ u_{26}\\ u_{36}\\ \begin{array}{c}{u}_{46}\\ u_{56}\\ u_{66}\\ u_{76}\end{array}\end{array}\mathrm{}\begin{array}{c}{u}_{17}\\ u_{27}\\ u_{37}\\ \begin{array}{c}{u}_{47}\\ u_{57}\\ \begin{array}{c}{u}_{67}\\ u_{77}\end{array}\end{array}\end{array}\right]\end{array}$$

Hence, the limiting probability of each state in each index is 1 × 7 row vector denoted by the following expressions: (13) $${\prod }_i\left(\mathit{SPI}\right)=\left[\begin{array}{cccc}{\prod }_1\left(ED\right),& {\prod }_2\left(SD\right),& {\prod }_3\left(MD\right),& \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{\prod }_4\left(ND\right),& {\prod }_5\left(MW\right)\end{array},& {\prod }_6\left(SW\right),\end{array}& {\prod }_7\left(EW\right)\end{array}\end{array}\right]$$(13) (14) $${\prod }_i\left(\mathit{SPEI}\right)=\left[\begin{array}{cccc}{\prod }_1\left(ED\right),& {\prod }_2\left(SD\right),& {\prod }_3\left(MD\right),& \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{\prod }_4\left(ND\right),& {\prod }_5\left(MW\right)\end{array},& {\prod }_6\left(SW\right)\end{array},& {\prod }_7\left(EW\right)\end{array}\end{array}\right]$$(14) (15) $${\prod }_i\left(\mathit{SPTI}\right)=\left[\begin{array}{cccc}{\prod }_1\left(ED\right),& {\prod }_2\left(SD\right),& {\prod }_3\left(MD\right),& \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{\prod }_4\left(ND\right),& {\prod }_5\left(MW\right)\end{array},& {\prod }_6\left(SW\right)\end{array},& {\prod }_7\left(EW\right)\end{array}\end{array}\right]$$(15)

These vectors (steady-states probabilities) are nothing but the long-term behavior of drought states.

On the same rationale of PWJADI, in this paper steady-states probability is utilized as a weight. However, the significance of the proposed model over PWJADI (Ali et al., 2019a) is that the proposed model is free from the assumption of 1^st order Markov chain. The next phase describes how these probabilities are configured in the aggregation model.

3.3. Phase 3 and Phase 4 the proposed model of weighting scheme for aggregation criterion- The long-averaged weighted joint aggregation criterion (LAWJAC)

In this section, we described phase 3 and phase 4. Consider the vectors of stationary distributions of drought classes in each index denoted by $\prod _i(\mathit{SPI}),$$\prod _i(\mathit{SPEI}),$ $\prod _i(\mathit{SPTI}).$ Here, each vector of stationary distributions describes the averaged long-term probability (or proportion) of drought classes in each index. That is, the steady-state probability of drought class corresponding to the drought index describes the visit of drought class in the long term. Hence, to adjust the errors of inaccurate determination of drought classes, and to aggregate the decisions, we proposed to define and rely on those drought categories that have a larger value of the corresponding stationary distribution of drought classes. Equation (16)(16) $$\mathit{LAWJAC}=\{\begin{array}{l}\mathit{SPI}\mathrm{}\mathrm{}\text{if}\mathrm{}\mathrm{}\prod _i(\mathit{SPI})>\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPEI}\mathrm{}\text{if}\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPTI},\mathrm{}\text{other}\mathrm{}\text{wise}\end{array}$$(16) presents the mathematical form of the aggregative model of SPI, SPEI, and SPTI. (16) $$\mathit{LAWJAC}=\{\begin{array}{l}\mathit{SPI}\mathrm{}\mathrm{}\text{if}\mathrm{}\mathrm{}\prod _i(\mathit{SPI})>\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPEI}\mathrm{}\text{if}\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPTI},\mathrm{}\text{other}\mathrm{}\text{wise}\end{array}$$(16)

At a single time point, the interpretation of the proposed aggregative model presented in equation (2)(2) $${\pi }_j=\underset{n\to \infty }{\lim }{P}_{ij}^{(n)}$$(2) is straightforward. For instance, if a station is identified in a Normal Dry condition by SPI index, while SPEI and SPTI indicate Severe Wet and Severe Dry conditions, respectively. Then by the systems of equation (16)(16) $$\mathit{LAWJAC}=\{\begin{array}{l}\mathit{SPI}\mathrm{}\mathrm{}\text{if}\mathrm{}\mathrm{}\prod _i(\mathit{SPI})>\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPEI}\mathrm{}\text{if}\prod _i(\mathit{SPEI})>\prod _i(\mathit{SPTI})\\ \mathit{SPTI},\mathrm{}\text{other}\mathrm{}\text{wise}\end{array}$$(16) , the choice of drought class is made for that drought index which has a higher value of corresponding steady-state (stationary distribution) probability. More loosely speaking, the system of LAWJAC selects those drought classes which have a greater value of average long-run probability (proportion) in a specific month. Accordingly, the assignment of weights is suggested for whole time series data of SPI, SPEI, and SPTI. And those drought classes which receive maximum values of weights are selected as an aggregative class. We called the new aggregative vectors of drought classes and quantitative values as a LAWJAC.

4. Comparison

In this research, comparative assessment and suitability of the proposed methods are addressed by including the PWJADI criterion. A brief note on PWJADI criterion is as follows:

Recently, Ali et al. (2019a) proposed a new drought aggregation criterion- the Probabilistic Weighted Joint Aggregative Drought Index (PWJADI). In PWJADI, the vectors of three commonly used drought indices namely SPI, SPEI, and SPTI are combined in a systematic way (see equation 17(17) $$\mathit{PWJADI}=\max (a_{ij}\in \mathit{SPI};b_{ij}\in \mathit{SPEI};c_{ij}\in \mathit{SPI})$$(17) ). (17) $$\mathit{PWJADI}=\max (a_{ij}\in \mathit{SPI};b_{ij}\in \mathit{SPEI};c_{ij}\in \mathit{SPI})$$(17)

In the above equation, $a_{ij},$ $b_{ij},$ and $c_{ij}$ are the switching probabilities. The probabilities are defined from the transition probability matric of the corresponding index (see Ali et al., 2019a). In the PWJADI procedure, transition probabilities of Markov Chain are used as a switching weight. These switching weights are extracted from the first-order transition probability matrix. In the literature, numerous studies are based on the second or n^th order Markov chain (Khadr, 2016). For example, Sen (1990) have used a second-order Markov chain for the characterization of dry and wet states. Lennartsson et al. (2008) reported that first-order Markov chain is inadequate, hence, a multiple-step Markov chain is required. Recently, Avilés et al. (2016) used the first- and second-order Markov chains for the prediction of the Drought Index (DI). Considering this literature, the constraint of first-order Markov condition reduces the scope of PWJADI. Therefore, the proposed method of weighting scheme is an alternative to the weighting scheme of PWJADI.

5. Application

Pakistan is listed in the three most water-stressed countries of the globe (Farooqi et al., 2005). In Pakistan, the annual average amount of precipitation falls below 250 mm. Therefore, water scarcity and recurrent occurrence of drought hazards are the main problems of the country (El Kharraz et al., 2012; Lal, 2018). Recently due to severe drought, several human deaths were occurred in the Tharpaker district (Rana, and Naim, 2014). Agricultural and livestock sectors are badly destroyed by the severe drought. In addition, a shortage of water has created crises in hydropower energy. Consequently, the country is bearing severe economic crises and a decline in Gross Domestic Product (GDP). These crises are occurring due to the big difference between demand and supply of energy. Although the construction of many water reservoirs is undergoing. However, due to global warming and climate change, Pakistan is facing many other challenges related to the current water reservoirs as well.

The proposed criterion is tested at three different meteorological stations located in the Northern region of Pakistan. The role of the Northern areas of Pakistan has significant importance in the overall country’s climatology. The location of this region lies in high precipitation clusters. From a drought monitoring perspective, these regions are crucial for advance water management. Due to high altitude, the rainfall in these regions affects overall country climate, streamflow and surface runoff. Figure 3 shows the location of the selected stations. For this research, time series data of precipitation and temperature are collected from Karachi Data Processing Center (KDPC http://www.pmd.gov.pk/rmc/RMCK/Services_Climatology.html). The data sets fulfill standard requirements of the World Meteorological Organization (WMO). Before dispatching data to us, issues related to the tabulation, removal of errors, adjusting outliers, and missing values are done by KDPC themselves. This dataset has been cited in our recent publications (see Ali, et al., 2019a, 2019b, 2017). Table 2 shows the summary statistic of the selected stations.

Fig. 3. Geographical locations of selected meteorological stations of Pakistan.

Table 2. Climatology of the selected stations.

Station	Months	Precipitation		Temperature (max)		Temperature (min)
		Mean	Kurtosis	Mean	Kurtosis	Mean	Kurtosis
Astore	1	40.079	0.994	3.033	−0.616	−7.157	−0.165
	2	45.545	0.779	4.418	0.805	−5.185	0.610
	3	73.862	1.558	9.097	−0.344	−0.812	0.106
	4	79.423	0.641	15.472	0.202	4.083	−0.191
	5	65.636	0.667	20.207	−0.172	7.784	0.089
	6	24.8	−0.772	24.647	−0.274	11.419	1.348
	7	25.266	2.243	27.123	−0.625	14.505	0.287
	8	25.838	0.544	26.769	−0.449	14.497	0.204
	9	25.551	15.810	23.251	−0.951	10.216	−0.271
	10	20.923	2.655	17.834	0.776	4.632	0.809
	11	17.192	2.153	11.631	0.103	−0.283	−0.559
	12	27.681	8.709	5.734	−0.111	−4.397	−0.559
Bunji	1	6.233	4.481	10.041	0.300	0.136	−0.602
	2	9.592	8.260	12.776	−0.057	2.892	0.548
	3	14.102	7.409	18.474	−0.545	7.786	−0.674
	4	24.956	2.342	24.467	−0.199	12.269	−0.380
	5	27.679	1.936	28.697	0.006	15.316	−0.152
	6	10.6	0.926	33.32	−0.313	19.2	1.050
	7	17.875	3.600	35.386	0.366	22.603	−0.210
	8	22.434	3.233	34.91	−0.091	21.733	0.429
	9	14.575	5.695	31.455	−0.245	16.674	1.986
	10	7.436	24.935	25.678	−0.092	10.592	−0.094
	11	2.941	9.026	18.751	−0.711	4.696	2.095
	12	4.945	12.467	12.179	0.427	1.156	0.007
Gupis	1	7.403	2.344	4.681	−0.426	−5.639	−0.561
	2	12.136	21.059	7.021	0.577	−3.281	0.463
	3	15.86	1.759	12.625	−0.346	1.776	−0.712
	4	40.606	11.641	19.044	−0.521	7.104	−0.675
	5	26.961	0.921	23.787	−0.008	10.939	−0.180
	6	17.859	3.845	28.856	0.064	14.972	0.024
	7	16.091	4.370	31.705	−0.348	17.587	0.491
	8	24.767	9.035	30.546	−0.480	16.763	0.099
	9	13.698	13.381	26.521	−0.116	12.458	0.040
	10	8.051	2.392	20.555	−0.179	6.545	−0.840
	11	2.483	1.166	14.23	0.153	1.046	−0.911
	12	4.475	10.181	6.932	0.500	−3.454	−0.202

6. Results and discussion

According to the concept of varying probability distributions, suitable probability distribution functions are examined for all the time scales of SPI, SPEI and SPTI indicators, independently. The suitability of probability function from the list of candidate distributions is decided under the minimum value of the Bayesian Information Criterion (BIC) criterion. That is, the CDF of those probability functions which have minimum BIC values is standardized. Table 3 provides BIC values against the candidate probability distributions for SPI, SPEI and SPTI indicators at a one-month time scale. We observed that 3 P Weibull distribution is the most suitable for SPI (BIC = −1036.51) and SPTI (BIC = −483.52) indicators at Astor station, while Trapezoidal distribution is the most appropriate fitted for SPEI indicator (BIC = − 710.05). Once again, the 3 P Weibull distribution has minimum values for SPI (BIC = 1030.98) and SPTI (BIC = −188.38) at Bunji station, whereas Johnson SB distribution is the most appropriate distribution for SPEI indicator (BIC = −1248.44). Accordingly, 4 P Beta distribution is the most appropriate for SPI (BIC = −788.69) and SPTI (BIC = −374.24) indicators, and Johnson SB distribution has the best candidacy score (BIC = −977.62) for SPEI indicator at Gupis station.

Table 3. BIC for Scale 1 SPI, SPEI and SPTI for Astore Bunji and Gupis.

	Astore			Bunji			Gupis
Distribution	SPI	SPEI	SPTI	SPI	SPEI	SPTI	SPI	SPEI	SPTI
2P Beta	−735.45	−565.00	−287.79	−722.44	−1019.68	−125.12	−399.44	−775.08	−254.43
3P Weibull	−1036.51	−700.54	−483.52	−1030.98	−1178.06	−188.38	−777.57	−910.74	−370.55
4P Beta	−1031.38	−700.28	−473.37	−1020.69	−1210.97	−181.74	−788.69	−823.15	−374.24
Arcsine	−853.92	−643.79	−274.88	−791.05	−1089.24	7.71	−577.40	−817.37	−151.87
Burr	−777.82	−656.97	−332.22	−753.84	−1152.29	158.50	−537.17	−905.52	−21.11
Cauchy	−906.22	−673.83	−357.62	−790.41	−1147.13	−2.03	−640.15	−884.69	−235.77
Chi	−772.18	−571.38	−253.75	−747.42	−1024.70	−89.65	−535.49	−779.77	−188.78
Chi-Square	−778.82	−570.07	−450.78	−771.56	−1024.59	−49.17	−549.29	−779.77	−146.95
Cosine	−827.71	−706.01	−182.80	−582.27	−1175.57	165.44	−409.73	−880.80	13.08
Curvilinear Trapezoidal	−770.35	−590.81	−369.28	−561.78	−1128.83	−40.20	−399.24	−809.67	−158.55
Exponential	−992.52	−569.99	−398.97	−727.50	−1024.70	30.67	−596.82	−779.77	−206.77
F-	−735.45	−565.75	−459.87	−863.89	−1019.68	2.54	−623.87	−775.23	−124.81
Gamma	−1024.97	−573.48	−482.01	−879.99	−1019.80	−116.91	−667.73	−775.23	−294.28
Generalized Extreme Value	−959.45	−700.89	−450.09	−886.80	−1216.78	−88.95	−744.99	−943.97	−338.99
Generalized normal	−974.59	−699.88	−472.53	−929.73	−1208.10	−124.81	−759.83	−938.39	−355.54
Gumbel	−942.14	−692.28	−382.34	−754.06	−1159.45	29.41	−607.56	−881.41	−205.28
Inverse Chi-Square	−847.29	−573.96	−260.44	−849.03	−1024.70	74.59	−615.75	−779.84	−126.69
Inverse Gamma	−894.39	−575.70	−393.13	−864.79	−1019.80	−75.36	−737.28	−775.30	−333.93
Inverse Gaussian	−848.87	−565.00	−355.22	−790.50	−1019.68	−12.65	−740.25	−775.08	−334.88
Johnson SB	−971.05	−640.08	−475.88	−925.95	−1248.44	−122.60	−457.52	−977.62	−351.09
Johnson SU	−969.88	−695.65	−467.89	−924.84	−1203.21	−120.21	−755.59	−936.90	−351.09
Laplace	−914.20	−683.10	−377.36	−756.94	−1150.62	10.91	−630.07	−886.00	−238.05
Logistic	−920.09	−697.50	−356.29	−752.10	−1166.09	37.57	−596.22	−895.15	−187.98
Log-normal	−970.31	−574.42	−474.91	−882.30	−1019.80	−105.32	−729.92	−775.30	−338.98
Normal	−917.65	−702.40	−348.30	−747.38	−1170.79	46.87	−585.58	−898.11	−172.80
Rayleigh	−934.02	−700.69	−365.81	−748.60	−1166.60	42.39	−591.38	−888.30	−182.47
Scaled/shifted t-	−915.53	−698.17	−359.42	−814.16	−1165.90	−11.77	−642.21	−893.57	−237.55
Skewed-normal	−951.45	−700.79	−378.12	−744.10	−1218.27	44.30	−590.96	−937.90	−185.69
Trapezoidal	−940.42	−710.05	−354.12	−712.78	−1224.00	80.19	−570.94	−946.94	−151.71
Triangular	−945.13	−685.96	−358.76	−664.59	−1228.23	76.39	−557.57	−947.69	−148.56
Uniform	−741.92	−664.32	−220.71	−622.63	−1084.39	121.33	−478.66	−884.52	−37.54
von Mises	−752.43	−568.01	−355.83	−752.74	−1019.80	46.82	−531.25	−775.31	−173.35

Figure 4 shows theoretical vs empirical histograms of all the selected probability distributions for SPI, SPEI and SPTI indicators at one-month time scales. We have observed that the histograms of SPTI indicators are more precise than SPI at Astor station. Comparative to SPI and SPTI, there is a significant deviation between empirical and theoretical plots of SPEI indicator.

Fig. 4. Probability distribution at scale one of SPI, SPEI and SPTI at Astore, Bunji and Gupis.

Analogous to one-month time scales, inferences are made to search appropriate probability distributions for other time scales of SPI, SPEI and SPI indices. Table 4 shows that the probability distribution functions which are observed throughout the screening process. Table 5 shows the BIC values associated with probability distributions for SPI, SPEI and SPTI. In all time scales, varying probability functions have varying levels of candidacy.

Table 4. Characteristics of some important distributions.

Sr. No	Distributions	Probability function	CDF	Range
1	3P Weibull	$P(x)=\frac{\alpha }{\beta }{\left(\frac{x-\gamma }{\beta }\right)}^{\alpha -1}\hspace{0.17em}\text{exp}\hspace{0.17em}\left[-{\left(\frac{x-\gamma }{\beta }\right)}^{\alpha }\right]$	F$(x)=1-\mathrm{}\hspace{0.17em}\text{exp}\hspace{0.17em}\left[-{\left(\frac{x-\gamma }{\beta }\right)}^{\alpha }\right]$	$\gamma \le x<+\text{∞}$
2	Gamma	$f(x)=\frac{{\left(x\right)}^{\alpha -1}}{{\beta }^{\alpha }\Gamma (\alpha )}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{-\left(x\right)}{\beta }\right)$	$F(x)=\frac{\Gamma (\alpha ,\frac{x}{\beta })}{\Gamma (\alpha )}$	$x>0$
3	Normal	$f(x)=\frac{1}{\sigma \sqrt{2\pi }}\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{-1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2\right]$	$F(x)=\frac{1}{2}\left[1+\text{erf}\hspace{0.33em}\left(\frac{x-\mu }{\sigma \sqrt{2}}\right)\right]$	$-\mathrm{\infty }<x<+\mathrm{\infty }$
4	Cosine	$f(x)=\frac{1}{2s}\left[1+\left.\text{cos}\left(\frac{x-\mu }{s}\pi \right)\right]\right.$	$F(x)=\frac{1}{2}\left[1+\left.\frac{x-\mu }{s}+\frac{1}{\pi }\text{sin}\left(\frac{x-\mu }{s}\pi \right)\right]\right.$	$\mathrm{\mu }-s\le x\le \mathrm{\mu }+s$
5	Skewed normal	$f(x)=\frac{2}{\omega \sqrt{2\pi }}\text{exp}\hspace{0.17em}\left(-\left(\frac{x-\mathrm{\xi }}{2{\omega }^2}\right)\right){\int }_{-\mathrm{\infty }}^{\alpha \left(\frac{x-\mathrm{\xi }}{\omega }\right)}\frac{1}{\sqrt{2\pi }}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\left(\frac{t^2}{2}\right)\right)dt$	F$(x)=\varphi \left(\frac{x-\mathrm{\xi }}{2\omega }\right)-2T\left(\frac{x-\mathrm{\xi }}{2\omega },\alpha \right)$ $T\left(h,\alpha \right)\mathrm{}is\mathrm{}\mathit{Owe}{n}^{\mathrm{\prime }}s\mathrm{}\mathit{Function}$	$-\mathrm{\infty }<x<+\mathrm{\infty }$
6	Trapezoidal	$f(x)=\left\{\begin{array}{l}\frac{2}{d+c-a-b}\frac{x-a}{b-a}\mathit{for}a\le x<b\\ \frac{2}{d+c-a-b}\mathit{for}b\le x<c\\ \frac{2}{d+c-a-b}\frac{d-x}{d-c}\mathit{for}c\le x\le d\end{array}\right\}$	$F\left(x\right)=\left\{\begin{array}{l}\frac{1}{d+c-a-b}\frac{1}{b-a}{\left(x-a\right)}^2\mathit{for}a\le x<b\\ \frac{1}{d+c-a-b}\left(2x-a-b\right)\mathit{for}b\le x<c\\ 1-\frac{1}{d+c-a-b}\frac{1}{d-c}{\left(d-x\right)}^2\mathit{for}c\le x\le d\end{array}\right\}$	$a\le x\le d$
7	Triangular	$f(x)=\left\{\begin{array}{l}\frac{2\left(x-a\right)}{\left(b-a\right)\left(c-a\right)}\mathit{for}a\le x<c\\ \frac{2}{b-a}\mathit{for}x=c\\ \frac{2\left(b-x\right)}{\left(b-a\right)\left(b-c\right)}\mathit{for}c<x\le b\end{array}\right\}$	F$(x)=\left(\begin{array}{l}\frac{(x-a)²}{\left(b-a\right)\left(c-a\right)}\mathit{for}a\le x<c\\ 1-\frac{(b-x)²}{\left(b-a\right)\left(b-c\right)}\mathit{for}x=c\\ 1\mathit{for}c<x\le b\end{array}\right)$	$a\le x\le b$
8	Johnson SU	$f\left(x\right)=\frac{\delta }{\lambda \sqrt{2\pi }\mathrm{}\sqrt{z^2+1}}\mathrm{e}\mathrm{xp}\left(-\frac{1}{2}{\left(\mathrm{\gamma }+\mathrm{\delta }\ln \left(z+\sqrt{z^2+1}\right)\right)}^2\right)$$z=\frac{x-\mathrm{\xi }}{\lambda }$	$F\left(x\right)=\varphi \left(\mathrm{\gamma }+\mathrm{\delta }\ln \left(z+\sqrt{z^2+1}\right)\right)$$\varphi \mathrm{}is\mathrm{}\mathit{Laplace}\mathrm{}\mathit{integral}$	$-\mathrm{\infty }<x<+\mathrm{\infty }$
9	Johnson SB	$f\left(x\right)=\frac{\delta }{\lambda \sqrt{2\pi }\mathrm{}z\left(1-z\right)}\mathrm{e}\mathrm{xp}\left(-\frac{1}{2}{\left(\mathrm{\gamma }+\mathrm{\delta }\ln \left(\frac{z}{1-z}\right)\right)}^2\right)$	$F\left(x\right)=\varphi \left(\mathrm{\gamma }+\mathrm{\delta }\ln \left(\frac{z}{1-z}\right)\right)$$\varphi \mathrm{}is\mathrm{}\mathit{Laplace}\mathrm{}\mathit{integral}$	$\mathrm{\xi }\le x\le \mathrm{\xi }+\mathrm{\lambda }$
10	Log normal	$f(x)=\frac{\hspace{0.17em}\text{exp}\hspace{0.17em}\left[\frac{-1}{2}{\left(\frac{In(x)-\mu }{\sigma }\right)}^2\right]}{\left(x\right)\sigma \sqrt{2\pi }}$	$f\left(x\right)=\frac{1}{2}+\frac{\mathit{erf}\left[\frac{In(x)-\mu }{\sqrt{2}\sigma }\right]}{2}$	$0<x<+\mathrm{\infty }$
11	Laplace	$f\left(x\right)=\frac{1}{2b}+\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{\left|x-\mu \right|}{b}\right)$	$F\left(x\right)=\left\{\begin{array}{l}\frac{1}{2}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{x-u}{b}\right)\mathrm{}\mathit{for}\mathrm{}x\le u\mathrm{}\\ 1-\frac{1}{2}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x-u}{b}\right)\mathit{for}x\ge u\end{array}\right\}$	$x\in R$
12	Logistic	$f\left(x\right)=\frac{\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x-\mu }{s}\right)}{s{\left(1+\hspace{0.17em}\text{exp}\hspace{0.17em}\left(\frac{x-\mu }{s}\right)\right)}^2}$	$F(x)=\frac{1}{1+\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x-\mu }{s}\right)}$	$-\mathrm{\infty }<x<+\mathrm{\infty }$
13	Gumbel	$f\left(x\right)=\frac{1}{\beta }\hspace{0.17em}\text{exp}\hspace{0.17em}\left[-(z+\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-z\right))\right]$$z=\left(\frac{x-\mu }{\beta }\right)$	$\text{exp}\hspace{0.17em}\left(-\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x-\mu }{\beta }\right)\right)$	$x\in R$
14	Rayleigh	$f\left(x\right)=\frac{x}{{\sigma }^2}\mathrm{e}\mathrm{xp}\left(\frac{{-x}^2}{{(2\sigma }^2)}\right)$	$F\left(x\right)=1-\exp \left(\frac{{-x}^2}{\left({2\sigma }^2\right)}\right)$	$0<x<+\mathrm{\infty }$
15	4P Beta	$f(x)=\mathrm{}\frac{{(x-\mathit{min})}^{\alpha -1}{(\mathit{max}-x)}^{\beta -1}}{\mathit{B}\left(\alpha ,\mathrm{}\beta \right){(\mathit{max}-\mathit{min})}^{\alpha +\beta -1}}$	$F(x)=\mathrm{}\frac{\mathit{B}\left(\frac{x}{\alpha },\mathrm{}\beta \right)}{\mathit{B}\left(\alpha ,\mathrm{}\beta \right)}$	$\mathit{min}<x<\mathit{max}$

Table 5. BIC for all scale of SPI, SPEI and SPTI indicators at Astore, Bunji and Gupis stations.

Astor				Bunji		Gupis
Index	Scales	Distribution	BIC	Distribution	BIC	Distribution	BIC
For SPI	1	3P Weibull	−1036.51	3P Weibull	−1030.98	4P Beta	−788.69
	3	Gamma	−1279.05	Gamma	−824.87	Gamma	−1264.95
	6	Gamma	−892.67	Skewed-normal	−1162.15	Gumbel	−1305.35
	9	Gamma	−896.09	Normal	−649.05	Johnson SU	−1518.96
	12	Cosine	−913.33	Laplace	−688.10	Johnson SU	−937.61
	24	Skewed-normal	−1294.42	Laplace	−843.74	Scaled/Shifted t	−1407.98
For SPEI	1	Trapezoidal	−710.05	Johnson SB	−1248.44	Johnson SB	−977.62
	3	Trapezoidal	−941.49	Johnson SB	−1323.78	Johnson SB	−1098.61
	6	Trapezoidal	−1440.21	Johnson SB	−1094.44	Trapezoidal	−1482.59
	9	Triangular	−1405.81	Trapezoidal	−976.73	Trapezoidal	−1513.12
	12	Trapezoidal	−1052.48	Logistic	−1158.15	Laplace	−1063.66
	24	Johnson SU	−1471.48	Gumbel	−1508.11	Laplace	−1474.12
For SPTI	1	3P Weilbull	−483.52	3P Weibull	−188.38	4P Beta	−374.24
	3	Johnson SU	−542.87	Gumbel	−190.69	Chi-Square	−432.18
	6	Log-normal	−721.00	Skewed-normal	−300.34	Johnson SU	−410.04
	9	Gamma	−725.22	Rayleigh	−380.83	Johnson SU	−463.21
	12	Triangular	−702.87	Laplace	−411.86	Johnson SU	−466.57
	24	Laplace	−564.00	Trapezoidal	−304.73	Scaled/Shifted t	−629.46

From the empirical analysis, we have observed that the deviations between theoretical and empirical histograms (see figure 4) are due to the ineffectiveness of the probability distribution functions for extreme values. Consequently, improper and subjective choice of probability distribution functions has the potential to produce imprecise drought indices values. Hence, it’s reasonable to say that the choice of probability functions using BIC values is purely subjective in nature. That is, there is a big gap in the search for more appropriate probability functions.

Our experimental results agree with the conclusion of past research (see Balint et al., 2013; Bayissa et al., 2018; Zhu et al., 2018). Therefore, by the rationale of our criterion, the use of long-term behavior of drought classes can overcome the effect of extreme values in reporting accurate drought class. In recent years, some authors are working with mixture probability distribution functions (Mallya et al., 2015), while some authors have recommended non-parametric functions-based standardization (Farahmand and AghaKouchak, 2015). Although these approaches can be used to reduce the risk of inaccurate drought classes, yet the solution of discrepancies among drought indices has not addressed.

After the estimation of drought indices, their temporal quantitative values are categorized by their severity level. Resultant qualitative vectors are then used to construct transition probability matrices using Markovchain (Spedicato et al., 2016) R package, independently and separately with respect to time scales and stations. Figure 5 shows the graphical display of the transition probability matrix of SPTI index at one-month time scale. We noted that, except SPTI-1 of Astore station, there is no single transition between Severe Wet to Severe Dry drought categories. The results related to other time scales are archived in the author’s gallery.

Fig. 5. Graphical representation of TPM of SPI, SPEI and SPTI at scale 1 in Astore, Bunji and Gupis.

After the computation of transition probability matrices, the long-term behavior of each drought category is quantified using steady-state probabilities. Table 6 shows the steady-state probability of SPI, SPEI, and SPTI drought indices at Astor, Gupis, and Bunji. Here, NA values show the absence of drought categories in the temporal vectors of drought classification state.

Table 6. Steady-states probabilities of various drought classes.

			ED	EW	MD	MW	ND	SD	SW	Sum
Astore		SPI	NA	0.02309	0.09254	0.09583	0.67654	0.06949	0.04252	1.0
	Scale-1	SPEI	NA	0.03013	0.16716	0.08661	0.66113	0.02135	0.03361	1.0
		SPTI	NA	0.01952	0.07995	0.11872	0.67887	0.05339	0.04954	1.0
		SPI	NA	0.02674	0.09626	0.07487	0.76292	NA	0.03922	1.0
	Scale-3	SPEI	NA	0.01756	0.16475	0.11997	0.63395	0.03761	0.02616	1.0
		SPTI	NA	0.01575	0.13215	0.10907	0.65843	0.01267	0.07193	1.0
		SPI	0.01905	0.02137	0.10311	0.10330	0.66614	0.04428	0.04275	1.0
	Scale-6	SPEI	NA	0.02151	0.15233	0.11470	0.63620	0.03226	0.04301	1.0
		SPTI	NA	0.01603	0.13401	0.09798	0.65198	0.04300	0.05700	1.0
		SPI	NA	0.07027	0.02162	0.08468	0.77117	NA	0.05225	1.0
	Scale-9	SPEI	0.04573	0.00178	0.08788	0.15350	0.62473	0.05424	0.03213	1.0
		SPTI	0.00727	0.01791	0.10650	0.08598	0.67892	0.04968	0.05374	1.0
		SPI	0.01131	0.01974	0.10654	0.11128	0.66767	0.05654	0.02692	1.0
	Scale-12	SPEI	0.02434	0.02316	0.09480	0.09800	0.68597	0.05771	0.01604	1.0
		SPTI	NA	0.02136	0.13068	0.09610	0.64246	0.06491	0.04449	1.0
		SPI	0.01111	0.03704	0.07037	0.05741	0.73889	0.07593	0.00926	1.0
	Scale-24	SPEI	0.01667	NA	0.12778	0.10556	0.66111	0.06481	0.02407	1.0
		SPTI	NA	NA	0.12407	0.11111	0.64074	0.07407	0.05000	1.0
		SPI	NA	0.02131	0.28419	0.07460	0.57549	NA	0.04441	1.0
	Scale-1	SPEI	NA	0.00534	0.16201	0.08852	0.68365	0.04629	0.01419	1.0
		SPTI	NA	0.01243	0.28419	0.09414	0.55773	NA	0.05151	1.0
		SPI	0.01248	0.01070	0.05526	0.07665	0.77897	0.01604	0.04991	1.0
	Scale-3	SPEI	NA	0.00713	0.16934	0.14617	0.63815	0.03209	0.00713	1.0
		SPTI	NA	0.02674	0.09626	0.07487	0.76292	NA	0.03922	1.0
		SPI	NA	0.03943	0.09140	0.05376	0.77240	NA	0.04301	1.0
	Scale-6	SPEI	NA	0.00713	0.16934	0.14617	0.63815	0.03209	0.00713	1.0
		SPTI	NA	0.02674	0.09626	0.07487	0.76292	NA	0.03922	1.0
Gupis		SPI	NA	0.06126	0.04144	0.08288	0.74955	NA	0.06486	1.0
	Scale-9	SPEI	0.01622	0.01622	0.11892	0.08649	0.68829	0.04505	0.02883	1.0
		SPTI	NA	0.07027	0.02162	0.08468	0.77117	NA	0.05225	1.0
		SPI	NA	0.06703	0.04891	0.07246	0.74457	NA	0.06703	1.0
	Scale-12	SPEI	NA	0.03623	0.10326	0.10507	0.69565	0.00543	0.05435	1.0
		SPTI	0.03986	0.00362	0.07428	0.10870	0.69746	0.01812	0.05797	1.0
		SPI	0.03704	NA	0.05741	0.13148	0.65926	0.04444	0.07037	1.0
	Scale-24	SPEI	NA	0.02963	0.10926	0.11111	0.67037	NA	0.07963	1.0
		SPTI	NA	0.04074	0.05741	0.08333	0.73704	NA	0.08148	1.0
		SPI	NA	0.01422	0.20121	0.08695	0.64801	NA	0.04961	1.0
	Scale-1	SPEI	NA	0.00706	0.18527	0.17128	0.59370	0.02673	0.01597	1.0
		SPTI	NA	0.01424	0.20472	0.07275	0.66037	NA	0.04793	1.0
		SPI	NA	0.02674	0.09626	0.07487	0.76292	NA	0.03922	1.0
	Scale-3	SPEI	NA	0.01756	0.16475	0.11997	0.63395	0.03761	0.02616	1.0
		SPTI	NA	0.01575	0.13215	0.10907	0.65843	0.01267	0.07193	1.0
		SPI	0.01905	0.02137	0.10311	0.10330	0.66614	0.04428	0.04275	1.0
	Scale-6	SPEI	NA	0.02151	0.15233	0.11470	0.63620	0.03226	0.04301	1.0
		SPTI	NA	0.01603	0.13401	0.09798	0.65198	0.04300	0.05700	1.0
Bunji		SPI	NA	0.07027	0.02162	0.08468	0.77117	NA	0.05225	1.0
	Scale-9	SPEI	0.04573	0.00178	0.08788	0.15350	0.62473	0.05424	0.03213	1.0
		SPTI	0.00727	0.01791	0.10650	0.08598	0.67892	0.04968	0.05374	1.0
		SPI	0.01131	0.01974	0.10654	0.11128	0.66767	0.05654	0.02692	1.0
	Scale-12	SPEI	0.02434	0.02316	0.09480	0.09800	0.68597	0.05771	0.01604	1.0
		SPTI	NA	0.02136	0.13068	0.09610	0.64246	0.06491	0.04449	1.0
		SPI	NA	NA	0.14815	0.12778	0.61111	0.04815	0.06481	1.0
	Scale-24	SPEI	0.05185	NA	0.06111	0.18704	0.64630	0.04815	0.00556	0.9
		SPTI	NA	NA	0.14444	0.13519	0.61111	0.02778	0.08148	1.0

To combine three vectors of drought indices (i.e. SPI, SPEI, and SPTI), steady-states probabilities of each drought category are employed as weights. Those drought classes that have maximum values of corresponding steady-states probabilities are extracted in a separate vector. The resultant vectors are denoted as the LAWJADC index. Accordingly, the PWJADI criterion is applied for all the time scales at all stations. Quantitative time series values of all the involved drought indices, LAWJADC, and PWJADI are presented in Figures 6–8.

Fig. 6. Temporal plot of SPI, SPEI, SPTI, PWJADI and LAWJADC for Astore.

Fig. 7. Temporal plot of SPI, SPEI, SPTI, PWJADI and LAWJADC for Bunji.

Fig. 8. Temporal plot of SPI, SPEI, SPTI, PWJADI and LAWJADC for Gupis.

In Astore station, the behavior of LAWJADC is very close to PWJADI at one, six, nine and twelve and 24-month time scale. However, a little deviation between PWJADI and LAWJADC is observed at the three-month time scale (see Figure 8). Whereas, some discrepancies are found in the 6-month and 24-month time scale at Bunji station (see Figure 9). However, homogenous temporal behavior is observed in the remaining time scales. Accordingly, there exists a little deviation in three-month time scales at Gupis station (see Figure 10).

Fig. 9. Bar plot counts of drought categories observed in SPI, SPEI, SPTI, PWJADI and LAWJADC for Astore.

Fig. 10. Bar plot counts of drought categories observed in SPI, SPEI, SPTI, PWJADI and LAWJADC for Bunji.

Insight deviations among the proportion of each drought, categories are explored using the bar plot. It is observed that the proportions of the Normal Drought category between PWJADI and LAWJADC are looks equally distributed in all the scales. However, the distribution of frequencies varies in Severe and Extreme drought case. Irrespective of PWJADI, LAWJADC has fetched some frequencies of Extreme Drought and Extreme Wet drought categories in nearly all the time scales (see Figures 9–11). These deviations advocate the use of LAWJADC where the assumptions of PWJADI seem to fail.

Fig. 11. Bar plot counts of drought categories observed in SPI, SPEI, SPTI, PWJADI and LAWJADC for Gupis.

To check the association among each drought index, Kappa measure is used to assess the correlation of the proposed method and PWJADI with all-time scales. Table 7 shows the Kappa values of PWJADI and LAWJADC with SPI, SPEI, and SPTI in all the time scales in all the three stations. In most of the data sets, LAWJADC has weak correlation between SPEI. This is due to the inappropriate fitness of the probability function. These findings are consistent with previous research (see Ali et al., 2019a). However, there are strong correlations among PWJADI, SPI, and SPTI. Hence, in case of violation of the assumption of Markov’s condition, the proposed method is a strong candidate for the aggregation of multiple indices.

Table 7. Kappa values showing association among various drought indices and criterion.

			Kappa					P-value
Station	Scale	Method	SPI	SPEI	SPTI	PWJADI	LAWJADC	SPI	SPEI	SPTI	PWJADI	LAWJADC
	Scale-1	LAWJADC	0.2850	0.0176	0.2220	0.8250	1.0000	0.0000	0.6030	0.0000	0.0000	0.0000
Astore		PWJADI	0.2970	0.0970	0.2430	1.0000	0.8250	0.0000	0.0048	0.0000	0.0000	0.0000
	Scale-3	LAWJADC	0.2650	0.1060	0.1090	0.4800	1.0000	0.0000	0.0002	0.0002	0.0000	0.0000
		PWJADI	0.4840	0.0381	0.1630	1.0000	0.4800	0.0000	0.2920	0.0000	0.0000	0.0000
	Scale-6	LAWJADC	0.3900	0.1730	0.2790	0.5560	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
		PWJADI	0.3210	0.4060	0.3870	1.0000	0.5560	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-9	LAWJADC	0.2180	0.0263	0.0615	0.6320	1.0000	0.0000	0.2580	0.0216	0.0000	0.0000
		PWJADI	0.4320	0.0189	0.1590	1.0000	0.6320	0.0000	0.5570	0.0000	0.0000	0.0000
	Scale-12	LAWJADC	0.1850	0.0029	0.4960	0.3920	1.0000	0.0000	0.9400	0.0000	0.0000	0.0000
		PWJADI	0.2960	0.3630	0.5310	1.0000	0.3920	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-24	LAWJADC	0.1730	0.3250	0.3140	0.0683	1.0000	0.0000	0.0000	0.0000	0.0880	0.0000
		PWJADI	0.4500	0.0638	0.0515	1.0000	−0.0683	0.0000	0.0717	0.2100	0.0000	0.0880
	Scale-1	LAWJADC	0.3430	0.1190	0.2850	0.9040	1.0000	0.0000	0.0002	0.0000	0.0000	0.0000
Bunji		PWJADI	0.2980	0.1710	0.2990	1.0000	0.9040	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-3	LAWJADC	0.1160	0.0689	0.2330	0.6160	1.0000	0.0000	0.0047	0.0000	0.0000	0.0000
		PWJADI	0.1530	0.1530	0.3640	1.0000	0.6160	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-6	LAWJADC	0.0727	0.1300	0.1840	0.1730	1.0000	0.0064	0.0000	0.0000	0.0000	0.0000
		PWJADI	0.5990	0.0363	0.6270	1.0000	0.1730	0.0000	0.3820	0.0000	0.0000	0.0000
	Scale-9	LAWJADC	0.2240	0.2060	0.2670	0.5780	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
		PWJADI	0.2930	0.2540	0.2950	1.0000	0.5780	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-12	LAWJADC	0.4890	0.5080	0.5670	0.6750	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
		PWJADI	0.5380	0.5680	0.7560	1.0000	0.6750	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-24	LAWJADC	0.2200	0.6010	0.2770	0.7300	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
		PWJADI	0.3920	0.5520	0.4070	1.0000	0.7300	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-1	LAWJADC	0.2060	0.1280	0.2250	0.9120	1.0000	0.0000	0.0004	0.0000	0.0000	0.0000
Gupis		PWJADI	0.2120	0.1830	0.2450	1.0000	0.9120	0.0000	0.0000	0.0000	0.0000	0.0000
	Scale-3	LAWJADC	0.1420	−0.0109	0.2040	0.6590	1.0000	0.0000	0.6680	0.0000	0.0000	0.0000
		PWJADI	0.2000	0.0145	0.4010	1.0000	0.6590	0.0000	0.6530	0.0000	0.0000	0.0000
	Scale-6	LAWJADC	0.3340	0.0581	0.2880	0.6560	1.0000	0.0000	0.0478	0.0000	0.0000	0.0000
		PWJADI	0.5500	0.0739	0.0739	1.0000	0.6560	0.0000	0.0461	0.0461	0.0000	0.0000
	Scale-9	LAWJADC	0.4020	0.1390	0.4960	0.5850	1.0000	0.0000	0.0001	0.0000	0.0000	0.0000
		PWJADI	0.6030	0.1490	0.7340	1.0000	0.5850	0.0000	0.0002	0.0000	0.0000	0.0000
	Scale-12	LAWJADC	0.7520	0.3060	0.0715	0.8780	1.0000	0.0000	0.0000	0.0782	0.0000	0.0000
		PWJADI	0.7100	0.3790	0.0264	1.0000	0.8780	0.0000	0.0000	0.5200	0.0000	0.0000
	Scale-24	LAWJADC	0.5860	0.1750	0.0809	0.8770	1.0000	0.0000	0.0000	0.0400	0.0000	0.0000
	Scale-24	PWJADI	0.5430	0.1930	0.0658	1.0000	0.8770	0.0000	0.0000	0.1140	0.0000	0.0000

7. Conclusion

Continuous drought monitoring is an emerging feature of hydrological research. Drought indices are the most commonly practiced methods for drought monitoring in many drought monitoring modules. However, the existence of numerous drought indices creates a chaotic situation for users, analyzers, and databases of hydrological statistics. This research proposed a new criterion-The Long Average Weighted Joint Aggregative Drought Criterion (LAWJADC) - for combining two or more drought indices. We found that LAWJADC has sufficient features for its deployment in the combinative criterion list of drought indices.

Acknowledgement

The authors are very grateful to the deanship of scientific research at King Khalid University, Abha, Saudi Arabia, for the financial support through the Research Program under project number R.G.P 2/67/41.

Disclosure statement

No potential conflict of interest was reported by the authors.