Modeling, prediction and trend assessment of drought in Iran using standardized precipitation index

In recent decades, increasing regard has been assigned to drought events and their tremendous detrimental influences on agriculture, economy, and environment, in scheduling water management programs (Zarei et al. 2016; Hao et al. 2017). As Iran is located in an arid and semi-arid region, prediction and modeling drought in the future is vital. Prediction of drought will lead to plans and administrations that are significantly effective for preparation against drought (Salhvand & Montazeri 2013). Drought is a costly and frequent phenomenon that usually occurs after an extraordinary water deficit due to low rainfall over a large geographical area. The frequency of drought occurrence increases as the region becomes more and more arid (Zarei et al. 2016; Bahrami et al. 2017). There exists no precise definition for drought, but in general, almost all the definitions give the concept of deficit of precipitation resulting in water shortages in applications related to water usage (Wilhite & Glants 1985).

The most well-known indices are Palmer Drought Severity Index, Reconnaissance Drought Index (RDI), Standardized Precipitation Evapotranspiration Index, Multivariate Standardized Drought Index, China Z Index (CZI), and Standardized Precipitation Index (SPI) among many indices that have been studied and utilized all around the world by many scientists (Salhvand & Montazeri 2013; Azadi et al. 2015; Tajbakhsh et al. 2015; Zarei & Mahmoudi 2017; Zarei 2018). The SPI index (McKee et al. 1993) has been utilized, due to being modest, flexible, with a precise approach to calculating and the least number of meteorological parameters needed (Choubin et al. 2016). This index can be used in different time scales (3, 6, 12, 24, 48-month scale). Various studies have been done on drought analysis using SPI (Shakiba et al. 2010; Ansari et al. 2010; Tabari et al. 2011; Paulo et al. 2012; Hanafi et al. 2012; Hasanzade et al. 2013; Zehtabian et al. 2013; Moreira et al. 2014; Shah et al. 2015; Huang et al. 2016).

Zarei et al. (2017) used SPI index and Markov chain method to monitor and predict monthly drought in Iran. The results showed that, in the most synoptic stations, normal, moderately dry, and severe dry classes of drought have the highest frequency of occurrence. Roshan et al. (2018) evaluated a new application for TOPSIS: monitoring drought and wet periods in Iran using RDI index. Sadeghi & Hazabvi (2017) evaluated the spatiotemporal variation of watershed health propensity through a reliability–resilience–vulnerability based drought index (case study: Shazand watershed in Iran).

Time series datasets are a sequence taken set of data collected successively in a prolonged period, in various time scales. The natural order of time series datasets is the most striking feature of time series analysis that distinguishes them from any other analysis. Understanding the previous trend of time series leads to mathematical modeling of the observed trend for monitoring the datasets and a precise prediction in the future as well (Shirmohammadi 2013). Nowadays, analysis of time series models is widely used. Time series analysis is employed for five main reasons: descriptive, explanation, forecasting, intervention analysis, and quality control. Time series have drawn researchers' attention throughout the world (Asakereh et al. 2013; Azad Talatapeh et al. 2013; Ghahraman 2013; Deng et al. 2015; Deng & Wang 2017).

The aims of this study are to focus on: (1) evaluation of occurrence of seasonal drought severity in Iran, using SPI index; (2) modeling and prediction of seasonal drought using time series models, from 2015 up to 2019 based on seasonal SPI index data series of 1967 to 2014; (3) considering the changes trend of drought (using nonparametric trend analysis techniques), termed as: (Spearman rho test) based on observed (1967–2014) and the combination of observed and predicted (1967–2019) seasonal SPI data series. Therefore, meteorological data of 38 synoptic stations in Iran with suitable spatial distribution and various climate conditions were collected.

