Development of a Non-stationary Standardized Precipitation Evapotranspiration Index (NSPEI) for Drought Monitoring in a Changing Climate

In a changing climate, drought indices as well as drought definitions need to be revisited because some statistical properties, such as the long-term mean, of climate series may change over time. This study aims to develop a Non-stationary Standardized Precipitation Evapotranspiration Index (NSPEI) for reliable and robust quantification of drought characteristics in a changing climate. The proposed indicator is based on a non-stationary log-logistic probability distribution, assuming the location parameter of the distribution is a multivariable function of time and climate indices, as covariates. The optimal non-stationary model was obtained using a forward selection method in the framework of the Generalized Additive Models in Location, Scale, and Shape (GAMLSS) algorithm. The Non-stationary and Stationary forms of SPEI (i.e., NSPEI and SSPEI) were calculated using the monthly precipitation and temperature data of 32 weather stations in Iran for the common period of 1964–2014. The results showed that almost at all the stations studied, the non-stationary log-logistic distributions outperformed the stationary ones. The AICs of the non-stationary models for 97% of the stations were lower than those of the stationary models. The non-stationary models at 90% of the stations were statistically significant at the 5% significance level. While SSPEI identified the long-term and continuous drought and wet events, NSPEI revealed the short-term and frequent drought/wet periods at almost all the stations of interest. Finally, it was revealed that NSPEI, compared to SSPEI, was a more reliable and robust indicator of drought duration and drought termination in vegetation cover during the severest drought period (the 2008 drought). Therefore, it was suggested as a suitable drought index to quantify drought impacts on vegetation cover in Iran.


