## Abstract

Drought as an exigent natural phenomenon, with high frequency in arid and semi-arid regions, leads to enormous damage to agriculture, economy, and environment. In this study, the seasonal Standardized Precipitation Index (SPI) drought index and time series models were employed to model and predict seasonal drought using climate data of 38 Iranian synoptic stations during 1967–2014. In order to model and predict seasonal drought ITSM (Interactive Time Series Modeling) statistical software was used. According to the calculated seasonal SPI, within the study area, drought severity classes 4 and 3 had the greatest occurrence frequency, while classes 6 and 7 had the least occurrence frequency. Results indicated that the best fitted models were Moving-Average or MA (5) Innovations and MA (5) Hannan-Rissenen, with 60.53 and 15.79 percentage, respectively. On the other hand, results of the prediction as well, indicated that drought class 4 with the highest percentages, was the most abundant class over the study area and drought class 7 was the least frequent class. According to results of trend analysis, without attention to significance of them, observed seasonal SPI data series (1967–2014), in 84.21% of synoptic stations had a negative trend, but this percentage changes to 86.84% when studying the combination of observed and predicted simultaneously (1967–2019).

## INTRODUCTION

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.

## MATERIALS

### Study area and data collection

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.

Station | Latitude (N) | Longitude (E) | Elevation (meter from sea level) | Average of precipitation (mm, year^{−1}) | Average of ET_{0} (mm.day^{−1})^{a} | Aridity index | Climate condition |
---|---|---|---|---|---|---|---|

Abadan | 30.367 | 48.250 | 6.25 | 158.48 | 6.82 | 0.065 | Arid |

Ahvaz | 31.333 | 48.667 | 22.34 | 234.76 | 6.79 | 0.096 | Arid |

Arak | 34.100 | 49.767 | 1,708.18 | 325.11 | 3.78 | 0.239 | Semi-arid |

Babolsar | 36.717 | 52.650 | −21.31^{b} | 922.00 | 2.55 | 1.004 | Humid |

Bandar Abbas | 27.220 | 56.367 | 10.19 | 175.51 | 5.30 | 0.092 | Arid |

Bandar Lenge | 26.533 | 54.833 | 23.16 | 134.70 | 5.40 | 0.069 | Arid |

Birjand | 32.867 | 59.200 | 1491.45 | 163.76 | 5.23 | 0.087 | Arid |

Bushehr | 28.983 | 50.833 | 20.12 | 254.17 | 5.12 | 0.138 | Arid |

Chabahar | 25.283 | 60.617 | 8.24 | 113.66 | 4.21 | 0.075 | Arid |

Dezful | 32.400 | 48.383 | 143.25 | 392.62 | 4.77 | 0.229 | Semi-arid |

Esfahan | 32.617 | 51.667 | 1,550.39 | 130.98 | 4.36 | 0.083 | Arid |

Fassa | 28.967 | 53.683 | 1,288.48 | 147.20 | 4.69 | 0.087 | Arid |

Ghazvin | 36.250 | 50.050 | 1,279.13 | 321.27 | 3.79 | 0.235 | Semi-arid |

Gorgan | 36.850 | 54.267 | 13.37 | 573.11 | 2.81 | 0.567 | Sub-humid |

Hamedan | 35.200 | 48.717 | 1,697.22 | 321.88 | 4.09 | 0.219 | Semi-arid |

Iran Shahr | 27.200 | 60.700 | 591.24 | 111.53 | 6.28 | 0.049 | Hyper-arid |

Kashan | 33.983 | 51.450 | 982.31 | 136.41 | 3.40 | 0.111 | Arid |

Kerman | 30.250 | 56.967 | 1,753.47 | 139.34 | 5.41 | 0.072 | Arid |

Kermanshah | 34.350 | 47.150 | 1,318.47 | 447.09 | 4.35 | 0.285 | Semi-arid |

Khoram Abad | 33.433 | 48.283 | 1,147.32 | 499.73 | 4.09 | 0.339 | Semi-arid |

Khoy | 38.550 | 44.967 | 1,103.14 | 292.54 | 3.09 | 0.263 | Semi-arid |

Mashhad | 36.267 | 59.633 | 999.35 | 255.87 | 3.96 | 0.179 | Arid |

Oroomieh | 37.533 | 45.083 | 1,315.15 | 332.33 | 3.20 | 0.288 | Semi-arid |

Ramsar | 36.900 | 50.667 | −20.25^{b} | 1208.89 | 2.24 | 1.499 | Humid |

Rasht | 37.250 | 49.600 | −6.21^{b} | 1354.60 | 2.32 | 1.622 | Humid |

Sabzevar | 36.200 | 57.717 | 977.35 | 196.10 | 5.34 | 0.102 | Arid |

Sanandaj | 35.333 | 47.000 | 1,373.35 | 445.60 | 3.96 | 0.313 | Semi-arid |

Semnan | 35.583 | 53.550 | 1,130.42 | 142.10 | 4.10 | 0.096 | Arid |

Shahre Kord | 32.283 | 50.850 | 2,048.12 | 332.43 | 3.54 | 0.261 | Semi-arid |

