Monthly drought prediction based on ensemble models

Drought is a natural hazard, which is a result of a prolonged shortage of precipitation, high temperature and change in the weather pattern. Drought harms society, the economy and the natural environment, but it is difficult to identify and characterize. Many areas of Pakistan have suffered severe droughts during the last three decades due to changes in the weather pattern. A drought analysis with the incorporation of climate information has not yet been undertaken in this study region. Here, we propose an ensemble approach for monthly drought prediction and to define and examine wet/dry events. Initially, the drought events were identified by the short term Standardized Precipitation Index (SPI-3). Drought is predicted based on three ensemble models i.e., Equal Ensemble Drought Prediction (EEDP), Weighted Ensemble Drought Prediction (WEDP) and the Conditional Ensemble Drought Prediction (CEDP) model. Besides, two weighting procedures are used for distributing weights in the WEDP model, such as Traditional Weighting (TW) and the Weighted Bootstrap Resampling (WBR) procedure. Four copula families (i.e., Frank, Clayton, Gumbel and Joe) are used to explain the dependency relation between climate indices and precipitation in the CEDP model. Among all four copula families, the Joe copula has been found suitable for most of the times. The CEDP model provides better results in terms of accuracy and uncertainty as compared to other ensemble models for all meteorological stations. The performance of the CEDP model indicates that the climate indices are correlated with a weather pattern of four meteorological stations. Moreover, the percentage occurrence of extreme drought events that have appeared in the Multan, Bahawalpur, Barkhan and Khanpur are 1.44%, 0.57%, 2.59% and 1.71%, respectively, whereas the percentage occurrence of extremely wet events are 2.3%, 1.72%, 0.86% and 2.86%, respectively. The understanding of drought pattern by including climate information can contribute to the knowledge of future agriculture and water resource management.