Introduction
Drought, as a recurrent and normal feature of the climate, affects the economy, society, and environment across the globe (Wilhite 2000). Historical evidence shows that Iran has suffered from long-term and devastating droughts during few centuries ago. For instance, the intensive famines during 1870-1872 and 1917-1919, which occurred due to drought occurrences, endangered water and food security throughout the country, and destroyed half of the population (De Planhol 2012). Bazrafshan et al. (2017) specified that during 1894-2010, Iran (with an average annual precipitation of 254 mm) had experienced 23 drought events, ranging from 1 to 10 years. The and 1998 droughts (with the total precipitation deficits of 176.1 and 180.4 mm, respectively) were identified as the most severe and extensive drought events.
Drought is conventionally defined as a temporary period of moisture deficit relative to the long-term average at a given location (Mishra and Singh 2010;Paulo and Pereira 2006;Van Loon 2015). This simple definition, however, cannot be used in the regional studies of drought. To solve this issue, drought indices are used. Using standardization of moisture deficiencies, drought indices enable us to compare different locations based on their drought characteristics (such as severity and duration) (Sen 2015;Tsakiris et al. 2007). Several drought indices have been presented to monitor meteorological, agricultural, and hydrological droughts (Wilhite and Glantz 1985). Some of most popular traditional drought indices are: Palmer Drought Severity Drought Index (PDSI) (Palmer 1965), Standardized Precipitation Index (SPI) (McKee et al. 1993), Reconnaissance Drought Index (RDI) (Tsakiris and Vangelis 2005), Standardized Precipitation Evapotranspiration index (SPEI) (Vicente-Serrano et al. 2010). While the SPI is only based on precipitation data, the other above-mentioned indices principally employ precipitation and evapotranspiration data.
Incorporating environmental changes in the context of drought indices is an important issue which is not regarded in the traditional drought indices . So far, efforts have been made to develop drought indices under climate non-stationarity. According to the literature, the first study in this context was carried out by Russo et al. (2013). They developed a non-stationary Gamma distribution in the SPI structure and called the new index the Standardized non-stationary Precipitation Index (SnsPI). Wang et al. (2015) proposed a time-dependent Standardized Precipitation Index (SPIt) in the Luanhe River Basin, China. Non-stationarity in SPIt was modeled by fitting a Gamma distribution with a time-dependent location parameter to the summer precipitation data. Non-stationary models were introduced as polynomial regression functions of time within the Generalized Additive Models in Location, Scale, and Shape (GAMLSS) framework. The results showed that in most situations, the non-stationary Gamma distribution has a better fit to the precipitation series compared to its stationary form. Bazrafshan and Hejabi (2018) developed a Non-stationary Reconnaissance Drought Index (NRDI) for drought monitoring under changing climate conditions in Iran. They considered the location parameter of the log-normal probability distribution as a polynomial regression function of time. Results showed that, in most stations, the non-stationary log-normal distribution was able to simulate the monotonic downward trend identified by the Mann-Kendall test, and outperformed the stationary one.
One of the important limitations of time-dependent drought indices such as SnsPI, SPIt, and NRDI is that the assumption of the linear relationship of distribution parameters with time is not always true in the future behavior of climate data, because this linear relationship may be part of a larger cycle that is not manifest in current data . Furthermore, the time-dependent non-stationary models can not accurately represent the variability in distribution parameters (Serinaldi and Kilsby 2015). Climate signals (i.e., teleconnections) such as ENSO and NAO have been proposed as the explanatory variables (or covariates) in some studies Rashid and Beecham 2019;Wang et al. 2015) to address this issue. In addition to the mentioned covariates, the anthropogenic effect is also incorporated in some studies (Kang and Jiang 2019;Wang et al. 2020) to model nonstationarity in hydrological drought behavior. The results showed that incorporating the climate signals and the anthropogenic effects into the framework of traditional drought indices has improved the reliability and robustness of drought characteristics.
Climate change, in the most optimistic case, will change the average of climatic elements such as precipitation and air temperature (IPCC 2014), and this change in the average would not be consistent with the existing definitions of droughts. Therefore, the definition of drought in a changing climate needs to be revisited. The main effects of future climate change are the regional increase or decrease in precipitation and the global increase in air temperature (Bazrafshan 2017;IPCC 2014). Assuming the stationary precipitation series and low variability of other climatic factors such as air temperature, the traditional SPI is capable of representing the effect of precipitation variations on drought (Russo et al. 2013). For non-stationary precipitation data, the development of a Non-stationary SPI (NSPI) for drought assessment has been proposed Rashid and Beecham 2019;Wang et al. 2015). Although the NSPI considers precipitation changes in a changing climate, it, like the SPI, is unable to demonstrate the effect of global warming on drought characteristics (Zarch et al. 2015). In such conditions, the SPEI may be a suitable option. However, SPEI calculations are also reliable under the assumption of climate stationarity. In a changing climate where both precipitation and evapotranspiration or one of them are non-stationary, the drought assessment and the statistical analyses based on the stationary assumption of data would not be reliable. Therefore, in the continuation and progression of the previous studies, there will be a need to develop the non-stationary form of SPEI for more robust and reliable monitoring of meteorological drought in a changing climate.
The aim and novelty of this study is to develop a Non-stationary Standardized Evapotranspiration Precipitation Index (NSPEI) under changing climate conditions in Iran. Nonstationary modeling of climate variables is carried out using the covariates, namely, time and the teleconnections (such as ENSO and NAO) affecting the climate of Iran. Both stationary and non-stationary indices are compared in terms of temporal and spatial drought characteristics. Also, the capability of the proposed drought index for quantifying the impact of the severest drought period (the 2008 drought) on vegetation cover in Iran is compared to the traditional SPEI.

Study Area
Geographically, Iran is located in West Asia in the range of 63°-44° East longitude and 29°-35° North latitude (Fig. 1). Most of Iran's climate is arid (64%) and semi-arid (20%), due to its location in the subtropical high-pressure subsidence zone. In addition, about 16% of its area is Mediterranean to very humid (Khalili 1997). The main controllers of Iran's climate are latitude, altitude, and distance from the seas and oceans. The long-term average annual precipitation in Iran is 254 mm, which varies between 13 and 2003 mm across the country. The long-term average temperature in the country varies between 1.6 • C (at the elevation of 3000 m) and 28 • C (on the southern coast) (Khalili 1997;Khalili and Rahimi 2018). There are a large number of studies that connect large-scale atmospheric-oceanic phenomena (i.e., teleconnections) to the occurrence of droughts and floods in different parts of the country in different seasons of the year Khalili 2006, 2008;Golian et al. 2015;Moradi 2004;Nazemosadat and Cordery 2000;Nazemosadat and Ghasemi 2004). Research in this context suggests that four teleconnections, namely El Niño-Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), Arctic Oscillation (AO), and North Caspian Sea Pattern (NCP), have major effects on climate variability in Iran.