Shiraz | 29.533 | 52.600 | 1,484.32 | 324.88 | 4.88 | 0.185 | Arid |

Tabass | 33.600 | 56.917 | 711.46 | 82.80 | 4.79 | 0.048 | Hyper-arid |

Tabriz | 38.083 | 46.283 | 1,361.29 | 272.80 | 4.18 | 0.181 | Arid |

Tehran | 35.683 | 51.317 | 1,190.36 | 239.76 | 4.85 | 0.137 | Arid |

Torbat Hydarieh | 35.267 | 59.217 | 1,450.19 | 270.39 | 3.89 | 0.193 | Arid |

Yazd | 31.900 | 54.283 | 1,237.28 | 58.89 | 5.50 | 0.030 | Hyper-arid |

Zabol | 31.033 | 61.483 | 489.31 | 56.21 | 8.62 | 0.018 | Hyper-arid |

Zahedan | 29.467 | 60.883 | 1,370.12 | 79.20 | 5.81 | 0.038 | Hyper-arid |

Zanjan | 36.683 | 48.483 | 1,663.18 | 306.68 | 3.47 | 0.246 | Semi-arid |

Station | Latitude (N) | Longitude (E) | Elevation (meter from sea level) | Average of precipitation (mm, year^{−1}) | Average of ET_{0} (mm.day^{−1})^{a} | Aridity index | Climate condition |
---|---|---|---|---|---|---|---|

Abadan | 30.367 | 48.250 | 6.25 | 158.48 | 6.82 | 0.065 | Arid |

Ahvaz | 31.333 | 48.667 | 22.34 | 234.76 | 6.79 | 0.096 | Arid |

Arak | 34.100 | 49.767 | 1,708.18 | 325.11 | 3.78 | 0.239 | Semi-arid |

Babolsar | 36.717 | 52.650 | −21.31^{b} | 922.00 | 2.55 | 1.004 | Humid |

Bandar Abbas | 27.220 | 56.367 | 10.19 | 175.51 | 5.30 | 0.092 | Arid |

Bandar Lenge | 26.533 | 54.833 | 23.16 | 134.70 | 5.40 | 0.069 | Arid |

Birjand | 32.867 | 59.200 | 1491.45 | 163.76 | 5.23 | 0.087 | Arid |

Bushehr | 28.983 | 50.833 | 20.12 | 254.17 | 5.12 | 0.138 | Arid |

Chabahar | 25.283 | 60.617 | 8.24 | 113.66 | 4.21 | 0.075 | Arid |

Dezful | 32.400 | 48.383 | 143.25 | 392.62 | 4.77 | 0.229 | Semi-arid |

Esfahan | 32.617 | 51.667 | 1,550.39 | 130.98 | 4.36 | 0.083 | Arid |

Fassa | 28.967 | 53.683 | 1,288.48 | 147.20 | 4.69 | 0.087 | Arid |

Ghazvin | 36.250 | 50.050 | 1,279.13 | 321.27 | 3.79 | 0.235 | Semi-arid |

Gorgan | 36.850 | 54.267 | 13.37 | 573.11 | 2.81 | 0.567 | Sub-humid |

Hamedan | 35.200 | 48.717 | 1,697.22 | 321.88 | 4.09 | 0.219 | Semi-arid |

Iran Shahr | 27.200 | 60.700 | 591.24 | 111.53 | 6.28 | 0.049 | Hyper-arid |

Kashan | 33.983 | 51.450 | 982.31 | 136.41 | 3.40 | 0.111 | Arid |

Kerman | 30.250 | 56.967 | 1,753.47 | 139.34 | 5.41 | 0.072 | Arid |

Kermanshah | 34.350 | 47.150 | 1,318.47 | 447.09 | 4.35 | 0.285 | Semi-arid |

Khoram Abad | 33.433 | 48.283 | 1,147.32 | 499.73 | 4.09 | 0.339 | Semi-arid |

Khoy | 38.550 | 44.967 | 1,103.14 | 292.54 | 3.09 | 0.263 | Semi-arid |

Mashhad | 36.267 | 59.633 | 999.35 | 255.87 | 3.96 | 0.179 | Arid |

Oroomieh | 37.533 | 45.083 | 1,315.15 | 332.33 | 3.20 | 0.288 | Semi-arid |

Ramsar | 36.900 | 50.667 | −20.25^{b} | 1208.89 | 2.24 | 1.499 | Humid |

Rasht | 37.250 | 49.600 | −6.21^{b} | 1354.60 | 2.32 | 1.622 | Humid |

Sabzevar | 36.200 | 57.717 | 977.35 | 196.10 | 5.34 | 0.102 | Arid |

Sanandaj | 35.333 | 47.000 | 1,373.35 | 445.60 | 3.96 | 0.313 | Semi-arid |

Semnan | 35.583 | 53.550 | 1,130.42 | 142.10 | 4.10 | 0.096 | Arid |

Shahre Kord | 32.283 | 50.850 | 2,048.12 | 332.43 | 3.54 | 0.261 | Semi-arid |