INTRODUCTION
In the last few decades, global warming has become an undeniable fact which causes climate irregularities i.e., extreme weather events and droughts (Bai et al., 2020;Horton et al., 2016;Mishra & Singh, 2010). Drought is the most complicated hazard among other natural disasters (e.g., flood, earthquake, tropical cyclones) and has a huge impact on society, the economy and natural environments, such as reducing air and water quality, and causing erosion, landscaping (dust), and ecological habitat damage (Bonsai & Wheaton, 2005). It is challenging to determine the characteristic of droughts such as beginning, ending, intensity and duration (Vicente-Serrano, Beguería & López-Moreno, 2010). Sharafati & Pezeshki (2020) proposed a robust approach to assess climate change impact variability on future extreme events over Dehbar catchment in Iran. They used the Lars-WG6 model to generate future climate weather variables such as temperature and rainfall depth using five coupled models with several emission scenarios. They concluded that trends in extreme rainfall depth and river discharge is increasing due to climate and inferred that the future extreme rainfall depth at more than 500-year return periods has more variability.
Several drought indices are used for the detection and characterization of drought, i.e., Standardized Precipitation Index (SPI) (Mckee, Doesken & Kleist, 1993), Standardized Anomaly Index (SAI) (Katz & Glantz, 1986), Standardized Precipitation Evapotranspiration Index (SPEI), (Vicente-Serrano, Beguería & López-Moreno, 2010), Standardized Precipitation Temperature Index (SPTI) (Zulfiqar et al., 2019;Zuliqar et al., 2017) and Reconnaissance Drought Index (RDI) (Tsakiris & Vangelis, 2005). SPI is a popular drought index and is used to identify and monitor drought in the past studies (Achour et al., 2020;Bai et al., 2020;Caloiero & Veltri, 2019;Ellahi et al., 2020). During the estimation phase of SPI, the selection of unsuitable probability distribution may provide biased values of SPI (Zulfiqar et al., 2019;Stagge et al., 2015). Zulfiqar et al. (2019) fitted several appropriate probability distributions in the estimation phase of SPI, SPEI and SPTI. The goodness of fit test was used with 5% level of significance to check the appropriateness of candidate probability distribution and the lowest value of Bayesian Information Criterion (BIC) was used for the selection of suitable probability distribution. Mishra & Desai (2005) used the Autoregressive Integrated Moving Average (ARIMA) and the Seasonal Autoregressive Integrated Moving Average (SARIMA) model on the SPI series  in Kansabati River of India to predict the drought for 2-3 month time scale, and also stated that predicted results are consistent with the observed data. Belayneh, Adamowski & Khalil (2016) used the SPI-3 and SPI-6 of Awash River basin in Ethiopia to forecast drought through Artificial Neural Network (ANN), Support Vector Machine (SVM), and Wavelet Neural Network (WNN) model. The ARIMA model was applied to SPI-12 and SPI-24 to forecast drought and the results of the forecast were compared with machine learning models (Belayneh, Adamowski & Khalil, 2016). Among these several models, the WNN performed better based on R square, MAE and RMSE. Choubin, Malekian & Golshan (2016) utilized the rainfall data from January 1967 to December 2009 of four stations to compute the SPI (1-12) months in Iran. They applied the Adaptive Neuro Fuzzy Inference System (ANFIS), M5P model tree and the Multilayer Perceptron (MLP) model to predict the SPI with the incorporation of climate signal. The accuracy results indicated that the performance of MLP model is better in comparison to other models.
Moradkhani & Meier (2010) used two procedures e.g., Linear Regression (LR) and Ensemble Streamflow Prediction (ESP), for the ensemble prediction. The use of a statistical approach with large scale climate variable is effective for the hydrological prediction (Moradkhani & Meier, 2010;Najafi, Moradkhani & Piechota, 2012). Choubin et al. (2019) applied the Multivariate Adaptive Regression Splines (MARS), M5 Model Tree and the Least Squares Support Vector Machine (LSSVM) to predict the streamflow pattern over the Mediterranean region of Turkey. They concluded that the performance of the LSSVM model is better in comparison to other models that utilized the climate information for streamflow modelling. Also, it is stated that the North Pacific (NP) and the East Central Tropical Pacific Sea Surface Temperature (ECTP-SST) damage the streamflow patterns. Choubin et al. (2014) used the rainfall data from January 1967 to December 2009 of four stations in Iran. They used eight most relevant climate indices in the prediction of drought. They observed that the Atlantic surface temperature had the inverse relationship with SPI and Atlantic Meridional Mode (AMM) had the highest correlation. It was concluded that the forecast performance of the Neuro-Fuzzy (NF) model is better as compared to the Stepwise Regression (SR) model at the two-year lag. AghaKouchak (2015) proposed a Multivariate Standardized Drought Index (MSDI) for the ensemble drought prediction in Africa. In their study, firstly the ESP procedure was applied on the monthly rainfall series and soil moisture series to predict the seasonal changes. The estimation of MSDI is based on the joint probability of accumulated predicted seasonal changes of rainfall and soil moisture from ESP procedure. The ensemble model estimates the amount of severe drought and provide information about their probability of occurrence and results were found to be consistent with the observed data. Bradley, Habib & Schwartz (2015) applied a bootstrap sampling resampling procedure on the daily flow time series (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) of the Blue Nile river of Ethiopia to produce the ensemble stream flow forecast. The results of the technique were compared with the traditional weighting approach. Moreover, the traditional procedure utilized the Gaussian kernel function for the weighting scheme, whereas the sampling resampling approach uses the concept of Bayesian updating (Smith & Gelfand, 1992). Both techniques assigned weights to each ensemble member in a dissimilar way, but the results of Bayesian updating was found to be better. Zhang, Dong & Chen (2019) proposed three ensemble models, i.e., Equal Ensemble Drought Prediction (EEDP), Weighted Ensemble Drought Prediction (WEDP) and Conditional Ensemble Drought Prediction (CEDP) to forecast the monthly droughts of 26 meteorological stations in Jiangxi province of China. As a result, the CEDP model was found to be better in parameter estimation and accuracy than other models. The choice of climate variables is based on the strong correlation between the forecast variable and climate variable (Najafi, Moradkhani & Piechota, 2012;Zhang, Dong & Chen, 2019;Sharafati, Nabaei & Shahid, 2020) and further climate variables were combined into one Integrated Climate Predictor (ICP) through linear regression model. The hydrological variables do not fulfill the requirement of independence and normality. Zhang & Singh (2006) fitted the joint frequency distribution of flood peak, volume and duration using the data of Amite River Basin, Louisiana, USA. Several techniques were applied to obtain the conditional returned periods such as Archimedean copula, Gumbel mix, and bivariate normal Box-Cox transformation. The copula distribution was found to be better in frequency estimation than others. The choice of suitable copula family is problematic, especially the hydrological variables provide imprecise information about the dependence structure (Renard & Lang, 2007). Therefore, the appropriate bivariate copula was fitted between the climate index and the rainfall series in Jiangxi province of China. The choice of appropriate copula family i.e., Frank, Clayton, Gumbel and AMH copula were based on the minimum value of BIC. Among these copula families, the performance of the Frank copula was better in comparison to others in the copula family (Zhang, Dong & Chen, 2019). Nabaei et al. (2019) have made a comprehensive spatial drought analysis to observe the joint return periods of droughts by using Copula. They applied Archimedean Copulas on various combinations of three characteristics of meteorological droughts (Severity (S), Peak (P) and Duration (D)) and concluded that the Gumbel Copula was the most appropriate for S-P according to the Copula Information Criterion, while the Clayton was most appropriate for S-D and P-D.
Changes in the weather pattern have become a serious problem for the biodiversity, hydrology, water resources, agriculture, forestry, human and livestock health (Farooqi, Khan & Mir, 2005). A drought analysis with the incorporation of climate information has not yet been undertaken in this study region. Therefore, the EEDP, WEDP and CEDP models have been used for describing the effect of climate indices on the drought condition. This study aims to describe and classify the short term wet and dry events. Further, these events will be examined, to find out the most occurring event.

