Identification of Hydrological Drought in Eastern China Using a Time-Dependent Drought Index

Long records (1960–2013) of monthly streamflow observations from 8 hydrological stations in the East Asian monsoon region are modeled using a nonstationarity framework by means of the Generalized Additive Models in Location, Scale and Shape (GAMLSS). Modeling analyses are used to characterize nonstationarity of monthly streamflow series in different geographic regions and to select optimal distribution among five two-parameter distributions (Gamma, Lognormal, Gumbel, Weibull and Logistic). Based on the optimal nonstationarity distribution, a time-dependent Standardized Streamflow Index (denoted SSIvar) that takes account of the possible nonstationarity in streamflow series is constructed and then employed to identify drought characteristics at different time scales (at a 3-month scale and a 12-month scale) in the eight selected catchments during 1960–2013 for comparison. Results of GAMLSS models indicate that they are able to represent the magnitude and spread in the monthly streamflow series with distribution parameters that are a linear function of time. For 8 hydrological stations in different geographic regions, a noticeable difference is observed between the historical drought assessment of Standardized Streamflow Index (SSI) and SSIvar, indicating that the nonstationarity could not be ignored in the hydrological drought analyses, especially for stations with change point and significant change trends. The constructed SSIvar is, to some extent, found to be more reliable and suitable for regional drought monitoring than traditional SSI in a changing environment, thereby providing a feasible alternative for drought forecasting and water resource management at different time scales.


Introduction
Drought is an insidious natural hazard that has damaging and costly impacts on agriculture, water resources, ecology, and society [1,2].Due to global climate change and rapid socio-economic development, drought disasters have occurred more frequently and extensively in recent decades [3].Drought is commonly defined as below-normal water availability [4] and can be subdivided into different types of drought related to the variables of the hydrological cycle, precipitation (meteorological drought), soil moisture (agricultural drought), and streamflow (hydrological drought) [5].
Hydrological drought is defined as a significant shortage of availability of water in all its forms appearing in the land phase of the hydrological cycle (e.g., streamflow, groundwater level and lake level) [6].It is traditionally detected using field observations of streamflow, surface water and groundwater levels, thereby providing direct evidence of any below-normal water availability [6,7].Below-normal water availability in rivers, lakes and reservoirs can cause water scarcity in combination with water demand, threating water supply and associated food production [8].Thus, it is important to identify hydrological drought characteristics and assess the effects of hydrological drought quantitatively.The Standardized Streamflow Index (SSI), proposed by Shukla and Wood [9], has been used as a useful index for characterizing hydrological drought.It is constructed using streamflow data based on the concept of the standardized precipitation index (SPI).Similar to SPI, with the advantages of computational simplicity, the SSI is capable of characterizing hydrological drought condition at different time scales [3,10], and it enables the severity of hydrological drought in different locations to be compared independently of the local characteristics [11].In addition, the SSI is sensitive to the factors and assumptions that govern probabilistic hydrology, since it is a probability-based drought index [6].
Generally, statistical inferences and statistical analyses for hydrologic time series have relied heavily on the assumption of stationarity in hydrology [12,13].Under the stationarity assumption, hydrological series keep their distributional properties invariant with time, implying lack of trends and shift [14].However, the stationarity assumption has been widely questioned and should no longer serve as a central, default assumption, as a result of global climatic change and man-induced disturbance [15][16][17].In recent years, hydrological nonstationarity has drawn considerable attention [14,18].Coulibaly and Baldwin [18] developed an optimal dynamic recurrent neural networks method to directly forecast hydrological series under nonstationarity conditions, and they found that neural networks are good alternatives for modeling the complex dynamics of the hydrological system.Villarini et al. [14] modeled a long record of seasonal rainfall and temperature in a nonstationarity framework to characterize non-stationarities in hydro-climatic variables.
In the traditional calculation of SSI values, the streamflow sample series are firstly fitted to a suitable stationary probability distribution on the basis of stationarity assumption [11].This means that historical features of streamflow series can be used to derive SSI values in the future.However, changing environments (e.g., climate change and human activities) might alter the statistical characteristics of hydroclimate time series [19], resulting in so-called nonstationarity.With respect to streamflow, there have been numerous studies on individual impacts of climate change (mainly changes in precipitation and temperature) and human activity (mainly water construction and building of dams) on it [20], implying clear violations of the stationarity assumption.Ignoring the nonstationarity would therefore most likely diminish the availability and validity of traditional SSI in hydrological drought analysis and could lead to the underestimation or overestimation of the drought severity [21].Thus, it is essential to incorporate the nonstationarity of streamflow sample series in constructing an appropriate variant SSI (SSI var ) under nonstationarity conditions, thereby providing significant information for evaluation and mitigation of risk of hydrological drought hazards and management of water resources.From this point of view, although some researchers have attempted to consider nonstationarity in developing drought index for drought monitoring, most research has only focused on mean time variance of precipitation time series for meteorological drought (e.g., [22,23]).However, to date, relatively few studies have addressed the stationarity or nonstationarity of streamflow series for hydrological drought using a nonstationarity framework in different geographic regions considering both trends and change points in the parameters.This constituted the major motivation of this study.
Several methods have been proposed to model nonstationarity time series in the previous literature (e.g., [24,25]), each having their own strengths and weakness.The Generalized Additive Models for Location, Scale and Shape (GAMLSS), proposed by Rigby and Stasinopoulos [26], has recently gained popularity in modeling nonstationarity time series in hydrology [27].This model provides a high degree of flexibility in addressing nonstationarity probabilistic modeling.In GAMLSS, the assumption that the variable of interest follows a distribution from the exponential family is relaxed, allowing the use of more general distributions, such as highly skewed or kurtotic distributions, which may be more appropriate for modeling the record of interest.This makes it an appealing framework for nonstationarity modeling of hydrometeorological variables to improve the existing drought index.
In this study, the GAMLSS was used to model streamflow with nonstationarity distribution to construct a variant SSI index (SSI var ) in eight catchments in the eastern region of China.The main objectives of this study are: (1) to analyze whether streamflow series are stationary in the eastern region of China, and (2) to construct an appropriate variant SSI that accounts for the changes in the parameters of the selected distribution under nonstationarity conditions.The results of this study could provide important information for the management of water resources and evaluation of hydrological drought hazards under a changing environment.