Shiraz | 29.533 | 52.600 | 1,484.32 | 324.88 | 4.88 | 0.185 | Arid |

Tabass | 33.600 | 56.917 | 711.46 | 82.80 | 4.79 | 0.048 | Hyper-arid |

Tabriz | 38.083 | 46.283 | 1,361.29 | 272.80 | 4.18 | 0.181 | Arid |

Tehran | 35.683 | 51.317 | 1,190.36 | 239.76 | 4.85 | 0.137 | Arid |

Torbat Hydarieh | 35.267 | 59.217 | 1,450.19 | 270.39 | 3.89 | 0.193 | Arid |

Yazd | 31.900 | 54.283 | 1,237.28 | 58.89 | 5.50 | 0.030 | Hyper-arid |

Zabol | 31.033 | 61.483 | 489.31 | 56.21 | 8.62 | 0.018 | Hyper-arid |

Zahedan | 29.467 | 60.883 | 1,370.12 | 79.20 | 5.81 | 0.038 | Hyper-arid |

Zanjan | 36.683 | 48.483 | 1,663.18 | 306.68 | 3.47 | 0.246 | Semi-arid |

^{a}ET0 calculated based on FAO Penman–Monteith (FAO-56).

^{b}Mentioned synoptic stations located on the coasts of the Caspian Sea. The elevation of this area is lower than sea level.

## METHODS

### Standardized Precipitation Index

*et al.*(1993) reported the SPI index to monitor and evaluate drought. Fitting a gamma probability density to the total precipitation is essential, in order to assess SPI indices. The gamma distribution function (г(

*α*)) is fitted to the dataset of precipitation involving a shape factor and a scale factor, termed as

*α*and

*β*, respectively, which require an estimation for each year and time scale, hence the maximum likelihood solutions are employed. If the amount of precipitation is represented by

*x*, the cumulative probability (

*G*(

*x*)), can be assessed by:

*x*is equal to zero, the gamma function is undefined and the precipitation distribution may contain zeros (Mckee

*et al.*1993), and the definition of cumulative probability (

*H(x)*) changes to (Zarei & Mahmudi 2017): where

*q*is the probability of a zero. The cumulative probability is then transformed to the standard normal random variable

*Z*with mean zero and variance of one, which is the value of SPI. SPI is categorized based on their range values as shown in Table 2.

SPI domain | Drought class | Class number |
---|---|---|

+2 ≥ | Extremely wet | 1 |

1.5 < SPI < 1.99 | Very wet | 2 |

1.0 < SPI < 1.49 | Moderately wet | 3 |

−0.99 < SPI < 0.99 | Near normal | 4 |

−1.0 < SPI < −1.49 | Moderately dry | 5 |

−1.5 < SPI < −1.99 | Very dry | 6 |

−2 ≤ | Extremely dry | 7 |

SPI domain | Drought class | Class number |
---|---|---|

+2 ≥ | Extremely wet | 1 |

1.5 < SPI < 1.99 | Very wet | 2 |

1.0 < SPI < 1.49 | Moderately wet | 3 |

−0.99 < SPI < 0.99 | Near normal | 4 |

−1.0 < SPI < −1.49 | Moderately dry | 5 |

−1.5 < SPI < −1.99 | Very dry | 6 |

−2 ≤ | Extremely dry | 7 |

### Time series models

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

### AR model

*t*

*+*1)th period is dependent on the present

*tth*period magnitude and those preceding values, then for such a series, the observed sequence

*X*

_{1},

*X*

_{2},,

*X*is used to fit an AR model (Equation (3)): where ,, … , are model parameters and coefficients,

_{t}*e*is random error and

_{t}*Z*is the random component of the data that follows a normal distribution with mean 0 (Jahandideh & Shirvani 2011; Dastorani

_{t}*et al.*2016).

### MA model

*e*is random error, and

_{t}*Z*is the random component of the data that follows a normal distribution with mean 0 (Jahandideh & Shirvani 2011; Dastorani

_{t}*et al.*2016).

### ARMA model

*t*(Jahandideh & Shirvani 2011; Dastorani

*et al.*2016; Zarei 2018).

### ARIMA model

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

### SARIMA model

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

### Stationarity of SPI data series

*z*is the initial value in time

_{i}*i*,

*λ*conversion parameter and

*z*is the converted value.

_{i}^{λ}### Time series model selection

*k*lag as a function for expressing time depending on a time series structure is as follows: where

*ρ*is the value of time series autocorrelation function with

_{k}*k*lags,

*z*and

_{i}*z*are the values of variables or time series data within phase time

_{i+k}*I*, and phase with lag time

*k*and is average value corresponding to variables.

*k*, the equation for PACF becomes: where is the value of time series PACF with the lag of

*k*.

*n*) and parameters used in modeling process (

*k*) and the abbreviation form of mean square error as presented in Equation (9): The less AICC index, the better the model fits the dataset.

### Prediction and validation of SPI index

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.

### Trend analysis of seasonal SPI data series

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

## RESULTS AND DISCUSSION

### SPI index

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

Drought class | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|

Abadan | 0 | 5.21 | 11.98 | 82.81 | 0 | 0 | 0 |