Study area
The study area consists of four meteorological stations of Pakistan (Multan, Bahawalpur, Barkhan and Khanpur). These selected stations are located in the southern part of Pakistan, mostly very hot and mildly cold areas. According to the Meteorological Department of Pakistan, climate change in recent years has led to severe droughts in Southern Pakistan. Agricultural sector plays a leading role in Pakistan's economy. These stations are popular based on the agriculture sector. So, these stations have been selected based on the agriculture sector and also located in the Southern region (see Fig. 1). For the statistical analysis, the monthly quantitative data of temperature and precipitation from January 1990 to December 2018 has been obtained from the Karachi Data Processing Center through the Pakistan Meteorological Department, Karachi. Initially, the monthly average precipitation and temperature are computed for all stations which are depicted in Fig. 2.
According to (PMD, 2019), many climate indices are the reason for climate change in Southern Pakistan. So, the monthly data of eleven climate indices (Sea Surface

Standardized precipitation index
McKee, Nolan & Kleist (1993) proposed a Standardized Precipitation Index (SPI) for defining and monitoring of wet and dry events i.e., beginning, ending and intensity. The SPI is used to measure the precipitation shortage from the long-term historical record of precipitation and represents the quantitative definition of droughts on multiple time scales i.e., 3, 6, 9, 12, 24 and 48 months. Among these time scales, 3-months SPI is used for the short term, 6-months SPI is used for medium-term and 48-months SPI is used for long term drought analysis (Wu et al., 2007). In this paper, the cumulative precipitation series such as 3 months' time scale is used to calculate SPI. This cumulative precipitation time series is utilized to examine the appropriate probability distribution. In the next stage, the appropriateness of candidate probability distributions is checked by Kolmogorov Simonov (KS) (Justel, Peña & Zamar, 1997) and Anderson Darling (AD) (Anderson & Darling, 1952) test. We used the Easy fit (Schittkowski, 2002) software to compare the goodness of fit measure of several probability distributions. The selection of suitable probability distribution is identified based on the lowest value of BIC. After this, the parameters of each fitted probability distribution are utilized to compute the Cumulative Distribution Function (CDF). In various cases, precipitation time series consist of zero values, due to this reason CDF could be undefined at y = 0. For example, we fit the gamma distribution over a precipitation time series and the CDF of gamma distribution could be undefined at y = 0. According to Zulfiqar et al. (2019) andMishra &Desai (2005), an equation is used for this purpose that is given below: where, q represents the probability of zero precipitation in time series, while F (y) is the CDF of a selected probability distribution.
where, m is the number of zeros and n represents the total number of observations in the precipitation time series. Following Mckee, Doesken & Kleist (1993) and Mishra & Desai (2005), applying the procedure of approximate conversion developed by Abramowitz & Stegun (1948) that converts the CDF into a standard normal variable Z with mean zero and unit variance. The SPI index is as follows: Then, it will be used, where, c 0 = 2.515517 c 1 = 0.802853 c 2 = 0.010308 d 1 = 1.432788 d 2 = 0.189269 d 3 = 0.001308. In both cases, these multipliers are used. According to Edwards & McKee (1997) and Mckee, Doesken & Kleist (1993), the classification of SPI is presented in Table 2.
The positive values of SPI are indicating wet conditions with greater than median precipitation. Whereas, the negative values of SPI are indicating dry conditions with lower than median precipitation. The Standardized Precipitation Index (SPI) is a versatile tool in drought monitoring and analysis. The climatologists around the world had extensively used it for monitoring droughts. SPI is a simple index and only one indicator precipitation is required for the computation. SPI is also used in operational monitoring systems in various countries around the world and is recognized by the World Meteorological Organization (WMO). So, the authors mostly consider SPI based on these features.

Ensemble streamflow prediction model
ESP model uses the three types of series such as temperature, precipitation and SPI. The aim of using the ESP model is to assess monthly dry conditions. The inclusion of dependence structure (between precipitation and SPI) in the ESP model can strengthen the statistical drought analysis and are useful for estimating inter-annual variability. The choice of the precipitation series is based on the strong relationship between SPI and precipitation. Due to strong correlation, which potentially allows authentic monthly prediction of precipitation that improves monthly SPI prediction. The detailed methodology of initial mean prediction and ensemble mean prediction is described in Liu & Hwang (2015). Moreover, the statistical analysis has been done by the R package 'locfit ' (Loader, 2011). In the ESP model, the prediction of temperature, precipitation and SPI is required. Hydrological time series include various components such as periodically, serial correlation and this feature can be described only by ARIMA and SARIMA models Mishra & Desai (2005). (1977) provide a new forecasting tool known as ARIMA model. It is a modification of the ARMA model, which deals with non-stationary series. According to Mishra & Desai (2005), the general ARIMA model is as follows,

Non-seasonal model (ARIMA). Geurts
where, ϕ(β) is polynomial of Autoregressive (AR) model of order p, θ (β) is polynomial of the moving average model of order q and ∇ d is d th differencing operator.
Seasonal model (SARIMA). The SARIMA model is a generalization of the ARIMA model. The major advantage of using the SARIMA model is that it deals with the non-stationary series as well as the seasonal series. According to Durdu (2010) and Mishra & Desai (2005), the general SARIMA model is represented as, where p and q are the order of non-seasonal part of AR and MA and drepresents the non-seasonal differencing parameter. Besides, P and Q are the order of seasonal part of the AR and MA model, D is a seasonal differencing parameter and srepresents the length of the season. The development of time series model includes three stages e.g., identification, estimation of parameters and diagnostic check. Initially, the normality assumption of the time series is checked. In case of non-normal series, Box-Cox transformation and log transformation are applied to satisfy the normality condition. After this, Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) were used to scrutinize the structure of transformed series and this information is helpful for fitting the appropriate model. So, the initial parameters have been estimated visually from the plot of ACF and PACF of the time series. Then, the parameters of the model have been done by using Maximum Likelihood Estimation (MLE). Various combinations of parameters were used for the fitting of an appropriate model. Among these appropriate models, the choice of a suitable model was based on the minimum value of AIC. Further, the normality condition of residuals was checked by histogram and normal probability plot. We use the R package 'forecast' for the time series modelling.