Iran is located geographically between 25° 03′ to 39° 47′ Northern latitude and 44° 05′ to 63° 18′ Eastern longitude, with an area of 1,648,198 square kilometers. The climate varies from sub-humid to hyper-arid in northern and southern parts of Iran, respectively (based on the United Nation Environmental Program aridity index). The annual precipitation of Iran varies from less than 100 to more than 2,000 millimeters but, in general, the average annual precipitation is about 250 millimeters. This implies there are various climates in Iran. The major amount of the annual precipitation in Iran falls in winter (Dinpashoh et al. 2004). In this study, precipitation data series required for calculating seasonal SPI index were collected from 38 synoptic stations (with different climate conditions) of the country (Figure 1), obtained from the meteorological organization service of Iran for the period 1967–2014 (http://www.irimo.ir/). Characteristics of synoptic stations are presented in Table 1.

In general, time series models have three major models: Auto-Regressive (AR), Moving-Average (MA), and Auto-Regressive–Moving-Average (ARMA). AR model divides into Yule-Walker and Burg, and MA model divides into Hannan-Rissanen and innovations. Times series models, like any other model, have some variables (Mirzavand & Ghazavi 2015). The number of variables in AR and MA models are demonstrated by p and q indices, respectively. Various species of models for time series data exist that illustrate different non-deterministic processes. Modeling of time series typically occur based on a linear technique, for instance AR, MA, ARMA, and ARIMA (ARMA with more than one difference) models have a linear base (Klose et al. 2004). These models have been widely employed in various branches of science, such as meteorology sciences, water sciences, and others (Shabani et al. 2013). In this research, to assess time series model ITSM (Interactive Time Series Modeling) software was used. This software is an interactive Windows-based menu-driven software for time series modeling and forecasting (Brockwell & Davis 1991).

Autoregressive Integrated Moving Average (ARIMA) models are one of the well-known linear models for time series modeling and predicting. ARIMA models originated from the synthesis of AR and MA models. ARIMA is used to model time series data behavior and to make predictions. ARIMA modeling uses correlational methods and could be used to model arrays that may not be observable in plotted data. In ARIMA, the future amount of a parameter is assumed to be a linear function of past observations and random errors (Jahandideh & Shirvani 2011).

A SARIMA model can be explained as ARIMA (⁠⁠,⁠,⁠) (⁠⁠,⁠,⁠)⁠, where (⁠⁠,⁠,⁠) is the non-seasonal component of the model and (⁠⁠,⁠,⁠) is the seasonal component of the model in which p is the order of non-seasonal auto-regression, d is the number of regular differencing, q is the order of nonseasonal moving average, P is the order of seasonal auto-regression, D is the number of seasonal differencing, Q is the order of seasonal moving average. The subscripted letter ‘s’ shows the length of the seasonal period. For example, in an hourly data time series s = 7, in a quarterly data s = 4, and in a monthly data s = 12 (Martinez et al. 2011).

In this stage, the best time series models that fitted the seasonal SPI data series from 1967 to 2014 were used to predict seasonal SPI data series from 2015 to 2019 using ARAR method. (The ARAR algorithm is basically the process that applies memory shortening transformation and fitting the AR model to the transformed data. It is used to predict the future data from existing sequence data. The algorithm was introduced by Brockwell & Davis (1991) and it consists of three phases throughout the process.) To evaluate the accuracy of selected model prediction Ljung–Box p-value (if p-value of Ljung–Box index was more than 0.05 randomness of data series is significant at the 5% level), residual ACF, residual PACF (according to these tests, in both charts it should not be more than 5% of maximum number of lags (40) out of the zero band) and correlation coefficient between predicted and observed seasonal SPI data during 2010–2014 were used.

In this section, after assessing the normality of seasonal SPI data series, the trends of these series were evaluated based on non-parametric statistical (Spearman's rho test) and parametric test (Pearson test). In this study, the changes trend of data series was evaluated in two stage: Stage 1, the changes trend of observed seasonal SPI data series (1967–2014) and Stage 2, changes trend of observed and predicted seasonal SPI data series (1967–2019).

In the present study, the SPI index in seasonal scale for 38 synoptic stations was calculated (Figure 2) based on the expressed approach by McKee et al. (1993). Results of SPI calculation showed that in all stations, drought severity with class 4 (normal class) and 3 (moderately wet) had the most occurrence frequency, and according to the results, extremely dry and very dry classes had the least occurrence frequency (Table 3).

After calculation of SPI indices, for stationarity of seasonal data series, Box–Cox power conversion was used to immobilize the variances of datasets of seasonal SPI indices. Difference conversion was used for removal of trend in data series and to exert the frequency impact in the SPI indices dataset. Then the set of non-stationary data was converted to stationary data and was ready for modeling (Figure 3).

In this stage, various kinds of time series models were fitted to SPI datasets and the best fitted model was chosen based on having the least AICC index (Table 4). To determine order of time series models (p and q), ACF and PACF were used (Figure 4).

According to the results of ACF and PACF, p and q orders equal to 5, 12, 13, 14, 16, 18, 22, and 25 were the best orders that were capable in the fitted time series models.

Results of time series model selection indicated that at Abadan, Arak, Bandar Lenge, Birjand, Dezful, Esfahan, Fasa, Qazvin, Gorgan, Iran Shahr, Kashan, Kermanshah, Khoramabad, Khoy, Oroomieh, Sabzevar, Sanandaj, Semnan, Tabriz, Tehran, Torbate Heydarieh, Zahedan, and Zanjan stations, MA (5) Innovations model was the best fitted time series model. At Ahvaz, Babolsar, Bushehr, Kerman, Sharekord, and Shiraz stations, MA (5) Hannan-Rissanen was the best fitted time series model. The only station fitted by MA (13) Innovations model is Bandar Abbas. Chabahar and Hamedam stations were fitted by AR (16) Burg model while at Rasht and Tabas stations, AR (18) Burg was the best fitted time series model. In Mashhad, Ramsar, Yazd and Zabol, the best fitted time series were AR (25) Burg, AR (12) Burg, AR (14) Burg, and AR (22) Burg, respectively.

In order to validate the model dataset of seasonal SPI indices, five years (20 seasons) observed data were available, and were predicted using ARAR model. To compare predicted and observed seasonal SPI data, correlation coefficient was used (Table 5). According to the results, in all stations correlation coefficients between observed and predicted data series of seasonal SPI from 2010 to 2014 were significant at 5% level. In stations with humid, sub-humid, semi-humid, and semi-arid climate condition and stations with high relative humidity such as Chabahar, Gorgan, and Orumieh stations, R was less than the correlation coefficient in stations with arid and hyper-arid climate condition with low relative humidity such as Zabol, Yazd, and Tabas stations.

In this stage, to predict seasonal SPI indices for five years after 2014 (from 2015 to 2019), the ARAR approach was employed. Predicted seasonal SPI indices and occurrence frequency are presented in Tables 6 and 7 and Figure 5.

Results of observed, validated, and predicted seasonal SPI data series are indicated in Figure 5 for some stations.

According to the results, in all the synoptic stations, drought severity class 4 (normal) had the most frequency (in 71% of stations occurrence frequency of class 4 was 100%, in 21% of stations occurrence frequency of class 4 was more than 80%, occurrence frequency of class 4 in other stations was more than 60%). In addition, classes 1 and 7 (extremely wet and extremely dry) had the least occurrence frequency (Tables 6 and 7).

In the next stage, capability of time series models (based on ARAR method) for prediction of seasonal SPI indices were evaluated using Ljung–Box p-value and residual ACF and PACF methods (Figure 6 and Table 8).

Results of residual ACF and residual PACF demonstrated that in all the stations, capability of time series models was suitable at the 5% significance level because in both charts, residual ACF and PACF, no more than 5% of maximum number of lags (40) was out of the 95% confidence band (Figure 5). According to the Ljung–Box p-value, in 82% of stations randomness of data series was suitable at 5% significance level.

Finally, trend of observed seasonal SPI data series from 1967 to 2014 and the whole observed and predicted seasonal SPI data series from 1967 to 2019 simultaneously was assessed statistically at 5% significance level, using Spearman's rho and Pearson, non-parametric and parametric tests, respectively (Table 9).

Based on the statistical test for observed seasonal SPI data series, in none of the stations except Gorgan, Sanandaj, and Tabriz was any significant trend observed. However, these three stations (Gorgan, Sanandaj, and Tabriz) had a significant decreasing trend in the values of their observed seasonal SPI data series. The same results were shown in the statistical test for the combination of observed and predicted seasonal SPI data series. Apart from having a significant trend, in all study stations except Chabahar, Esfahan, Qazvin, Kashan, Ramsar, and Tabas, a slight descending trend was observed in seasonal SPI data series. On the other hand, results of parametric and nonparametric test for the whole observed and predicted seasonal SPI data series from 1967 to 2014 indicated a gradual negative trend in the seasonal SPI data series in all stations, except Chabahar, Esfahan, Kashan, Ramsar, and Tabas. Decreasing trend in SPI index is due to the reduction of precipitation, which can be dangerous for different sections such as agriculture, ecosystems, water resources, and others.

Since drought affects human activities, agriculture, economy, environment, water resources management, evaluation and recognition of drought behavior will help society to overcome this disastrous phenomenon. In the present study, the seasonal SPI data series of 38 synoptic stations of Iran were calculated, modeled, and predicted. Finally, change trend of calculated seasonal SPI data series, based on non-parametric and parametric statistical test was evaluated. Results of this research demonstrated that in the majority of stations throughout the study, MA (5) Innovation and MA (5) Hannan-Rissenen were fitted better than other models with minimum AICC index. Results of predicted seasonal SPI indices indicated that drought severity class 4 (normal) had the most and classes 1 and 7 (extremely wet and extremely dry) had the least occurrence frequency. Results of validation indicated that the model prediction performed very cautiously and the lower band and upper band of the predicted SPI index and validated SPI index varies from the observed SPI index. Based on changes trend of data series test, in the vast majority of the synoptic stations in the study area, with various drought severity classes, a slight negative trend was seen in the seasonal SPI data series. In this study, changes trend of drought based on observed data (1967–2014) and changes trend of drought based on observed and predicted data (1967–2019) were evaluated, and the effect of SPI index from 2014 to 2019 in trend of changes in drought severity resulted (in other works mainly changes trend of observed data were evaluated). Zarei et al. (2016) employed the RDI index and time series model for modeling and prediction of the seasonal drought in Tehran. Their results demonstrated that MA (5) Hannan-Rissenen was the best fitted model, contrary to the obtained result of this study where MA (5) Innovation fitted the dataset as well as possible. The difference may be because of the various indices. Jahandideh & Shirvani (2011) used SPI and time series models to forecast drought in Fars province (Shiraz and Fasa synoptic stations) using precipitation data from 1967 to 2005. According to the results, the SARIMA model with the minimum AICC index was selected as the best model. However, the current research rejects their results in some cases. The difference may be because of the various time periods. According to the fact that the main occupation of Iranian people is agriculture and livestock farming, drought and its side effects have more detrimental effects on Iran. Therefore, drought assessment and drought prediction can be more helpful in reducing drought impacts. Finally, it is suggested that other research should be done using other models to predict drought index to confirm the results of this study. According to the results of validation, time series model prediction implemented very cautiously and the lower band and upper band of the predicted SPI index were close to the average of observed SPI data series.