Ahvaz | 0.52 | 2.60 | 16.67 | 80.21 | 0 | 0 | 0 |

Arak | 0.52 | 1.04 | 12.50 | 64.06 | 21.88 | 0 | 0 |

Babolsar | 0.52 | 6.77 | 10.42 | 64.06 | 11.46 | 5.21 | 1.56 |

Bandar Abbas | 2.60 | 5.73 | 9.38 | 82.29 | 0 | 0 | 0 |

Anzali | 2.08 | 4.69 | 9.90 | 64.58 | 11.98 | 5.21 | 1.56 |

Lenge | 1.56 | 4.69 | 13.02 | 80.73 | 0 | 0 | 0 |

Birjand | 0.52 | 6.25 | 10.42 | 82.81 | 0 | 0 | 0 |

Bushehr | 1.04 | 3.13 | 14.58 | 81.25 | 0 | 0 | 0 |

Chabahar | 1.56 | 6.25 | 8.33 | 83.85 | 0 | 0 | 0 |

… | … | … | … | … | … | … | … |

Shiraz | 0.52 | 5.21 | 12.50 | 62.50 | 19.27 | 0 | 0 |

Tabas | 1.04 | 5.21 | 8.33 | 85.42 | 0 | 0 | 0 |

Tabriz | 0.52 | 3.65 | 8.85 | 71.35 | 6.25 | 5.73 | 3.65 |

Tehran | 0 | 3.13 | 10.94 | 66.15 | 8.85 | 10.94 | 0 |

Torbat | 0 | 4.17 | 10.94 | 64.06 | 20.83 | 0 | 0 |

Yazd | 1.56 | 3.65 | 12.50 | 82.29 | 0 | 0 | 0 |

Zabol | 1.04 | 6.77 | 8.33 | 83.85 | 0 | 0 | 0 |

Zahedan | 1.56 | 4.17 | 11.98 | 82.29 | 0 | 0 | 0 |

Zanjan | 0 | 1.56 | 12.50 | 68.75 | 6.25 | 7.29 | 3.65 |

Drought class | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|

Abadan | 0 | 5.21 | 11.98 | 82.81 | 0 | 0 | 0 |

Ahvaz | 0.52 | 2.60 | 16.67 | 80.21 | 0 | 0 | 0 |

Arak | 0.52 | 1.04 | 12.50 | 64.06 | 21.88 | 0 | 0 |

Babolsar | 0.52 | 6.77 | 10.42 | 64.06 | 11.46 | 5.21 | 1.56 |

Bandar Abbas | 2.60 | 5.73 | 9.38 | 82.29 | 0 | 0 | 0 |

Anzali | 2.08 | 4.69 | 9.90 | 64.58 | 11.98 | 5.21 | 1.56 |

Lenge | 1.56 | 4.69 | 13.02 | 80.73 | 0 | 0 | 0 |

Birjand | 0.52 | 6.25 | 10.42 | 82.81 | 0 | 0 | 0 |

Bushehr | 1.04 | 3.13 | 14.58 | 81.25 | 0 | 0 | 0 |

Chabahar | 1.56 | 6.25 | 8.33 | 83.85 | 0 | 0 | 0 |

… | … | … | … | … | … | … | … |

Shiraz | 0.52 | 5.21 | 12.50 | 62.50 | 19.27 | 0 | 0 |

Tabas | 1.04 | 5.21 | 8.33 | 85.42 | 0 | 0 | 0 |

Tabriz | 0.52 | 3.65 | 8.85 | 71.35 | 6.25 | 5.73 | 3.65 |

Tehran | 0 | 3.13 | 10.94 | 66.15 | 8.85 | 10.94 | 0 |

Torbat | 0 | 4.17 | 10.94 | 64.06 | 20.83 | 0 | 0 |

Yazd | 1.56 | 3.65 | 12.50 | 82.29 | 0 | 0 | 0 |

Zabol | 1.04 | 6.77 | 8.33 | 83.85 | 0 | 0 | 0 |

Zahedan | 1.56 | 4.17 | 11.98 | 82.29 | 0 | 0 | 0 |

Zanjan | 0 | 1.56 | 12.50 | 68.75 | 6.25 | 7.29 | 3.65 |

### Stationarity of SPI series

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

### Selection of the best time series model

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

Station | AICC Index | Method of time series model | Non-significant θ or φ | Fitted model |
---|---|---|---|---|

Abadan | 302.88 | MA (5) Innovations | Θ2 and θ3 | X(t) = Z(t) −1.043 Z(t-1) −0.9407 Z(t-4) +1.002 Z(t-5) |

Bushehr | 240.062 | MA (5) Hannan-Rissanen | Θ2 and θ3 | X(t) = Z(t) −0.8941 Z(t-1) −0.9832 Z(t-4) +0.9123 Z(t-5) |

Bandar Abbas | 385.08 | MA (13) Innovations | Θ2, Θ3, Θ6, Θ7, Θ8, Θ9, Θ10, Θ11 and Θ12 | X(t) = Z(t) −1.364 Z(t-1) −1.038 Z(t-4) + 1.480 Z(t-5) −0.1564 Z(t-13) |