Traditional Standardized Streamflow Index (SSI)
The Standardized Streamflow Index (SSI) was used as the useful counterpart for depicting hydrological aspects of drought.This index has the same theoretical background, transforming monthly streamflow into z-scores.The procedure of traditional SSI calculation is statistically similar to the well-known Standardized Precipitation Index (SPI).
The long-term streamflow series were first fitted to the lognormal distribution (Table 1), as suggested by Nalbantis and Tsakiris [6].Once the distribution is determined, the SSI (in z-scores) can easily be calculated by following the classical approximation of Abramowitz and Stegun [28].For example where F(x) is the cumulative distribution function.The constants are C 0 = 2.515517, C 1 = 0.802853, C 2 = 0.010328, d 1 = 1.432788, d 2 = 0.189269, d 3 = 0.001308.Table 1 shows the range of SSI values along with their classifications [6].As a nonparametric statistical test that allows detection of changes in the mean when the change point time is unknown, the Pettitt test [29] has been suggested by Villarini et al. [12] for analyzing the change point in this study.Mathematically, when a sequence of random variables is divided into two segments represented by X 1 , . . ., X m and X m+1 , . . ., X n , if each segment has a common distribution function, then the change point is identified at m.The Pettitt test uses a version of the Mann-Whitney statistic (U t,n ).The breakpoint is defined to be where |U t,n | reaches its maximum value: The significance level associated with K n is calculated as This test is valid for continuous variables, and its null hypothesis is the absence of a change point.If p < 0.05, a significant change point exists.