Equal ensemble drought prediction model
The concept of EEDP model is based on the ESP model. In the EEDP model, the historical series that is assumed to have equal weights in the future. The climate series contains the n years for a study area so that the number of ensemble members are corresponding to the historical years in the time series and can be generated.
Therefore, the ensemble members of the target month are as follows: EEDP model assigns equal weights (1/n) to each ensemble member of the target month because there is no information about the climate change to assign divergent weights (Belayneh, Adamowski & Khalil, 2016).

Weighted ensemble drought prediction model
One way to represent the impact of climate change is to weight the climate index (Bradley, Habib & Schwartz, 2015). The WEDP model deals with climate index to assign weights to all ensemble members. The aim of assigning weights to each ensemble member is to reflect the change of season when predicted (Belayneh, Adamowski & Khalil, 2016). There are two types of procedure for assigning weights such as Traditional Weighting (TW) and Weighted Bootstrap Resampling (WBR) procedure.
Weighted ensemble drought prediction model algorithm. The procedure for generating ensemble members in WEDP model is similar to that of the EEDP model Eq. (9). According to Zhang, Dong & Chen (2019), WEDP model algorithm is divided into the following steps: • The appropriate climate index should be chosen based on the strong correlation with the target month of precipitation (P i,j,k ).
• Assign weights to climate indices through linear a regression model to simplify them into one integrated climate predictor (ICP) ϕ.
• Calculate the weight of each ensemble member based on the similarity of climate indices between ICP of observed and predicted year.
Traditional weighting procedure. The impact of climate indices can be signified by an unequal weighting of the ensemble members. Traditional weighting procedure utilizes the Gaussian kernel K (y) to assign a weight to each ensemble member which is based on the similarity of climate indices between the ICP of predicted and observed year (Belayneh, Adamowski & Khalil, 2016;Najafi, Moradkhani & Piechota, 2012;Zhang, Dong & Chen, 2019). A Gaussian kernel function is used to calculate the climate distance between the ICP of predicted and observed year. Therefore, the weight w i for ith ensemble member is as follows: where, ϕ i is the ICP of observed year, whereas ϕ n+1 is ICP of predicted year, which is obtained from a linear regression model. The Gaussian kernel function is defined as: The ' h' is the bandwidth kernel parameter and h = σ ε (σ ε represents the variance of error between the ICP of observed year and ICP of the predicted year).
Weighted bootstrap resampling procedure. Smith & Gelfand (1992) proposed the simplest resampling approach known as Weighted Bootstrap Resampling (WBR) procedure. The procedure of WBR is based on random variates via acceptance and rejection method. The procedure of WBR is as follows: Suppose ϕ i i = 1,2,3,..,n are samples from g and f .
and then, where, q i is the ith ratio between f and g . Therefore, w i , i = 1,2,...,n are the n weights which are corresponding to ' n' ensemble members and by definition n i=1 ω i = 1. • Use the copula function to fit the bivariate distribution between the ICP (ϕ) and target month of precipitation (P i,j,k ), i.e., C(F ϕ (ϕ),F P i,j,k (P i,j,k )).