Data
In this study, three different datasets were used: 1) monthly precipitation and air temperature data recorded at 32 weather stations ( Fig. 1) across Iran, 2) large-scale climate indices (teleconnections) including the monthly data of ENSO, NAO, AO, and NCP, and 3) Normalized Difference Vegetation Index (NDVI) data (Pinzon and Tucker 2014). The groundbased station data covers the common period of 1964 to 2014, and the satellite-based data covers the period of 1981 to 2014.

Standardized Precipitation Evapotranspiration Index (SPEI)
Calculation of the traditional SPEI (Vicente-Serrano et al. 2010) includes several steps, as follows: (i) Calculation of a simple climatic water balance according to the following equation:]  Table 1 where D i is the difference between precipitation ( P i ) and atmospheric evaporative demand ( AED i ) in the month number i in the time series. The method of Thornthwaite (1948) was employed to calculate AED in this study according to some previous studies (Mavromatis 2007;Vangelis et al. 2013;Zarei and Mahmoudi 2020).
(ii) Selection of timescale: In this research, SPEI was calculated on a 12-month timescale due to the management of Iran's water resources on a water-year scale (i.e. from October of the current year to September of the next year) (Raziei et al. 2008). (iii) Moving aggregation of consecutive D-values at a given timescale: where x t,k is the cumulative D-values of the k-month timescale (here k = 12 ) at the time t ( t = 1, 2, … , n ) and is called the x-series, hereafter. (iv) Fitting a three-parameter log-logistic distribution function to the x-series as below: where , and are the parameters of scale, shape, and location of the distribution, respectively, and x is the cumulative series of D-values at a given timescale. The parameters of the log-logistic distribution function (Eq. (3)) are estimated using the L-moment method (Hosking 1990). The log-logistic cumulative distribution function of x is defined by the expression: (v) Calculation of SPEI: Using the approximation of Abramowitz and Stegun (1965), the cumulative probability of any value of x (i.e.,F(x) ) is transformed to the standard normal variate (i.e., SPEI), as follows: in which, and P = 1 − F(x) . If P > 0.5 , 1 − P is substituted for P in Eq. (6) and the sign of SPEI is reversed. The constants in the SPEI equation are: Since the traditional SPEI is calculated under stationary condition; henceforth, we will refer to it as the Stationary SPEI (SSPEI) from now on. The negative values of SSPEI demonstrate the dry periods (droughts), and the positive values of SSPEI imply the wet periods. (1)