Hamedan | 309.058 | AR(16) Burg | – | X(t) = −0.9065 X(t-1) −0.6757 X(t-2) - 0.7055 X(t-3) −1.399 X(t-4) −1.236 X(t-5) - 0.9266 X(t-6) −1.031 X(t-7) −1.395 X(t-8) −1.011 X(t-9) −0.7304 X(t-10) −0.8228 X(t-11) −0.8379 X(t-12) − 0.4322 X(t-13) −0.2229 X(t-14) −0.2219 X(t-15) −0.2780 X(t-16) + Z(t) |

Ramsar | 489.954 | AR(12) Burg | – | X(t) = −0.8678 X(t-1) −0.6830 X(t-2) −0.4751 X(t-3) −1.081 X(t-4) −0.9604 X(t-5) −0.7525 X(t-6) −0.4821 X(t-7) −0.7516 X(t-8) −0.5709 X(t-9) −0.4051 X(t-10) −0.2051 X(t-11) −0.2698 X(t-12) +Z(t) |

Yazd | 348.922 | AR(14) Burg | – | X(t) = −0.8125 X(t-1) −0.5461 X(t-2) −0.4024 X(t-3) −1.230 X(t-4) −0.9642 X(t-5) −0.7978 X(t-6) −0.4816 X(t-7) −0.9558 X(t-8) −0.7731 X(t-9) −0.6326 X(t-10) −0.2937 X(t-11) −0.5050 X(t-12) −0.4143 X(t-13) −0.2123 X(t-4) +Z(t) |

Station | AICC Index | Method of time series model | Non-significant θ or φ | Fitted model |
---|---|---|---|---|

Abadan | 302.88 | MA (5) Innovations | Θ2 and θ3 | X(t) = Z(t) −1.043 Z(t-1) −0.9407 Z(t-4) +1.002 Z(t-5) |

Bushehr | 240.062 | MA (5) Hannan-Rissanen | Θ2 and θ3 | X(t) = Z(t) −0.8941 Z(t-1) −0.9832 Z(t-4) +0.9123 Z(t-5) |

Bandar Abbas | 385.08 | MA (13) Innovations | Θ2, Θ3, Θ6, Θ7, Θ8, Θ9, Θ10, Θ11 and Θ12 | X(t) = Z(t) −1.364 Z(t-1) −1.038 Z(t-4) + 1.480 Z(t-5) −0.1564 Z(t-13) |

Hamedan | 309.058 | AR(16) Burg | – | X(t) = −0.9065 X(t-1) −0.6757 X(t-2) - 0.7055 X(t-3) −1.399 X(t-4) −1.236 X(t-5) - 0.9266 X(t-6) −1.031 X(t-7) −1.395 X(t-8) −1.011 X(t-9) −0.7304 X(t-10) −0.8228 X(t-11) −0.8379 X(t-12) − 0.4322 X(t-13) −0.2229 X(t-14) −0.2219 X(t-15) −0.2780 X(t-16) + Z(t) |

Ramsar | 489.954 | AR(12) Burg | – | X(t) = −0.8678 X(t-1) −0.6830 X(t-2) −0.4751 X(t-3) −1.081 X(t-4) −0.9604 X(t-5) −0.7525 X(t-6) −0.4821 X(t-7) −0.7516 X(t-8) −0.5709 X(t-9) −0.4051 X(t-10) −0.2051 X(t-11) −0.2698 X(t-12) +Z(t) |

Yazd | 348.922 | AR(14) Burg | – | X(t) = −0.8125 X(t-1) −0.5461 X(t-2) −0.4024 X(t-3) −1.230 X(t-4) −0.9642 X(t-5) −0.7978 X(t-6) −0.4816 X(t-7) −0.9558 X(t-8) −0.7731 X(t-9) −0.6326 X(t-10) −0.2937 X(t-11) −0.5050 X(t-12) −0.4143 X(t-13) −0.2123 X(t-4) +Z(t) |

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.

### Time series model validation

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.

Station | R^{a} | Station | R^{a} | Station | R^{a} |
---|---|---|---|---|---|

Abadan | 0.713 | Gorgan | 0.623 | Sanandaj | 0.883 |

Ahvaz | 0.856 | Hamedan | 0.962 | Semnan | 0.790 |

Arak | 0.858 | Iran Shahr | 0.785 | Sharekord | 0.906 |

Babolsar | 0.777 | Kashan | 0.868 | Shiraz | 0.874 |

Bandar Abbas | 0.828 | Kerman | 0.845 | Tabas | 0.892 |

Bandar Lenge | 0.812 | Kermanshah | 0.824 | Tabriz | 0.774 |

Birjand | 0.913 | Khoram Abad | 0.812 | Tehran | 0.901 |

Bushehr | 0.868 | Khoy | 0.673 | Torbate Heydarieh | 0.920 |

Chabahar | 0.481 | Mashhad | 0.857 | Yazd | 0.817 |

Dezful | 0.756 | Oroomieh | 0.801 | Zabol | 0.824 |