Conditional ensemble drought prediction model
• Generate conditional precipitation events 1,000 times given the ICP (ϕ) of the observed year through the bivariate copula distribution.
• Use the conditional precipitation events i.e., P 1 i,j,k ,P 2 i,j,k ,...,P 1000 i,j,k to generate the ensemble members of the target month.
Copula ensemble. The basic aim of using the appropriate copula family is to define the dependency structure between ICP and P i,j,k . The bivariate copula distribution between ICP (ϕ) and target month of precipitation (P i,j,k ) is as follows: These, F ϕ (ϕ) and F P i,j,k (P i,j,k ) are the marginal distribution of copula function. Marginal distribution has ranged between 0 and 1 because, F ϕ (ϕ) and F P i,j,k (P i,j,k ) are the CDF of ICP and P i,j,k respectively.
According to Zhang, Dong & Chen (2019), the conditional density function is as follows: where, f (ϕ,P i,j,k ) is the joint density distribution between ICP and P i,j,k . The density f (ϕ) is the marginal density distribution of ICP. Therefore, the 1,000 times conditional precipitation events are obtained from conditional density Eq. (15), i.e., P 1 i,j,k ,P 2 i,j,k ,...,P 1000 i,j,k .
Finally, the conditional ensemble members are generated by the SPI approach described in section 2.2 on the conditional precipitation events Eq. (16), i.e.,

Copula
In 1959, the general notation of copula was first introduced (Ward, Glorig & Sklar, 1959). A major feature of the Copula function is to examine the dependence structure between the random variables. Following Kao & Govindaraju (2007) and Zhang & Singh (2006), the one-parameter Archimedean copula is expressed as, The notation ψ −1 (·) represents the copula generator and it is a convex decreasing function. The notation φ is the parameter of Copula function and φ(1) = 0. In the present study, four families of Archimedean copulas such as Clayton, Joe, Frank and Gumbel are used.

