Seasonal drought forecasting in arid regions, using different time series models and RDI index

Early indications of possible drought can help to plan drought mitigation strategies and measures. Modeling and predicting the behavior of the phenomenon in a standard time period is the best way to address damaging this natural phenomenon. Drought damages environmental factors and agriculture, vegetation, humans and wildlife, as well as local economies (Azarakhshi et al. 2011; Dastorani & Afkhami 2011; Nohegar et al. 2013; Zarei et al. 2013). This phenomenon begins naturally and ends slowly and its damage will be very serious and long-lasting, especially in arid and semi-arid regions. It can result in long-term harmful effects on water resources and the environment (Zarei et al. 2016; Zarei 2018).

In recent decades, many techniques have been used as suitable tools for modeling and forecasting meteorological information such as drought (Soltani et al. 2007; Shamshirband et al. 2015). Chun et al. (2013) evaluated and predicted the impact of climate change on drought in the UK using ARIMA models and the generalized linear model (GLM) approach. Results indicated that the drought pattern in the 2080s is less certain than for the 1961–1990 period, based on the Shannon entropy, but droughts are expected to be more clustered and intermittent. Jahandideh & Shirvani (2011) used time series models in Fars Province to forecast drought based on the standardized precipitation index. According to the results of this research, the SARIMA model with the minimum of Akaike's information criterion bias-corrected (AICC) was selected as the best model.

Mossad & Ali Alazba (2015) developed several ARIMA models for drought forecasting in a hyper-arid climate using the SPEI index. The results reveal that all developed ARIMA models demonstrate the potential ability to forecast drought over different time scales. This study recommends that ARIMA models can be very useful tools for drought forecasting. Djerbouai & Souag-Gamane (2016) investigated the ANN models for drought forecasting in the Algerois basin in Algeria in comparison with traditional stochastic models (ARIMA and SARIMA models). Results of this study showed that despite the linear property of SARIMA models, they yielded satisfactory performances, with respect to various model efficiencies, for SPI-12 forecasts at one-month lead time. Paul et al. (2017) used time series models (robust statistics) to assess the trend of rainfall data at Rajahmundry city located in lower Godavari basin, India. Phuong et al. (2018) evaluated the spatiotemporal variability of annual and seasonal rainfall time series in Ho Chi Minh City for the period 1980–2016. The results of trend estimation also indicated higher increasing rates of rainfall in the dry season compared to the rainy season at most stations. Other researchers have also used time series models to assess drought characteristics and conditions of climate parameters (Barua et al. 2010; Hong-yan et al. 2015; Wang et al. 2016; Zarei & Moghimi 2017; Bahrami et al. 2018; Zarei 2018).

Proper management of drought requires more precise monitoring, modeling and forecasting of this phenomenon using strong drought indices. In this research, the reconnaissance drought index (RDI) was used due to a number of advantages over other drought indices and because it can be used to quantify most types of drought events. This index was proposed by Tsakiris & Vangelis (2005) and for calculation of this index the ratio of cumulative values of precipitation to potential evapotranspiration was used. In addition, this index has certain advantages compared to indices based on precipitation, because it is more representative of the deficient water balance conditions. Since RDI resolved more climatic parameters, such as evapotranspiration, which had an important role in water resource losses in the Iranian basins, it was worthwhile to consider RDI in drought monitoring in Iran (Asadi Zarch et al. 2011). Therefore, for the regions with no sufficient data to calculate other indicators, the RDI index is the best and most comprehensive index.

The aims of this research are modeling and predicting of seasonal drought (calculated based on the RDI index, as a strong index for monitoring drought) in southern Iran using different time series models (ARIMA models). There is an integrated system from the model establishment, verification, and forecasting.

Figure 1 shows the study area boundaries enclosing a 594,996.54 km² area approximately between 25°17′N and 31°11′N latitudes and between 50°49′E and 62°20′E longitudes, in addition to the distribution of 16 synoptic stations used in this study area. These stations have a good distribution with sufficient length of meteorological data (31 years), which can spatially and temporally support drought modeling and forecasting studies. The central and northern areas are highlands and mountains, while the southern and western areas are mainly flat. The elevation of selected stations varied from 5 (at Jask) to 2,030 m (at Abadeh) from free sea surface level. The average annual precipitation varied from 53.1 (at Zabol) to 330.6 mm (at Shiraz) and the average annual potential evapotranspiration varied from 1,416.2 (at Jask) to 2,839.7 mm (at Bushehr). More detailed characteristics of the 16 surveyed stations are presented in Table 1.

Since the RDI considers the proportion of precipitation to PET, at first, precipitation data and the data that was required to compute PET (maximum and minimum temperature, maximum and minimum relative humidity, wind speed, and sunshine) were gathered, and the missing data was substituted with the corresponding long-term mean (Dinpashoh et al. 2011). ET₀ was estimated using the PM-56 equation (Allen et al. 1998) at the monthly time scale for each synoptic station.