Non-stationary SPEI (NSPEI)
While the SSPEI assumes all parameters of the log-logistic distribution fitted to the x -series (the cumulative D-values) are constant (stationary), the NSPEI fits an optimal non-stationary model to the series. The non-stationary statistical model assumes that one or more parameters of the log-logistic parameter vary as a function of some external covariates. The external covariates selected in this study are the time variable and the teleconnections (including ENSO, NAO, AO, and NCP) influencing the climate conditions of Iran (Alizadeh-Choobari et al. 2018;Dezfuli et al. 2010;Khalili 2006, 2008). Many researchers believe that climate change, in most probable cases, changes the average (the location parameter) of climate variables on a regional or global scale (Cheng and AghaKouchak 2014; Kwon and Lall 2016;Li et al. 2015;Russo et al. 2013;Wang et al. 2015). Therefore, the non-stationary models used in the study assume that the location parameter of the log-logistic distribution is a multivariable function of the covariates as below: are the climate indices at time t − l in which l is the lag-times 0, 1, … , 12 months. The other parameters of the log-logistic distribution are assumed to be constant. The optimal combination of the variables in Eq. (7) for each month of the year at a 12-month timescale is selected using the forward selection method in the framework of the Generalized Additive Models in Location, Scale, and Shape (GAMLSS) algorithm (Rigby and Stasinopoulos 2005). The GAMLSS is an advanced distributionbased method for (semiparametric) regression analysis in which the response variable is assumed to have a parametric distribution, and the parameters of this distribution can vary depending on covariates (Stasinopoulos et al. 2008).
In this study, GAMLSS is applied to model a non-stationary log-logistic distribution using a (semi) parametric relationship between the location parameter (t) and each combination of the covariates through the monotonic link function ( g(.) ), as below: where is a parameter vector of size J , X is a fixed, known design vector of size J , j is a normal-distributed vector of random effects parameters, and Z j is a non-parametric additive function of j . All parameters are estimated using the maximum (penalized) likelihood estimator (Stasinopoulos et al. 2008).
Equation (8) is evaluated for different combinations of 1 to 5 variables, corresponding to the five covariates used in this study. The best non-stationary model is the one with the least AIC value (AICmin). If the differences between the AIC of the best model and the AICs of some of the remaining ones are less than 2, the optimal model is that with a number of covariates lower than the best one (Burnham and Anderson 2002). The optimum model is also considered for the explained variability in data using the statistical technique of the analysis of variance (i.e., the ANOVA table). The optimum model is statistically significant where the p-value corresponding to the F-statistic in the ANOVA table is less than the given significance level of . It is worth noting that, with the support of the GAMLSS package incorporated into the R software, the computation processes for evaluation of the non-stationary models were accomplished (Stasinopoulos et al. 2008).
In order to calculate NSPEI, the cumulative probability corresponding to any values in the x-series is estimated from the non-stationary log-logistic distribution and then is transformed inversely to the standard normal variate (here called NSPEI). The NSPEI values less than zero imply dry conditions (droughts) and those equal to or more than zero indicate wet conditions.

Vegetation Condition Index
Knowing that precipitation and evapotranspiration directly impact on the growth and development of vegetation, it is expected that a drought index such as SPEI, which is based on the two variables, should be able to detect changes in vegetation caused by meteorological conditions. In this study, using the Vegetation Condition Index (VCI) (Kogan 1995), temporal variations in vegetation over the study area were assessed in association with SSPEI and NSPEI to determine the necessity of incorporating the non-stationary framework into the SPEI calculations. The VCI is an NDVI-based remote sensing drought index which is computed using the following formula: where NDVI is the smoothed monthly NDVI, and NDVI min and NDVI max are the multi-year minimum and maximum NDVI values, respectively, in each grid cell in each month of the year. VCI ranges from 0 to 100%, showing extremely unfavorable to optimal conditions for vegetation. Unlike the NDVI, which represents the compound effect of weather and ecology in drought monitoring, the VCI is able to separate the ecology component from the NDVI data and to quantify only the effect of meteorological droughts on vegetation (Rahimzadeh Bajgiran et al. 2008;Shamsipour et al. 2008).
In this study, the VCI was calculated from the NDVI averaged over a 3 × 3 window, including the nine pixels around the position of each weather station (Rahimzadeh Bajgiran et al. 2008). Then, a 12-month moving average method, like the SSPEI and NSPEI, was applied to the VCI series. A drought event in vegetation happens when the VCI is less than 40% (Ghaleb et al. 2015).