Performance measurement tools
Accuracy measure tools. Various tools are used to evaluate the accuracy of predictions (Chen et al., 2014;Pai & Lin, 2005). There are several accuracy measurements tools such as Mean Absolute Error (MAE), Mean Square Error (MSE) and Root Mean Square Error (RMSE), which are scale-dependent measurement tools. Mean Absolute Percentage Error (MAPE) is widely used for measuring the accuracy of predictions of time series models, it is independent of scale. MAPE measures the magnitude of error in term of percentage. According to Parmar & Bhardwaj (2014) and Swanson, Tayman & Bryan (2011), the accuracy measurement tools are described as follows. i.e., where, y obs,i is the observed value, y pred,i is the predicted value and n is the length of time series. Normalize Root Mean Square Error (NRMSE) is a non-dimensional form of RMSE and it is the standardized disaggrement between observed and predicted value. It is used to test the predictive accuracy between observed and predicted value. According to Zhang, Dong & Chen (2019), NRMSE is defined: where, y obs,max and y obs,min are the maximum and minimum value in the observed time series. The smallest value of MAE, MSE, RMSE, MAPE and NRMSE is an indication of better prediction.
Uncertainty measure tool (average bandwidth). In probabilistic prediction, Chen et al. (2014) presented Average Bandwidth (AB) to measure the prediction uncertainty. The AB provides very important information about the prediction uncertainty., AB is defined as follows: where, Q 1−a/2 i and Q a/2 i is an indication of quantiles result of ith the year and a is the level of significance. Therefore, the smallest value of AB is an indication of better prediction.

RESULTS
Our first objective is to define and monitor the short-term wet and dry events, which is done by SPI-3. In the estimation phase of SPI-3, several distributions are fitted on 3 months aggregated precipitation series (P-3). A detailed description of fitting suitable probability distribution is illustrated in Table 3.

Notes.
The b, r and k are the scale, location and rate parameter respectively. As well as both a and c are the shape parameter.
After this, the cumulative distribution function (CDF) is computed. Later, the CDF of P-3 is utilized for the construction of SPI-3 by following the procedure described in section 2.2.
According to Table 2, all meteorological stations are distressed from various categories of wet and dry events e.g., normal, moderate, severe and extreme. But the occurrence percentage of extreme drought events have appeared in the Multan, Bahawalpur, Barkhan and Khanpur is 1.44%, 0.57%, 2.59% and 1.71% respectively. Whereas, the occurrence percentage of extremely wet events are 2.3%, 1.72%, 0.86% and 2.86% for Multan, Bahawalpur, Barkhan and Khanpur, respectively.
In the ESP procedure, the prediction of Temperature (T), Precipitation (P), P-3 and SPI-3 are needed. In this regard, the T, P, P-3 and SPI-3 are divided into two sets known as a training and testing set. The entire series of T, P, P-3 and SPI-3 are from the year 1990 to 2018. The length of the series in a training set is from the year 1990 to 2011. Whereas, the length of the series in the testing set is from the year 2012 to 2018. The prediction accuracy is measured by several tools i.e., MSE, RMSE, MAE and MAPE. Hence, the detailed summary of time series model fitting is described in Table 4.
The selection of P-3 is based on the correlation between P-3 and SPI-3 for all stations. Therefore, the correlation between P-3 and SPI-3 has been found strong in all stations i.e., Multan (0.90), Bahawalpur (0.93), Barkhan (0.93) and Khanpur (0.77). Due to this strong correlation, which permits the reliable monthly predictions of P-3 potentially, which improves the monthly SPI-3 prediction. The selection of the target and the current month has been based on average low precipitation and high temperature of all meteorological stations (Fig. 2). Thus, resulting from the time series models, the observed and predicted values of SPI-3, T and P-3 kept as the current month (April). Whereas, the values of SPI-3, T and P-3 of May, June and July kept as target month. After this, the ensemble members of the target month (i.e., May, June and July) are generated by applying the procedure of initial mean prediction described in section 2.3.1. Therefore, seven ensemble members are generated because the test set consists of seven years of observation, i.e., (2012)(2013)(2014)(2015)(2016)(2017)(2018). In the EEDP model, the equal weights i.e., (1/7 = 0.1428) are distributed to all ensemble members. After this, the ensemble means the prediction of EEDP model is obtained by taking the average of ensemble member. The detailed results of EEDP model are presented in Table 5. Table 4 Two-time series models are applied to the P, P-3, T and SPI-3 of all meteorological stations. The term non-seasonal represents the parameters of the ARIMA model while the term seasonal represents the parameters of the SARIMA model.  Table 5 summarizes the statistical results of accuracy and uncertainty of EEDP. The MAE, MSE, RMSE and NRMSE indicate the amount of error between observed and predicted SPI-3 for all the target months of meteorological stations. The prediction uncertainty has been tested by AB at the 5% level of significance. The smallest value of AB indicates less uncertainty in the prediction for all target months of meteorological stations.
Eleven climates indices (sea surface temperature) have been used for weighting the ensemble members of target months (i.e., May, June and July) which is described in Table 2. Initially, the correlation is computed between the climate index and P 4,5,6 . Then, three highly correlated climate indices are selected for each target month. According to Bradley, Habib & Schwartz (2015), select only one climate index to represent climate information. Thus, the LR model has been used to simplify the climate indices into one ICP. Our purpose of climate index prediction is to construct the ICP of observed and predicted year. Further, the estimation of weighting parameters in the LR model is based on the least sum of squared of error between the climate index and P 4,5,6 . Two weighting procedures i.e., TW and WBR are used in WEDP model to assign the weights to seven ensemble members. The detailed results of WEDP are presented in Table 6. Table 6 showed that the numerical values of MAE, MSE, RMSE and NMSE of WBR procedure are lower as compared to TW procedure for all target months. While the AB indicates that the uncertainty in the prediction of SPI-3 by WBR procedure is higher as compared to TW procedure for some target months.
According to Zhang, Dong & Chen (2019), the log-normal distribution is fitted on the marginal distribution of observed ICP and P 4,5,6 in the construction of the CEDP model. The parameters estimation of log-normal distribution is done by MLE. Their specified distribution is tested by KS and AD test at the 5% level of significance which satisfies Table 7 Log-normal distribution is fitted on the ICP and P 4,5,6 . The CV represents the critical value of KS and AD test. Notes. b indicates scale parameter and r is the location parameter the specified distribution. Hence, the description of the log-normal fitting is presented in Table 7. However, the estimated parameters are used for the computation of CDF of ICP and P 4,5,6 . Besides, four copula families are used i.e., Frank, Joe, Gumbel and Clayton copula for fitting the bivariate distribution between the CDF of observed ICP and P 4,5,6 . The parameter estimation of bivariate copula distribution is done by MLE. Also, the choice of suitable copula family is based on the lowest value of BIC. Therefore, the detailed numerical description of copula fitting is presented in Table 8.
After this, 1,000 times conditionally precipitation events of the P 4,5,6 given the observed ICP are generated. These conditional precipitation events are transformed into SPI by fitting the suitable probability distribution. Their proper fitting is tested through the KS and AD test at the 5% level of significance. Among these appropriate candidate probability distributions, the suitable probability distribution is selected based on the lowest value of BIC. The description of suitable fitted probability distribution on conditional precipitation events is described in Table 9. Lately, CDF of the suitable fitted probability distribution is transformed into 1,000-time conditional ensemble drought member. Finally, the CEDP is obtained by taking the average of conditional ensemble drought member. The prediction accuracy is tested by several tools i.e., MSE, RMSE, MAE and NRMSE. Moreover, the prediction uncertainty is checked at the 5% level of significance. The detailed description of the CEDP model is described in Table 10.
From Table 10, it can be observed that numerical values of prediction accuracy (i.e., MAE, MSE, RMSE and NRMSE) and uncertainty (i.e., AB) indicate better prediction.