Esfahan | 0.910 | Ramsar | 0.809 | Zahedan | 0.664 |

Fasa | 0.869 | Rasht | 0.857 | Zanjan | 0.860^{a} |

Ghazvin | 0.851 | Sabzevar | 0.885 |

Station | R^{a} | Station | R^{a} | Station | R^{a} |
---|---|---|---|---|---|

Abadan | 0.713 | Gorgan | 0.623 | Sanandaj | 0.883 |

Ahvaz | 0.856 | Hamedan | 0.962 | Semnan | 0.790 |

Arak | 0.858 | Iran Shahr | 0.785 | Sharekord | 0.906 |

Babolsar | 0.777 | Kashan | 0.868 | Shiraz | 0.874 |

Bandar Abbas | 0.828 | Kerman | 0.845 | Tabas | 0.892 |

Bandar Lenge | 0.812 | Kermanshah | 0.824 | Tabriz | 0.774 |

Birjand | 0.913 | Khoram Abad | 0.812 | Tehran | 0.901 |

Bushehr | 0.868 | Khoy | 0.673 | Torbate Heydarieh | 0.920 |

Chabahar | 0.481 | Mashhad | 0.857 | Yazd | 0.817 |

Dezful | 0.756 | Oroomieh | 0.801 | Zabol | 0.824 |

Esfahan | 0.910 | Ramsar | 0.809 | Zahedan | 0.664 |

Fasa | 0.869 | Rasht | 0.857 | Zanjan | 0.860^{a} |

Ghazvin | 0.851 | Sabzevar | 0.885 |

^{a}Correlation coefficient is significant at 5% level.

### Model prediction

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.

Season | Station | Abadan | …. | Tehran | Yazd | Zabol | Zahedan | Zanjan |
---|---|---|---|---|---|---|---|---|

year | Predicted seasonal SPI index | |||||||

Winter | 2015 | 0.368 | … | 0.522 | 0.870 | 0.513 | 0.307 | 0.498 |

2016 | 0.008 | … | 0.580 | −0.118 | −0.046 | 0.448 | 0.443 | |

2017 | −0.422 | … | 0.463 | −0.618 | −0.502 | −0.493 | −0.897 | |

2018 | 0.714 | … | 0.525 | 0.171 | 0.373 | 0.102 | 0.661 | |

2019 | 0.348 | … | 0.534 | 0.958 | 0.906 | 0.785 | 0.407 | |

Spring | 2015 | −0.159 | … | −0.102 | −0.046 | 0.084 | 0.427 | 0.257 |

2016 | −0.410 | … | −0.188 | −0.600 | −0.290 | −0.257 | −0.960 | |

2017 | 0.506 | … | −0.206 | 0.386 | 0.174 | −0.299 | 0.698 | |

2018 | 0.518 | … | −0.152 | 0.836 | 0.646 | 0.631 | 0.400 | |

2019 | −0.140 | … | −0.154 | −0.013 | 0.029 | 0.224 | 0.282 | |

Summer | 2015 | −0.368 | … | −0.102 | −0.537 | −0.365 | −0.435 | −0.769 |

2016 | 0.467 | … | −0.884 | 0.303 | 0.301 | −0.184 | 0.569 | |

2017 | 0.422 | … | −0.775 | 0.814 | 0.704 | 0.450 | 0.350 | |

2018 | −0.085 | … | −0.620 | −0.016 | 0.125 | 0.420 | 0.266 | |

2019 | −0.247 | … | −0.567 | −0.505 | −0.224 | −0.247 | −0.718 | |

Autumn | 2015 | 0.506 | … | 0.386 | 0.318 | 0.256 | −0.044 | 0.556 |

2016 | 0.427 | … | 0.353 | 0.757 | 0.627 | 0.573 | 0.308 | |

2017 | −0.139 | … | 0.309 | −0.027 | 0.087 | 0.332 | 0.221 | |

2018 | −0.231 | … | 0.315 | −0.452 | −0.244 | −0.202 | −0.670 | |

2019 | 0.452 | … | 0.307 | 0.292 | 0.275 | −0.155 | 0.529 |

Season | Station | Abadan | …. | Tehran | Yazd | Zabol | Zahedan | Zanjan |
---|---|---|---|---|---|---|---|---|

year | Predicted seasonal SPI index | |||||||

Winter | 2015 | 0.368 | … | 0.522 | 0.870 | 0.513 | 0.307 | 0.498 |

2016 | 0.008 | … | 0.580 | −0.118 | −0.046 | 0.448 | 0.443 | |

2017 | −0.422 | … | 0.463 | −0.618 | −0.502 | −0.493 | −0.897 | |

2018 | 0.714 | … | 0.525 | 0.171 | 0.373 | 0.102 | 0.661 | |

2019 | 0.348 | … | 0.534 | 0.958 | 0.906 | 0.785 | 0.407 | |

Spring | 2015 | −0.159 | … | −0.102 | −0.046 | 0.084 | 0.427 | 0.257 |

2016 | −0.410 | … | −0.188 | −0.600 | −0.290 | −0.257 | −0.960 | |