Trend Analysis
A nonparametric Mann-Kendall (MK) statistical test (Kendall 1970;Mann 1945) was employed to evaluate the trend orientation in x-series (Eq. (2)) at all stations of interest. The benefit of this test over others is that it analyzes the presence of a monotonic (linear or nonlinear) trend in the time series. If the existence of trends in x-series is confirmed, it strengthens the case for using non-stationary probability distributions in the context of drought indices. Considering the serially independent data u i , i = 1, 2, ..., n , in which n is the length of data, the MK test statistic Z is defined as: where S is sum of the sign function values according to the following formula: , S < 0 and Var(S) is the variance of S which is defined as: where g is the number of tied groups and t p is the number of observations in the p-th group.
The null hypothesis (no monotonic trend) against the alternative hypothesis (monotonic trend is present) is rejected if |Z| > Z 1− ∕2 in which Z 1− ∕2 follows the standard normal variate at the significance level of . A positive (negative) significant value of Z indicates an upward (downward) monotonic trend. Note that the lag-1 autocorrelations in the data series were removed prior to the data being further analyzed for trend detection according to the approaches presented in Hamed and Rao (1998) and Yue et al. (2002). Figure 2 displays the boxplot of the Spearman correlation coefficients between the 12-month aggregated D-values (i.e. the x-series) and each of the four selected climate indices (ENSO, AO, NCP, and NAO) used in this study for the 0-to 12-month lag-times. As can be seen in the figure, while the x-series has negative correlations with ENSO, its correlations with the other three climate indices are mostly positive. The x-series has significant correlations with ENSO in more than half of the stations for the 0-to 1-month lag-times and, with increasing the lag-times to 12 months, the correlations become nonsignificant. In a number of stations, correlations of x-series with AO, NCP, and NAO are significant for lag-times longer than five months, shorter than eight months, and 0-1 and longer than seven months, respectively. The results showed the relative significance of ENSO, compared to the other three indices, in describing the behavior of the x-series at the stations of interest. Table 1 shows the step-by-step forward selection of the covariates (namely t, ENSO, AO, NCP, and NAO) for estimation of the location parameter of the log-logistic distribution at the stations of interest for the 12-month sub-period of October-September. The selection process for suitable covariates was carried out in five consecutive steps, according to Eq. (7). The variables were chosen for their significance in estimating the location parameter; the most important variable is set in the first position, and the other variables occupy the second to fifth position. In this selection method, five non-stationary models with a combination of 1 to 5 variables were constructed to estimate the location parameter at the studied stations.