DISCUSSION
It is challenging to define the features of drought i.e., beginning, ending, intensity and duration (Vicente-Serrano, Beguería & López-Moreno, 2010). Initially, the SPI-3 is used for the definition and monitoring of wet and dry events because the SPI is used in operational monitoring systems in various countries and also acknowledged by the World Meteorological Organization (WMO). It is observed from the temporal plots (Fig. 3), most of the meteorological stations have severely been damaged by various categories of wet and dry events. The highest percentage of extremely dry events is in Barkhan (2.59%) as compared to Khanpur (1.71%), Multan (1.44%) and Bahawalpur (0.57%). In contrast, the high percentage of extremely wet events are found in Khanpur (2.86%) as compared to Multan (2.3%), Bahawalpur (1.72%) and Barkhan (0.86%). Drought is a complicated phenomenon, so one indicator (i.e., precipitation) may be inadequate to explain the characteristics of drought (AghaKouchak, 2015). The prediction of drought based on the SPI may be insufficient to provide information to overcome drought, so that premature action can be taken. Furthermore, various variables, such as precipitation, climate indices and temperature can be involved in drought conditions in several ways and should also be taken into account in the evaluation and prediction of drought (Hao, Hao & Singh, 2016). In this regard, ensemble models were used which are based on the concept of ESP model that incorporates the weather information in the prediction. In this paper, three ensemble models are used for the drought prediction. Further, two weighting procedures i.e., TW and WBR are discussed in the WEDP model. Among these three models, the EEDP model does not use climate information whereas WEDP and CEDP models incorporate climate information. However, we assigned equal weights to seven ensemble members. The prediction of drought with the incorporation of climate indices is essential and effective. Eleven climates indices are selected based on correlation with target months of precipitation (Table 11). Then highest correlated climate indices are further simplified into ICP by LR model as the LR model contemplates the importance of climate index to simplifies into climate index. Najafi, Moradkhani & Piechota (2012) were utilized the Spearmen correlation and PCA to consider the predictors. There is one problem that the PCA obtain the information from climate indices to provide in several components. So, it is very hard to fit the joint distribution.
Two weighting procedures, i.e., TW and WBR, used the climate indices in the ensemble prediction. Then the drought is predicted by EEMD, WEDP and CEDP model. It is found that the accuracy and uncertainty which is obtained by TW procedure is like the EEDP model (Tables 5 and 6). Whereas, the results of accuracy which is obtained by WBR is almost better as compared to TW procedure and EEDP model (Tables 5 and 6). But the results of uncertainty indicate similarity to TW procedure and EEMD model. Prediction uncertainty is related to variability. The high variability in prediction may cause of reducing the reliability. The AB is based on the difference between lower and upper quantile. These lower and upper quantiles are associated with the confidence level. According to Xiong et al. (2009), the bandwidth of the prediction bounds is as narrow as possible as it captures the most important information about the prediction uncertainty. In the CEDP model, four Archimedean copula families are used to explain the dependency structure between climate indices and precipitation. The choice of best copula family is based on the lowest value of BIC, and it is observed that Joe copula has been used for most of the times (Table 8). As seen, the results of accuracy and uncertainty indicate CEDP model is better as compared to WEDP and EEDP (Tables 5, 6 and 10). The achieved results of CDEP model indicate that the climate indices are correlated with weather patterns for four meteorological stations. Table 11 The CIM represents the selected climate index month, TM represents the target month and the given correlation is representing the relationship between the P 4,5,6 and CIM. These described parameters are of the ARIMA model. In TW procedure, we always need the parameter of Gaussian kernel function for computing the weights. The estimation of Gaussian kernel parameter decreases the accuracy of prediction and required more investigation in the computation while the WBR procedure only utilizes the ICP of observed year. The ratios are obtained between ICP of observed year and corresponding ensemble member. The traditional bootstrap resampling procedure is dissimilar from WBR procedure (Burnham & Efron, 1983) because the weight is obtained by dividing the ratio further. The EEDP and WEDP model have generated the ensemble members corresponding to the historical year in series while the generation of ensemble members is not restricted to the historical year in CEDP model. Therefore, the CEDP model is more convenient as compared to the EEDP and WEDP model. But, there is a limitation in the CEDP model. The ensemble prediction by CEDP model required the accurate fit of bivariate copula distribution and probability distribution for conditional precipitation events. In contrast, ensemble prediction will not be accurate. Additionally, these techniques can be used to further improve the study. To deal with stations which are seasonally influenced by droughts, the seasonal drought analysis can be done instead of monthly drought prediction.

CONCLUSION
The present analysis consists of four meteorological stations i.e., Multan, Bahawalpur, Barkhan and Khanpur. Initially, SPI-3 was selected to define and classify the drought for all stations. It was observed that the highest percentage of extremely dry events were in Barkhan. In contrast, the high percentage of extremely wet events were found in Khanpur.
Later, three ensemble models are used for the monthly drought prediction. Further, two weighting procedures i.e., TW and WBR were included in the WEDP model. Among these three models, the performance of the CEDP model in terms of accuracy and uncertainty is better in comparison to EEDP and WEDP for all meteorological stations. The performance of the CEDP model indicates that the climate indices are correlated with a weather pattern of four stations as the construction of the CEDP model based on the correlated structure of climate index and precipitation.
• Ijaz Hussain conceived and designed the experiments, performed the experiments, prepared figures and/or tables, and approved the final draft.
• Muhammad Faisal conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft.
• Muhammad Ismail, Alaa Mohamd Shoukry, Elsayed Elsherbini Elashkar and Showkat Gani conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, and approved the final draft.

Data Availability
The following information was supplied regarding data availability: The R code and data from the Bahawalpur stations are available in the Supplemental Files.