The reconnaissance drought identification and assessment index was used in this research (RDI) because of the possible role of ET₀ in the detection of drought events. This index was proposed by Tsakiris & Vangelis (2005), utilizing the ratios of precipitation over Reference Crop Evapotranspiration (ET₀) for different time scales to be representative of the desired region. Calculating RDI was carried out using equations as follows:

The ARIMA model (the integrated ARMA) is a broadening of the class of ARMA that includes differencing (an important technique in data transformation; it attempts to de-trend to control autocorrelation and achieve stationary time series).

ARIMA modeling generally involves three stages as follows. First stage: model identification by specifying the type of the model (AR, MA, ARMA, or ARIMA) and its order. This identification is sometimes undertaken by looking at plots of the sample autocorrelation function (ACF) and sample partial autocorrelation function (PACF) and sometimes it is employed by an autofit procedure fitting many different possible model structures and orders and using a goodness-of-fit statistic to select the best model. Second stage: estimate the coefficients of the model by minimizing the sum of squared residuals. Third stage: model diagnostics. In this stage, it is very important to check that the residuals of the candidate model are random and normally distributed and the estimated parameters are statistically significant. It should be noted that the best model of all which fit the data is the one which has the fewest parameters.

In the cross-validation procedure, the data set split into two sets: training sample and prediction. The training sample is used to develop a model for prediction and the prediction set is used to evaluate the rationality and predictive ability of the selected model. This validation procedure is the statistical practice of splitting a sample of data into two subsets so that the analysis is initially performed on one subset and the other subset is retained for subsequent use in confirming and validating the initial analysis. For fitting the ARMA model, it must be at least likely that the data are in fact a realization of an ARMA process and, in particular, a realization of a stationary process. In the stationary time series, the statistical properties such as mean, variance, autocorrelation and so on are all constant over time. In order to obtain a stationary time series, a sequence of mathematical transformations (Box–Cox transformation, mean subtraction, and the differencing) was used.

When the variability of the data increases or decreases with the level, this transformation is useful. The variability can be made nearly constant by a suitable choice of ⁠. For instance, the variability of a set of positive data whose standard deviation increases linearly can be stabilized by choosing (Brockwell & Davis 2002).

The lowest order of differencing is the correct amount of differencing that resulted in time series which fluctuate around a well-defined mean value and whose ACF plot decays rapidly to zero, either from above or below. Thus, at every stage of differencing, the plots of the sample ACF and the sample PACF were checked to see where the ACF/PACF were out of the bounds ±1.96/⁠. The stationary series is the series with a sample ACF that decays fairly rapidly. If the ACF plot has a polynomial trend, it shows that the series still has some trends. The periodicity of ACF shows that the series has seasonality and some more differencing for the data should be applied.

The residual ACF/PACF of the models and the randomness of the residuals should be checked. For the sample with large n, the sample autocorrelations of an independent and identically distributed (iid) sequence (⁠⁠) with zero mean and finite variance are approximately iid with normal distribution N (0,1/n) (Amei et al. 2012). The consistency of the observed residuals with iid noise was considered by examining the sample correlations of the residuals and rejecting the iid noise hypothesis if more than two or three out of 40 falls outside the bounds ±1.96/ or if one falls far outside the bounds (Brockwell & Davis 2002).

The Ljung–Box test (Ljung & Box 1978) was used to check whether the residuals of a fitted model are iid in ARIMA modeling or not. This test was based on the autocorrelation plot and tests the overall independence based on a few of the time lags. The Ljung–Box test is defined as follows:

• H₀: The sequence data are iid

• H_a: The sequence data are not iid

The test statistic is ⁠, where is the estimated autocorrelation at lag-k and is equal to ⁠, n is sample size, m is the number of lags being tested and are the residuals after a model has been fitted to a series Z₁, …, Z_n. For large sample sizes (n), the distribution of is approximately under the null hypothesis, where p + q is the number of parameters of the fitted model. The hypothesis of iid is rejected if at level and the sequence data do have autocorrelations significantly different from zero and a new search for a fitted ARMA model for a mean-corrected data set should be followed.

The best ARIMA models that fitted to the data were used to forecast future values of the time series from the observed values.

After computation of three-month RDI for 16 synoptic stations of the study area (southern Iran), ARIMA models were fitted to these data and then optimized by eliminating non-significant coefficients. Then, using various indicators, the goodness of fitted models was evaluated. For this purpose, heterogeneity and randomness of the residuals, model validity in the forecast, and comparison between ACF/PACF of data and the fitted model were considered.

where RDI_i is computed RDI for the subset of observed data (2008–2010), is the arithmetic mean of RDI_i, is the forecasted RDI for the subset of observed data (2008–2010), and is the arithmetic mean of ⁠.

In general, high values of NSE (up to 100%) and correlation coefficient (R-squared value) and small values for RMSE and MAE indicate a good model.