Temporal Trend Analysis
The Mann-Kendall (M-K) statistical test [30,31] is commonly used for general use in temporal trend analysis due to its robustness for non-normally distributed and censored data, which are frequently encountered in hydro-climatic time series [32].However, the presence of serial correlation can complicate the identification of trends, in that a positive serial correlation can increase the expected number of false positive outcomes for the Mann-Kendall test.Any serial correlation should be removed before conducting the M-K trend test.Following Zhang et al. [33], the significant trend for the streamflow records was identified using the following steps: (1) compute the lag1 serial correlation coefficient (ρ 1 ); (2) if ρ 1 < 0.1, the M-K test was applied to the streamflow series directly; otherwise (3) the M-K test was applied to the preprocessed time series ( . The 95% confidence level was used to evaluate the significance of trends.

GAMLSS Model
The Generalized Additive Models for Location, Scale and Shape (GAMLSS), proposed by Rigby and Stasinopoulos [26], have been used in our study to model the streamflow series with nonstationarity probability distribution.The GAMLSS assume a parametric distribution for the response variable y (streamflow in our study) and model the parameters of the distribution as linear and/or nonlinear, parametric and/or additive nonparametric functions of explanatory variables [12].
A brief description of the GAMLSS model is provided below, with more detailed information available from Rigby and Stasinopoulos [26] and Stasinopoulos and Rigby [34].In a GAMLSS model, it is assumed that there are independent random variables y i , for i = 1, . . ., n, which are fitted to a distribution function of f y i θ i conditional on θ i = (θ i 1 , . . ., θ i p ), a vector of p distribution parameters accounting for location, scale, and shape.Generally, p is less than or equal to four, since one, two, three and four parameter distribution families could provide enough flexibility for most applications in hydrology [35].The GAMLSS allow for a general distribution function, including highly skewed and/or kurtotic continuous or discrete distributions.The distribution parameters are related to the design matrix of explanatory variables by monotonic link function g k (•), for k = 1, . . ., p. Similar to Villarini et al. [12], five commonly used two-parameter distributions, i.e., Gamma distribution (GA), Lognormal distribution (LOGNO), Gumbel distribution (GU), Weibull distribution (WEI) and Logistic distribution (LO), were used in this study to model the streamflow series (Table 2).In this study, taking the time as the explanatory variables, the linear functions relating t and parameters θ 1 (for mean µ) and θ 2 (for variance σ) were constructed in the form of: where β 1 and β 2 denote the vectors of coefficients of the linear function respectively, and the a 1 and a 2 is the corresponding constant term.The Akaike Information Criterion (AIC) [36] was used to select the optimal nonstationarity model with the highest goodness-of-fit, and the model with the minimum AIC value was selected.Additionally, to further assess the performance of the selected optimal model, the residuals of each model was checked by the statistics of the Filliben correlation coefficient (denoted by F r ) [37], together with the worm plot [38] as a visual inspection of diagnostic plots of the residuals.AIC and F r are calculated as AIC = −2 ln(ML) + 2k ( 6) where ML is the maximum likelihood function of models and k is the number of independently adjusted parameters within the model; S i are the ordered residuals; n is the length of the observation period and B i are the standard normal order statistic medians.
Then, based on the GAMLSS, nonstationarity models with different probability distributions, trends in the parameters and the changes point in the time series of streamflow were compared to select an optimal model to construct a newly Standardized Streamflow Index (SSI) for drought monitoring and quantification.Analysis and calculations related to the GAMLSS model in this study were performed with R-based GAMLSS package [34].

Construction of the Time-Dependent Standardized Streamflow Index (SSI var )
In our study, to construct a time-dependent Standardized Streamflow Index (SSI var ), the different nonstationarity probability distributions with their parameters changing over time were compared by using GAMLSS, as described in Section 2.3.Considering the long-term trends and abrupt changes can reveal the presence of the nonstationarity in streamflow series.The change points and trends were also included in the analysis by GAMLSS and AIC scores were used to select the optimal models for the Standardized Streamflow Index.Here, synthetic experiments were designed to construct the time-dependent SSI (SSI var ) drought index.Firstly, to select the optimal model for the hydrological stations with no change point, four different models have been analyzed that: The constructed SSI var is calculated in the following steps: (1) calculate the monthly average streamflow (Y t ) for a given k-month scale (3-month scale and 12-month scale in this study) with respect to time t; (2) within the GAMLSS framework, the nonstationarity model is developed by fitting the streamflow series (Y t ) to a nonstationarity optimal distribution which is selected by minimizing AIC, then obtaining the corresponding cumulative probability F(x) k .The parameters of the distribution is described as a linear function of time; (4) the cumulative probability F(x) k is converted to a standard normal deviate (with zero mean and unit variance), which also using the approximate conversion provided by Abramowitz and Stegun [28].This standard normal value is the SSI var index.
Positive values of the SSI var imply wet conditions, while negative values indicate dry conditions.Since the SSI var are developed based on the concept of SSI, the classifications listed in Table 2 can be also used for the SSI var as well.Compared with the traditional SSI, the proposed SSI var recognize the nonstationarity in the streamflow series and thus it has the ability to capture and model the nonstationarity to provide more reasonable and satisfactory for drought monitoring and drought analysis.

Study Area and Dataset
The eastern region of China (the area to the east of 100 • E), which is affected by the East Asian monsoon, is usually referred to as the East Asian monsoon (EAM) region.Rainfall mainly occurs in the summer season (June, July and August) in the EAM region.The EAM region is a critically important zone in China, having the densest population and exhibiting the fastest economic development in recent decades.In our study, we selected eight different catchments in the EAM region of China as the case study areas (Figure 1), including Songhua River Basin (SHRB), Hai River Basin (HRB), Yellow River Basin (YRB), Huai River Basin (HURB), Yangtze River Basin (YARB), Pearl River Basin (PRB).
Monthly streamflow records from 1960 to 2013, from eight hydrological stations in different catchments, were analyzed.The data were provided by the Hydrology Bureau of the corresponding River Conservancy Commission.Locations of these stations are shown in Figure 1.Information on the data, such as the length of streamflow series and the drainage areas of hydrological stations, is given in Table 3.

Change Point Analysis
Pettitt's test was used to detect change points in annual streamflow for the eight selected hydrological stations.The results showed that the annual streamflow in three stations (Luanxian, Huaxian and Tiane) have change points with a significance level of 5% (Figure 2), while the remaining five stations showed no significant change points.The change points for Luanxian, Huaxian and Tiane occurred in 1979, 1990 and 1986, respectively (Figure 2).

Change Point Analysis
Pettitt's test was used to detect change points in annual streamflow for the eight selected hydrological stations.The results showed that the annual streamflow in three stations (Luanxian, Huaxian and Tiane) have change points with a significance level of 5% (Figure 2), while the remaining five stations showed no significant change points.The change points for Luanxian, Huaxian and Tiane occurred in 1979, 1990 and 1986, respectively (Figure 2).

Trend Analysis
Autocorrelation analysis was done first to detect the significant serial effects before trend detection.Apart from Waizhou (lag1 = −0.073)and Boluo (lag1 = −0.066), the lag1 serial correlation coefficient 1 r for the other six stations were all beyond the threshold value (0.1) (Table 4), indicating significant effects on trend analysis.The significant autocorrelations in the six stations were removed before the M-K trend test.The results estimated by M-K test are exhibited in Table 4.Among the eight stations, two stations showed an increasing trend and the other six stations demonstrated a decreasing trend.The annual streamflow series at three stations (Luanxian, Huaxian and Tiane) showed significantly decreasing trends.As for the other five stations, though the M-K trends were insignificant, they also showed slight decreasing or increasing trends.

Trend Analysis
Autocorrelation analysis was done first to detect the significant serial effects before trend detection.Apart from Waizhou (lag1 = −0.073)and Boluo (lag1 = −0.066), the lag1 serial correlation coefficient ρ 1 for the other six stations were all beyond the threshold value (0.1) (Table 4), indicating significant effects on trend analysis.The significant autocorrelations in the six stations were removed before the M-K trend test.The results estimated by M-K test are exhibited in Table 4.Among the eight stations, two stations showed an increasing trend and the other six stations demonstrated a decreasing trend.The annual streamflow series at three stations (Luanxian, Huaxian and Tiane) showed significantly decreasing trends.As for the other five stations, though the M-K trends were insignificant, they also showed slight decreasing or increasing trends.

Modeling with GAMLSS
It can be seen from Table 5 that the nonstationarity models performed best for fitting the streamflow series at different time scales for the eight hydrological stations.At both 3-month and 12-month time scales, the gamma distribution was the optimal distribution for half of the selected stations (4 stations), followed by the lognormal and Weibull distribution.With respect to testing the stationarity, no models were found to be stationary with time-independent parameters (Table 5).Among the 16 selected models at different time scales, 12 models exhibited nonstationarity in θ 1 (for mean), and 4 models exhibited nonstationarity in θ 2 (for variance).No models of nonstationarity in mean and variance were selected.The statistic of F r (Table 5) also indicated that the selected nonstationarity model was an adequate fit for the streamflow series at different time scales at the respective stations.The averaged streamflow series in August was taken as an example to show the fitting of the GAMLSS model.Fitting of averaged streamflow series in August at different time scales for the eight stations using GAMLSS is shown in Figure 3 (3-month scale) and Figure 4 (12-month scale).The vast majority of the points were within the 0.05 and 0.95 quantiles (Figures 3 and 4), indicating that the selected models were able to capture the variability of the data at both 3-month scale and 12-month scale.The centile curves of the GAMLSS modeling in eight stations exhibited temporal trend, showing that the fitted models captured non-linear behaviors associated with the streamflow series.For the three stations (Luanxian, Huaxian and Tiane) with change points, the 5% and 95% percentile curves were significantly impacted by the abrupt changes.Visual inspections of worm plots were also conducted to check the residuals.All the worm points were within the 95% confidence intervals (Figures 5 and 6), indicating consistency between the selected nonstationarity model and the observed data at different time scales.Thus, the results indicated the GAMLSS-based modeling was able to capture the temporal variability of the streamflow series.
. Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 3-month scale.The red points are the observed streamflow at 3-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.

Figure 3.
Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 3-month scale.The red points are the observed streamflow at 3-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.

Figure 3.
Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 3-month scale.The red points are the observed streamflow at 3-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.
Figure 5. Worm plots for the residuals from the GAMLSS models for streamflow series as illustrated in Figure 3.For a good fit, the data points should be aligned, preferably along the red solid line, but within the 95% confidence intervals indicated by the two grey dashed line.Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.
Figure 5. Worm plots for the residuals from the GAMLSS models for streamflow series as illustrated in Figure 3.For a good fit, the data points should be aligned, preferably along the red solid line, but within the 95% confidence intervals indicated by the two grey dashed line.

Construction of the SSIvar
Based on the GAMLSS model, both the time-dependent SSI (SSIvar) and traditional SSI were employed to identify drought characteristics in the selected eight basins during 1961-2013.It can be seen from Figure 7 that noticeable differences existed in the values of cumulative probability F(x) for calculating SSI and SSIvar.The deviations of F(x) for SSI and SSIvar for different hydrological stations were varying in degree (Figure 7).Taking Luanxian and Wangjiaba stations as examples, the matching points were distributed around 1:1 line in Wangjiaba while the points were deviated from 1:1 line in Luanxian.The root mean square error (RMSE) of Luanxian (0.189) was larger than that in Wangjiaba (0.052).This may be because the streamflow series in Luanxian has change point and significant variation trend.During the procedure of SSI calculation, the cumulative probability F(x) was transformed to a standard normal deviation with a zero mean and unit variance, which is the value of SSI.Thus, the deviations of F(x) for calculating SSI and SSIvar could result in the differences between SSI and SSIvar.

Construction of the SSI var
Based on the GAMLSS model, both the time-dependent SSI (SSI var ) and traditional SSI were employed to identify drought characteristics in the selected eight basins during 1961-2013.It can be seen from Figure 7 that noticeable differences existed in the values of cumulative probability F(x) for calculating SSI and SSI var .The deviations of F(x) for SSI and SSI var for different hydrological stations were varying in degree (Figure 7).Taking Luanxian and Wangjiaba stations as examples, the matching points were distributed around 1:1 line in Wangjiaba while the points were deviated from 1:1 line in Luanxian.The root mean square error (RMSE) of Luanxian (0.189) was larger than that in Wangjiaba (0.052).This may be because the streamflow series in Luanxian has change point and significant variation trend.During the procedure of SSI calculation, the cumulative probability F(x) was transformed to a standard normal deviation with a zero mean and unit variance, which is the value of SSI.Thus, the deviations of F(x) for calculating SSI and SSI var could result in the differences between SSI and SSI var .8 and 9, it can be seen that similar differences between SSI and SSIvar were observed in four basins at different time scales.For the selected basins, the drought severity shown by SSIvar was first more severe than traditional SSI, and then it was subsequently milder than traditional SSI.The stations with change points and significant trend changes (i.e., Luanxian, Huanxian and Tiane) have larger deviations between SSI and SSIvar than the other stations.For instance, for the basin above Luanxian station at 3-month scale (Figure 8a), obviously lower values were detected by SSIvar during 1961-1968 compared to traditional SSI, while notably higher values were found for SSIvar during 1992-2013.For the basin above Huaxian station at 12-month scale (Figure 9b), we can also observe that the values of SSI were higher than SSIvar during 1961-1990 (e.g., moderate drought was detected by SSIvar-12 in May 1972 (−1.43), while slight drought was characterized by SSI-12 (−0.76)), while being lower than SSIvar during 1991-2013.As the SSI is a standardized variate, the average values of the SSI and the standardized deviation must equal 0 and 1, respectively.Meanwhile, the two-sample Kolmogorov-Smirnov (KS) test [39] was used to test whether two indexes cons from same distribution.The KS test statistically showed that SSI and SSIvar are statistically different in Luanxian, Huaxian and Tiane stations at 95% confidence level (p < 0.05) (Figures 8 and 9).It is worth noting that the boxplots of two SSI series show that the symmetry of the SSIvar is superior to that of the SSI at different time scales (Figures 8 and 9), and it is closer to the standard normal distribution (zero mean and unit variance).Thus, compared to the traditional SSI, the SSIvar is capable of modeling the nonstationarity of streamflow series and is more appropriate for hydrological drought assessment within a changing environment.8 (3-month scale) and Figure 9 (12-month scale).From Figures 8 and 9, it can be seen that similar differences between SSI and SSI var were observed in four basins at different time scales.For the selected basins, the drought severity shown by SSI var was first more severe than traditional SSI, and then it was subsequently milder than traditional SSI.The stations with change points and significant trend changes (i.e., Luanxian, Huanxian and Tiane) have larger deviations between SSI and SSI var than the other stations.For instance, for the basin above Luanxian station at 3-month scale (Figure 8a), obviously lower values were detected by SSI var during 1961-1968 compared to traditional SSI, while notably higher values were found for SSI var during 1992-2013.For the basin above Huaxian station at 12-month scale (Figure 9b), we can also observe that the values of SSI were higher than SSI var during 1961-1990 (e.g., moderate drought was detected by SSIvar-12 in May 1972 (−1.43), while slight drought was characterized by SSI-12 (−0.76)), while being lower than SSI var during 1991-2013.As the SSI is a standardized variate, the average values of the SSI and the standardized deviation must equal 0 and 1, respectively.Meanwhile, the two-sample Kolmogorov-Smirnov (KS) test [39] was used to test whether two indexes cons from same distribution.The KS test statistically showed that SSI and SSIvar are statistically different in Luanxian, Huaxian and Tiane stations at 95% confidence level (p < 0.05) (Figures 8 and 9).It is worth noting that the boxplots of two SSI series show that the symmetry of the SSI var is superior to that of the SSI at different time scales (Figures 8 and 9), and it is closer to the standard normal distribution (zero mean and unit variance).Thus, compared to the traditional SSI, the SSI var is capable of modeling the nonstationarity of streamflow series and is more appropriate for hydrological drought assessment within a changing environment.

Conclusions
In this study, streamflow series at eight different stations in the EAM region of China were modeled with nonstationarity distributions using GAMLSS.Based on the constructed nonstationarity distributions, the SSI var was developed and used to estimate regional drought characteristics in the selected basin.The main conclusions can be summarized as follows: (1) Generalized Additive Models in Location, Scale and Shape (GAMLSS) provide a flexible and useful framework for modeling distributions of streamflow series considering both trend and change point.In particular, they provide the capability for modeling the non-stationarities in streamflow records.(2) Based on the selected optimal distribution, the developed SSI var is capable of taking the nonstationarity of streamflow series into account; thus, it is likely to be more reliable and suitable than the traditional SSI for drought assessment in a changing environment.The differences between the SSI var and SSI indicate that the presence of nonstationarity should be considered in regional drought assessment.The SSI var is proven to be a feasible alternative for drought forecast and water resource management under changing environment.
(1) stationary model; (2) nonstationarity model in parameter for mean; (3) nonstationarity model in parameter for variance; (4) nonstationarity model in mean and variance.Likewise, for the streamflow series with change points, the study period could be divided into two periods by the change point and four different models are discussed: (1) stationary model; (2) change point only exists in mean; (3) change point only exists in variance; (4) change point exists in both mean and variance.

Figure 1 .
Figure 1.The study area in EAM region.

Figure 1 .
Figure 1.The study area in EAM region.

Figure 2 .
Figure 2. Pettitt's test for detecting a change in the mean of annual streamflow.Horizontal lines represent the significance (solid line represent 1% and dotted line represent 5%).

Figure 2 .
Figure 2. Pettitt's test for detecting a change in the mean of annual streamflow.Horizontal lines represent the significance (solid line represent 1% and dotted line represent 5%).

Figure 4 .
Figure 4. Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.

Figure 4 .
Figure 4. Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.

Figure 4 .
Figure 4.Centile curve plots for assessing the performance of the optimal nonstationarity model using time as covariate at 12-month scale.The red points are the observed streamflow at a 12-month scale in August.The black line is the 50% centile curve, the dark grey region is the area between the 25% and 75% centile curves, and the light grey region is the area between the 5% and 95% centile curves.

Figure 5 .
Figure 5. Worm plots for the residuals from the GAMLSS models for streamflow series as illustrated in Figure3.For a good fit, the data points should be aligned, preferably along the red solid line, but within the 95% confidence intervals indicated by the two grey dashed line.

Figure 6 .
Figure 6.Worm plots for the residuals from the GAMLSS models for streamflow series as illustrated in Figure4.For a good fit, the data points should be aligned, preferably along the red solid line, but within the 95% confidence intervals indicated by the two grey dashed line.

Figure 6 .
Figure 6.Worm plots for the residuals from the GAMLSS models for streamflow series as illustrated in Figure4.For a good fit, the data points should be aligned, preferably along the red solid line, but within the 95% confidence intervals indicated by the two grey dashed line.

Figure 7 .
Figure 7.Comparison of the cumulative probability F(x) between SSI and time-dependent SSI (SSIvar) for eight hydrological stations during 1961-2013.Taking Luanxian, Huaxian, Danjiangkou and Tiane as examples, to show the influence of change of parameter clearly, the SSIvar and traditional SSI during 1960-2013 in four basins are shown in Figure 8 (3-month scale) and Figure 9 (12-month scale).From Figures8 and 9, it can be seen that similar differences between SSI and SSIvar were observed in four basins at different time scales.For the selected basins, the drought severity shown by SSIvar was first more severe than traditional SSI, and then it was subsequently milder than traditional SSI.The stations with change points and significant trend changes (i.e., Luanxian, Huanxian and Tiane) have larger deviations between SSI and SSIvar than the other stations.For instance, for the basin above Luanxian station at 3-month scale (Figure8a), obviously lower values were detected by SSIvar during 1961-1968 compared to traditional SSI, while notably higher values were found for SSIvar during 1992-2013.For the basin above Huaxian station at 12-month scale (Figure9b), we can also observe that the values of SSI were higher than SSIvar during 1961-1990 (e.g., moderate drought was detected by SSIvar-12 in May 1972 (−1.43), while slight drought was characterized by SSI-12 (−0.76)), while being lower than SSIvar during 1991-2013.As the SSI is a standardized variate, the average values of the SSI and the standardized deviation must equal 0 and 1, respectively.Meanwhile, the two-sample Kolmogorov-Smirnov (KS) test[39] was used to test whether two indexes cons from same distribution.The KS test statistically showed that SSI and SSIvar are statistically different in Luanxian, Huaxian and Tiane stations at 95% confidence level (p < 0.05) (Figures8 and 9).It is worth noting that the boxplots of two SSI series show that the symmetry of the SSIvar is superior to that of the SSI at different time scales (Figures8 and 9), and it is closer to the standard normal distribution (zero mean and unit variance).Thus, compared to the traditional SSI, the SSIvar is capable of modeling the nonstationarity of streamflow series and is more appropriate for hydrological drought assessment within a changing environment.

Figure 7 .
Figure 7.Comparison of the cumulative probability F(x) between SSI and time-dependent SSI (SSI var ) for eight hydrological stations during 1961-2013.

Table 1 .
Drought classifications based on the traditional SSI.

Table 2 .
Summary of the distributions used to model the streamflow series at different time scales in our study.

Table 3 .
Information on the hydrological stations considered in this study.

Table 3 .
Information on the hydrological stations considered in this study.

Table 4 .
Results of analysis of trends in streamflow series for the eight hydrological stations at an annual scale.

Table 5 .
Summary of results for the GAMLSS models of eight hydrological stations in EAM region (the critical values of the Filliben correlation coefficient is F α = 0.978 and the F r bigger than F α , indicating that the nonstationarity model passes the goodness-of-fit test).