2017 | 0.506 | … | −0.206 | 0.386 | 0.174 | −0.299 | 0.698 | |

2018 | 0.518 | … | −0.152 | 0.836 | 0.646 | 0.631 | 0.400 | |

2019 | −0.140 | … | −0.154 | −0.013 | 0.029 | 0.224 | 0.282 | |

Summer | 2015 | −0.368 | … | −0.102 | −0.537 | −0.365 | −0.435 | −0.769 |

2016 | 0.467 | … | −0.884 | 0.303 | 0.301 | −0.184 | 0.569 | |

2017 | 0.422 | … | −0.775 | 0.814 | 0.704 | 0.450 | 0.350 | |

2018 | −0.085 | … | −0.620 | −0.016 | 0.125 | 0.420 | 0.266 | |

2019 | −0.247 | … | −0.567 | −0.505 | −0.224 | −0.247 | −0.718 | |

Autumn | 2015 | 0.506 | … | 0.386 | 0.318 | 0.256 | −0.044 | 0.556 |

2016 | 0.427 | … | 0.353 | 0.757 | 0.627 | 0.573 | 0.308 | |

2017 | −0.139 | … | 0.309 | −0.027 | 0.087 | 0.332 | 0.221 | |

2018 | −0.231 | … | 0.315 | −0.452 | −0.244 | −0.202 | −0.670 | |

2019 | 0.452 | … | 0.307 | 0.292 | 0.275 | −0.155 | 0.529 |

Drought class | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|

Abadan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Arak | 0 | 0 | 0 | 95 | 5 | 0 | 0 |

Bandar abbas | 0 | 0 | 5 | 95 | 0 | 0 | 0 |

Anzali | 0 | 0 | 20 | 60 | 20 | 0 | 0 |

Esfahan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Fasa | 0 | 0 | 20 | 65 | 10 | 5 | 0 |

Qazvin | 0 | 0 | 0 | 85 | 15 | 0 | 0 |

Hamedan | 0 | 0 | 0 | 80 | 20 | 0 | 0 |

Kerman | 0 | 0 | 5 | 95 | 0 | 0 | 0 |

Kermanshah | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Orumieh | 0 | 0 | 0 | 80 | 20 | 0 | 0 |

Ramsar | 0 | 0 | 5 | 85 | 10 | 0 | 0 |

Rasht | 0 | 5 | 10 | 65 | 20 | 0 | 0 |

Semnan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Shahrekord | 0 | 0 | 15 | 85 | 0 | 0 | 0 |

Zanjan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Drought class | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|

Abadan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Arak | 0 | 0 | 0 | 95 | 5 | 0 | 0 |

Bandar abbas | 0 | 0 | 5 | 95 | 0 | 0 | 0 |

Anzali | 0 | 0 | 20 | 60 | 20 | 0 | 0 |

Esfahan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Fasa | 0 | 0 | 20 | 65 | 10 | 5 | 0 |

Qazvin | 0 | 0 | 0 | 85 | 15 | 0 | 0 |

Hamedan | 0 | 0 | 0 | 80 | 20 | 0 | 0 |

Kerman | 0 | 0 | 5 | 95 | 0 | 0 | 0 |

Kermanshah | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Orumieh | 0 | 0 | 0 | 80 | 20 | 0 | 0 |

Ramsar | 0 | 0 | 5 | 85 | 10 | 0 | 0 |

Rasht | 0 | 5 | 10 | 65 | 20 | 0 | 0 |

Semnan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

Shahrekord | 0 | 0 | 15 | 85 | 0 | 0 | 0 |

Zanjan | 0 | 0 | 0 | 100 | 0 | 0 | 0 |

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

Station | Ljung p-value | Station | Ljung p-value | Station | Ljung p-value |
---|---|---|---|---|---|

Abadan | 0.702 | Gorgan | 0.133 | Sanandaj | 0.033 |

Ahvaz | 0.026 | Hamedan | 0.090 | Semnan | 0.743 |

Arak | 0.571 | Iran Shahr | 0.162 | Sharekord | 0.192 |

Babolsar | 0.065 | Kashan | 0.937 | Shiraz | 0.193 |

Bandar Abbas | 0.237 | Kerman | 0.672 | Tabas | 0.199 |

Bandar Lenge | 0.517 | Kermanshah | 0.326 | Tabriz | 0.155 |

Birjand | 0.055 | Khoram Abad | 0.874 | Tehran | 0.756 |

Bushehr | 0.041 | Khoy | 0.166 | Torbate eydarieh | 0.042 |

Chabahar | 0.096 | Mashhad | 0.020 | Yazd | 0.653 |

Dezful | 0.129 | Oroomieh | 0.533 | Zabol | 0.273 |

Esfahan | 0.012 | Ramsar | 0.201 | Zahedan | 0.011 |

Fasa | 0.019 | Rasht | 0.037 | Zanjan | 0.295 |

Ghazvin | 0.609 | Sabzevar | 0.147 |

Station | Ljung p-value | Station | Ljung p-value | Station | Ljung p-value |
---|---|---|---|---|---|

Abadan | 0.702 | Gorgan | 0.133 | Sanandaj | 0.033 |