Selection of the Optimal Non-stationary Models
For example, as shown in Table 1, at Abadan station, according to the AIC values in steps (1)-(5), the bivariate model with covariates t and NCP(t − 5) has the lowest AIC value and is considered the best non-stationary model to estimate the location parameter. However, the univariate model in step (1), which has only the t-covariate, has been selected as the optimal model because the difference between the AIC of the optimal model and that of the best model is not significant and the number of variables used in the optimal model is less than that of the best model. Also, the selected optimal model has a AIC of 720 which is less than the AIC of the stationary model (743.3). The non-stationary optimal models outperformed the stationary model at almost all stations except the Kerman station. This result is valid for all other sub-periods, in addition to the sub-period of October-September. Table 2 shows the mathematical formulae of the optimal non-stationary models to estimate the location parameter of the log-logistic distribution at the studied stations for the 12-month sub-period of October-September. As mentioned above, a time-dependent model is sufficient to model the non-stationarity of the location parameter at the Abadan station (#1). If time is a strong covariate in the model, it should also be verified by the non-parametric Mann-Kendall (MK) test. As can be seen in Table 2, the MK statistic at this station is significant at the 99% confidence level. The negative slope of t in the optimal model ( t = 2949.7 − 10.9t ) and the negative value of the MK statistic, both of which suggest a clear downward trend in the x -series at the station. As a general rule, it may be concluded that the larger the MK statistic, the stronger the effect of the time covariate on non-stationarity in the location parameter (Bazrafshan and Hejabi 2018).
The five covariates used in this study were compared to the times selected in the optimal non-stationary models across the studied stations. ENSO had the highest frequency, and the frequencies of NCP, t, NAO, and AO were ranked in the next positions, respectively. Notice that the importance of t in the optimal non-stationary models was higher than in NAO and AO. This result confirmed the findings of the other researchers (Golian et al. 2015;Nazemosadat and Cordery 2000;Nazemosadat and Ghasemi 2004) who emphasized the importance of ENSO on the precipitation pattern over Iran.

Table 1
Forward selection of covariates for non-stationary modeling of the location parameter of log-logistic distribution for the 12-month sub-period of October-September.
nAIC indicates the Akaike Information Criterion (AIC) for each step of the selection process of a non-stationary model. sAIC shows the AIC for the stationary assumption of the log-logistic distribution parameters Station Selected variables in Forward Selection nAIC sAIC Step(1) Step (2) Step (3) Step (4) Step (5) Step (1) Step (2) Step (3) Step (4) Step ( Table 1 (continued)

nAIC sAIC
Step(1) Step (2) Step (3) Step (4) Step (5) Step (1) Step (2) Step (3) Step (4) Step ( Table 2 The optimal non-stationary model (along with its corresponding AIC and p-value) fitted to the x-series of October-September for estimation of the location parameter ( (t)) of log-logistic distribution at the stations of interest. The MK test gives the trend analysis of the x-series

Comparison of SSPEI and NSPEI Time Series
The temporal behaviors of SSPEI and NSPEI for the sample stations in the arid (Yazd), semi-arid (Tabriz), Mediterranean (Gorgan), and per-humid (Bandar Anzali) climates were compared, as shown in Fig. 3. According to the figure, there was a slight difference between the time series of SSPEI and NSPEI in the extreme climates (arid and very humid). On the other hand, there was not much coincidence between the two indices in the semi-arid and Mediterranean climates. Especially in the Tabriz station (with a semi-arid climate), this difference was large in the final and early parts of the record period relative to the middle part. The large difference between the results of the two indices at the Tabriz station and, to some extent, at the Gorgan station seems to be attributed to the central role of covariates in non-stationary models. Focusing on the Tabriz station reveals that time covariate was the most important factor responsible for non-stationarity in all 12 sub-periods models. At the Gorgan station, however, the time covariate played a key role in the models of 3 out of 12 sub-periods. In the other two stations (Bandar Anzali and Yazd), the time covariate was not included in the optimal In-depth analysis of the results of other stations showed (figure not presented) that wherever the covariates had a clear impact on P and AED (the two components needed for the calculation of the indices), there was a substantial difference between SSPEI and NSPEI. As a result, the difference between the time series of SSPEI and NSPEI had no association with the stations' climate in the study area. Figure 4 shows the temporal distribution of drought/wet events severity and duration for SSPEI and NSPEI for all selected stations from 1964 to 2014. On the basis of SSPEI, the figure clearly shows the occurrence of a wet period during the first half of the 1990s and a drought period from 1998 to 2014 (except 2004) at about all stations of interest. However, according to NSPEI, the long-term drought period of 1998-2014 has been broken into several minor droughts and wet periods, especially after 2004. For the decades before the 1990s, although a dominant pattern of drought/wet period is not seen in the figure, the drought events in the first decade of record period (about 1964-1974) and the wet events in the second and third decades (about 1975-1990) are more common. Therefore, two key points from the figure may be inferred: 1) NSPEI modifies the severity of drought events, especially the long-term drought period of 1998-2014; and 2) NSPEI breaks the long-term drought/wet periods into several short-term periods. For example, at the Abadan station (the station shown with the number 1), while SSPEI identifies a 17-year (1964-1980) wet period, NSPEI breaks this continuous wet period into several minor droughts and wet periods. Figure 5 displays the spatiotemporal distribution of SSPEI and NSPEI over Iran during the water years (October-September) of 1964-2014. As shown in the figure, the difference between SSPEI and NSPEI in some years is relatively high. Focusing on drought events, NSPEI monitored more intensive and extensive droughts than SSPEI, especially in the water years 1965years , 1969years , and 1972years . However, in 1998years , 1999years , 2007years , and 2010, the intensity and extent of droughts quantified by SSPEI is higher than NSPEI. Among them, the two water years, 1998 and 2010, are significantly different in terms of the zonation maps of SSPEI and NSPEI. In 1998, while the SSPEI indicates that almost the entire country is under threat of drought, the NSPEI indicates that the eastern half of the country is experiencing above-normal conditions. Considering the 2010, the difference between the maps of SSPEI and NSPEI is higher than 1998. On the basis of SSPEI, almost all regions of the country were affected by drought in 2010, while NSPEI calculated above-normal conditions in the majority of the country in the same year.

The Vegetation Response to the Drought Indicators
The applicability of remote sensing data (such as NDVI and VCI) for drought monitoring in Iran has been considered by several researchers (Rahimzadeh Bajgiran et al. 2008;Shamsipour et al. 2008). The results of the studies conducted in Iran showed that the temporal and spatial characteristics of drought could be detected and mapped using the indices derived from satellite images.
In this study, the response of vegetation cover to the severest drought (the 2008 drought) monitored by SSPEI (as the base/traditional drought index) was investigated and compared to NSPEI. The monthly averages of VCI, SSPEI, and NSPEI throughout the studied stations were calculated for the period of September 2007 to November 2009, as displayed in Fig. 6. As shown in the figure, VCI identified a dry period in vegetation from May 2008 to March 2009, with the lowest value of 16% in November 2008. The ground-based drought indices (SSPEI and NSPEI) simultaneously discerned a dry period beginning in February 2008, about three months prior to when the satellite-based drought index (VCI) reported it. Although the two indicators were almost identical in terms of identifying the onset of the 2008 drought, they differed with regard to characteristics such as the length of the drought period, the intensity and time of the drought peak, and the termination of the drought event. Even several months after the drought termination identified by NSPEI and VCI, SSPEI does not appear to reach the above-normal condition, particularly in the latter case.After Fig. 6 The SSPEI, NSPEI, and VCI averaged over all stations from the September 2007 to November 2009. The reference line is at zero for the left y-axis and at 40 for the right y-axis the drought termination, both VCI and NSPEI remain above their normal condition, but SSPEI is still below its normal condition. On the other hand, in November 2008, when VCI is in its lowest value, NSPEI detects a moisture signal which is much stronger than that of SSPEI. In response to the signal four months later, VCI returns to normal.. Thus, one can conclude the superiority of NSPEI, as a better indicator of vegetation response to meteorological drought than SSPEI and the necessity for the development of non-stationary indicators for reliable and robust estimation of drought events.

Conclusion
In this study, a Non-stationary Standardized Precipitation Evapotranspiration Index (NSPEI) was developed for reliable and robust monitoring of meteorological drought characteristics under a changing climate in Iran. The NSPEI introduces a non-stationary, instead of stationary, log-logistic probability distribution in the mathematical structure of traditional SPEI SSPEI). The non-stationary distribution assumes that its location parameter is a multivariate function of external covariates, including time and several climate indices. The optimal non-stationary function was determined using a forward selection method in the GAMLSS algorithm. The NSPEI was constructed at 32 weather stations across Iran for the period of 1964-2014. The proposed NSPEI was compared to the SSPEI in different aspects. The key findings of this study are: -The non-stationary log-logistic distribution fitted the water deficits and surpluses series (i.e., the P minus AED series) better than the stationary log-logistic distribution at almost all chosen stations. -The location parameter of non-stationary log-logistic distribution is mostly influenced by the two teleconnections, ENSO and NCP. The time covariate had more relative importance than the other two teleconnections (i.e., NAO and AO) in simulating the temporal variations of location parameter. -The SSPEI monitored the drought events with durations and severities larger than those identified by NSPEI. Also, the long-term drought periods (e.g., the 1998-2014 droughts except 2004) quantified by SSPEI, were broken into several discrete drought periods when using NSPEI. -There were clear differences between SSPEI and NSPEI in terms of the spatial maps of drought and wet events during the record period of 1964-2014. -During the severest 2008 drought, the temporal variations of the average vegetation cover (quantified by average VCI) over the studied stations were more closely related to those of the average NSPEI than the average SSPEI. NSPEI successfully identified the 2008 drought characteristics in vegetation cover in terms of both the drought duration and termination. Although the drought index developed in this work was successful in modeling nonstationarity in drought behavior, it still has several limitations that may be addressed in future research: -In this study, only non-stationarity in the log-logistic distribution's location parameter was modeled. It is suggested that non-stationarity in the other parameters (scale and shape) of the aforesaid distribution be investigated in order to improve the results of quantifying drought and wet events and their consequences on vegetation.
-The non-stationary model used in this study was linear. It is suggested that nonlinear models be examined within the context of the non-stationary algorithm and their capabilities be compared to linear models.