The results of drought class probabilities relative to the three-month RDI for 16 synoptic stations of this study region for observed data (time period of 1980–2010) are shown in Table 4. Generally, the drought classes varied between ‘very wet’ to ‘severely dry’, and the number of non-zero probabilities related to the ‘severely dry’ class is considerably higher than other drought classes that indicated the drought occurrence in this time period at this study area. In the Bam, Bandar Abbas, Bandar Lengeh, Jask, Chabahar and Zabol stations, drought class probabilities varied between ‘very wet’ to ‘moderately dry’, which indicated the lowest severity of drought in these stations. The Abadeh, Shiraz, Fasa, Sirjan, Kerman, Shahre Babak and Saravan stations showed the drought class probabilities between ‘moderately wet’ to ‘severely dry’. Therefore, in the latter stations, the highest severity of drought occurred. In Bushehr and Zahedan stations, drought class probabilities varied between ‘moderately wet’ and ‘moderately dry’ which indicated the normal conditions related to drought.

Initially, in order to obtain stationary time series, Box–Cox transformation (to achieve constant variance), mean subtraction, and the differencing (to eliminate trend (k = 1) and periodicity (d = 4)) were used. Figure 2 shows the computed three-month RDI time series for observed data (1980–2010) for the synoptic meteorological stations of Abadeh, Bushehr and Zabol (results of other stations are not shown due to space). Figure 3 shows the stationary time series for the mentioned stations.

Determination of the order of p and q in the AR, MA, and ARMA models was examined by using the ACF/PACF plots (Figure 4). According to Figure 4, up to 5% of residual ACF/PACF/s are out of zero range (the dotted lines) and the sample and model ACF/PACF are close together indicating the appropriateness of the model. After selecting the best values of p and q, different models were fitted to the computed three-month RDI time series for observed data (1980–2010) in all of the stations and, according to the AICC statistic, the model with the minimum value of AICC was selected as the best model. The results of model diagnostics and fitting shown in Table 5 indicate the best models with the minimum value of AICC at each station. Results indicated that the best performances in seasonal drought forecasting were as follows: the AR model based on the Yule–Walker method in 6.3% of cases; the AR model based on the Burg method in 12.5% of cases; the MA model based on the Hannan–Rissanen method in 25% of cases; the MA model based on the Innovations method in 37.5% of cases; the ARMA model based on the Hannan–Rissanen method in 25% of cases; and the ARMA model based on the Innovations method in 12.5% of cases.

The AR coefficients (Ø_i) are presented in Table 6.

We checked the residual ACF/PACF of the models and the randomness of the residuals. According to the P-values obtained for the Ljung–Box statistic in various lags (Table 5), homogeneity and randomness of residuals (p > 0.05) and the reliability of forecasts at the level of 95% can be deduced.

To validate the models with minimum AICC statistic that were obtained for different stations (Table 5), using the subset of observed data (1980–2007), the values for the 2008–2010 time period (approximately 10% of the time period used for modeling) were forecast and compared with values observed in this time period. The results of this comparison are presented in Table 8. According to the Pearson test, a significant correlation between observed and forecasted values is noticed at the significance level of 0.01. This comparison is also shown in Figure 5, including all of the synoptic meteorological stations of this study area. Also, Table 9 shows the results of NSE, RMSE and MAE measures to evaluate the goodness of fitted models for forecasting. High values of NSE and correlation coefficient (R-squared value) and small values for RMSE and MAE indicate a good model.

Drought predicting plays an important role in the planning and management of water resource systems by significantly reducing drought impacts. In this study, different statistical models were evaluated to select the best time series models for forecasting droughts based on the RDI index in southern Iran. For comparison and evaluation of forecasting models, various performance measures (NSE, RMSE, MAE and R²) were used. Drought classes of observed data (1980–2010) varied between ‘very wet’ to ‘severely dry’. The occurrence frequency of the ‘severely dry’ class was considerably higher compared to other drought classes. This result indicates the drought occurrence in this time period in this study area. Drought class probabilities indicated the lowest severity of drought in Bam, Bandar Abbas, Bandar Lengeh, Jask, Chabahar and Zabol stations and the highest severity of drought in Abadeh, Shiraz, Fasa, Sirjan, Kerman, Shahre Babak and Saravan stations.

Between different models (AR, MA, ARMA and ARIMA) fitted to the computed three-month RDI time series for observed data (1980–2010) in all of the stations, the model with the minimum value of AICC was selected as the best model. Results indicated that in 18.8% of cases for the AR models, in 50.0% of cases for the MA models and in 31.3% of cases, the ARMA models had the best performance in seasonal drought forecasting. For evaluating the forecast performance of all fitted models, different measures of goodness of fit (NSE, RMSE and MAE) were computed. The results verified the goodness of fitted models.

The results of this study are applicable for water resources managers, decision makers and governments to face this phenomenon well prepared. They can achieve this by proper management of water resources consumption, especially in the agricultural sector, which is the main consumer of water resources in Iran, in order to control or prevent the effects of drought.