Ahvaz | 0.026 | Hamedan | 0.090 | Semnan | 0.743 |

Arak | 0.571 | Iran Shahr | 0.162 | Sharekord | 0.192 |

Babolsar | 0.065 | Kashan | 0.937 | Shiraz | 0.193 |

Bandar Abbas | 0.237 | Kerman | 0.672 | Tabas | 0.199 |

Bandar Lenge | 0.517 | Kermanshah | 0.326 | Tabriz | 0.155 |

Birjand | 0.055 | Khoram Abad | 0.874 | Tehran | 0.756 |

Bushehr | 0.041 | Khoy | 0.166 | Torbate eydarieh | 0.042 |

Chabahar | 0.096 | Mashhad | 0.020 | Yazd | 0.653 |

Dezful | 0.129 | Oroomieh | 0.533 | Zabol | 0.273 |

Esfahan | 0.012 | Ramsar | 0.201 | Zahedan | 0.011 |

Fasa | 0.019 | Rasht | 0.037 | Zanjan | 0.295 |

Ghazvin | 0.609 | Sabzevar | 0.147 |

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.

### Trend analysis of seasonal SPI data series

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

Station | Observed | Observed and predicted | Station | Observed | Observed and predicted |
---|---|---|---|---|---|

Abadan | −0.061 | −0.044 | Khoramabad | −0.111 | −0.109 |

Ahvaz | −0.099 | −0.068 | Khoy | −0.094 | −0.081 |

Arak | −0.080 | −0.087 | Mashhad | −0.054 | −0.056 |

Babolsar | −0.002 | −0.002 | Orumieh | −0.121 | −0.119 |

Bandar Abbas | −0.076 | −0.048 | Ramsar | 0.022 | 0.028 |

Lenge | −0.104 | −0.051 | Rasht | −0.065 | −0.050 |

Birjand | −0.083 | −0.063 | Sabzavar | −0.046 | −0.061 |

Bushehr | −0.014 | −0.003 | Sanandaj | −0.169^{a} | −0.172^{a} |

Chabahar | 0.002 | 0.070 | Semnan | −0.004 | −0.019 |

Dezful | −0.046 | −0.030 | Shahrekord | −0.014 | −0.010 |

Esfahan | 0.035 | 0.050 | Shiraz | −0.031 | −0.023 |

Fasa | −0.065 | −0.054 | Tabas | 0.100 | 0.120 |

Qazvin | 0.000 | −0.049 | Tabriz | −0.155^{a} | −0.135^{a} |

Gorgan | −0.232^{a} | −0.231^{a} | Tehran | −0.003 | −0.036 |

Hamedan | −0.063 | −0.113 | Torbat | −0.067 | −0.073 |

Iranshahr | −0.082 | −0.061 | Yazd | −0.069 | −0.042 |

Kashan | 0.000 | 0.004 | Zabol | −0.076 | −0.034 |

Kerman | −0.081 | −0.063 | Zahedan | −0.066 | −0.039 |

Kermanshah | −0.114 | −0.104 | Zanjan | −0.063 | −0.058 |

Station | Observed | Observed and predicted | Station | Observed | Observed and predicted |
---|---|---|---|---|---|

Abadan | −0.061 | −0.044 | Khoramabad | −0.111 | −0.109 |

Ahvaz | −0.099 | −0.068 | Khoy | −0.094 | −0.081 |

Arak | −0.080 | −0.087 | Mashhad | −0.054 | −0.056 |

Babolsar | −0.002 | −0.002 | Orumieh | −0.121 | −0.119 |

Bandar Abbas | −0.076 | −0.048 | Ramsar | 0.022 | 0.028 |

Lenge | −0.104 | −0.051 | Rasht | −0.065 | −0.050 |

Birjand | −0.083 | −0.063 | Sabzavar | −0.046 | −0.061 |

Bushehr | −0.014 | −0.003 | Sanandaj | −0.169^{a} | −0.172^{a} |

Chabahar | 0.002 | 0.070 | Semnan | −0.004 | −0.019 |

Dezful | −0.046 | −0.030 | Shahrekord | −0.014 | −0.010 |

Esfahan | 0.035 | 0.050 | Shiraz | −0.031 | −0.023 |

Fasa | −0.065 | −0.054 | Tabas | 0.100 | 0.120 |

Qazvin | 0.000 | −0.049 | Tabriz | −0.155^{a} | −0.135^{a} |

Gorgan | −0.232^{a} | −0.231^{a} | Tehran | −0.003 | −0.036 |

Hamedan | −0.063 | −0.113 | Torbat | −0.067 | −0.073 |

Iranshahr | −0.082 | −0.061 | Yazd | −0.069 | −0.042 |

Kashan | 0.000 | 0.004 | Zabol | −0.076 | −0.034 |

Kerman | −0.081 | −0.063 | Zahedan | −0.066 | −0.039 |

Kermanshah | −0.114 | −0.104 | Zanjan | −0.063 | −0.058 |

^{a}Correlation coefficient of trends is statistically significant at 5% level.

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.

## CONCLUSIONS

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.