Assessment of hydrological drought based on nonstationary runoff data

Drought is a kind of universal natural disaster which can bring serious ecological, environmental and social consequences. It may inhibit the growth of vegetation, accelerate the degradation of grassland and the disappearance of oasis, aggravate the reduction of rivers and lakes, and lead to the ecosystems more vulnerable (Zhang et al. 2010; Cheng 2012; Fang et al. 2018a; Han et al. 2018a). Thus, drought has attracted wide attention from the researchers. The accurate and timely drought assessment is crucial for devising plans and carrying out the necessary mitigation measures especially under the background of climate warming.

Being a robust tool for drought assessment, many drought indices have been developed during the last decades (e.g., Nalbantis & Tsakiris 2009; Vicente-Serrano et al. 2012; Bloomfield & Marchant 2013; Tabari et al. 2013). One of the most widely used indices worldwide in the hydrological drought assessment is the standardized runoff index (SRI; Luciano et al. 2012; Jiang et al. 2019; Zhou et al. 2019), which is proposed by Shukla & Wood (2008) based on the theory of the standard precipitation index (SPI; McKee et al. 1993). It remains the most broadly accepted index due to being simple to calculate, the limited data required, and the different time scales of hydrological drought assessed. However, the stationarity assumption of the SRI is often violated in the context of global warming and strong human disturbance (Xiong & Guo 2004; Villarini et al. 2009; Liu et al. 2019). For example, Xiong & Guo (2004) found that the annual mean runoff series in the Yangtze River basin in China exhibited a clear decreasing trend from 1882 to 2001. Villarini et al. (2009) found that the runoff series in the continental United States had significant change points during the 20th century. These mean that many runoff series are not stationary any more, and the SRI with the assumption of the runoff data to be analyzed keeping stationary is no longer suitable for fitting such nonstationarity data (Wang et al. 2019b). Therefore, a modified SRI, namely, a nonstationary standardized runoff index (NSRI), will be proposed in this study with incorporating the nonstationary characteristic of runoff series.

NSRI is defined similarly with the SRI, but using a nonstationary Gamma distribution by incorporating covariates. Seeking appropriate covariates is the first, also one of the vital steps in constructing the nonstationary drought index (Vu & Mishra 2019). A single covariate is usually considered in previous literatures. For example, time or climatic variables are usually used to describe the precipitation changes in developing the nonstationary meteorological drought index (Russo et al. 2013; Li et al. 2015; Wang et al. 2015a; Javad & Somayeh 2018). However, as known, the runoff changes are always affected by the recharge sources. The recharge sources for the upper reach of the Heihe River basin include precipitation, glacial and snow melting, frozen soil thawing and groundwater (Ding et al. 1999; Wang et al. 2011). These sources are affected by meteorological elements like temperature, relative humidity and wind speed. Therefore, these meteorological elements are considered as covariates herein. In addition, evapotranspiration is an important process of water transfer in the hydrosphere and atmosphere, playing a vital role in the hydrological cycle and the runoff processes (Han et al. 2018b; Aa et al. 2019). The antecedent runoff (runoff in the previous month) contributes to the runoff changes at that month especially in non-flood seasons (Ding et al. 1999). Consequently, the evapotranspiration and the antecedent runoff are also taken into account as alternative covariates. As the actual evapotranspiration data are not available and not easy to calculate, potential evapotranspiration is used here. In summary, the six variables, precipitation, temperature, relative humidity, wind speed, potential evapotranspiration, and antecedent runoff, are incorporated in developing the NSRI.

This study, taking the upper reach of the Heihe River basin (the second largest inland river basin in northwest of China) as the study area, first tests the stationarity of runoff series over the study area and then develops an NSRI to assess hydrological droughts under nonstationary conditions by using GAMLSS. GAMLSS is one of the techniques for simulating nonstationary time series and capable of providing a high degree of flexibility for solving nonstationary modelling (e.g., Villarini et al. 2009; Zheng et al. 2018). It can describe the linear or nonlinear relationship between any statistical parameter of the random variable sequence and the explanatory variable, and the description of the random variable distribution function has a wide range of types (Stasinopoulos & Rigby 2007). Many hydrologists have applied this technique to analyze the nonstationary hydrological sequences (Villarini et al. 2009; Zheng et al. 2018). Six variables are considered as covariates in total and three nonstationary models are developed in terms of different combinations of covariates. To verify the reliability of the NSRI, it is finally compared with the traditional SRI according to the historical droughts.

The outline of this paper is as follows. The study area, data description, and the methodology of SRI and NSRI by using GAMLSS are covered in the section ‘Materials and methods’. The modelling results and the choice of different combinations of covariates are discussed in the section ‘Results and discussion’, followed by the conclusions in the section ‘Conclusion’.

The Heihe River originates from Qilian County in the northeastern part of Qinghai Province. The river basin, between 97°37′ E–102°06′ E and 37°44′ N–42°40′ N, is located in the middle of Hexi Corridor in the northwest of China, with an area of 13,40,000 km². The circulation of the westerly belt in the middle and high latitudes and the polar air mass make the climate across the basin dry, strong wind and sunshine, large temperature differences, less precipitation and great evaporation, resulting in the scarcity of water resources in this basin (Gao & Zhao 2010).

The basin is divided into three reaches, the upper reach, the middle reach and the lower reach. The annual precipitation for the upper reach is about 350 mm, with 70% occurring from June to October. The water resources over the basin are mainly generated from the upper reach and consumed in the middle and lower reaches. The water resources over the basin is 2.62 billion m³, supporting for 1.5 million people and 384,000 ha of farmland irrigation (Ling 2014). The amounts and variations of water resources in the upper reach have great effects on agriculture, economic development and ecological environment for the middle and lower reaches. Thus, it is particularly necessary to carry on the researches about hydrological droughts in the upper reach.

The monthly runoff data at three hydrological stations (Qilian, Zhamushike and Yingluoxia) and the meteorological data at two national meteorological stations (Qilian and Yeniugou) in the upper reach are used. The datasets span from 1961 to 2014. Runoff data are used to calculate the hydrological drought index, whereas meteorological data, including precipitation, temperature, relative humidity and wind speed, are used as covariates to establish the NSRI. The meteorological data are also used for calculating the potential evapotranspiration by employing the Penman–Monteith method (Allen et al. 1998), which is suggested by the Food and Agriculture Organization of the United Nations. The locations of the stations are shown in Figure 1 and the detailed information are displayed in Table 1. These datasets were obtained from the China Meteorological Data Network (http://data.cma.cn/) and the Hydrology Yearbook.

Shukla & Wood (2008) proposed the SRI, analogous to the SPI (McKee et al. 1993). Many studies have shown the efficacy of SRI in depicting the hydrological droughts (e.g., Luciano et al. 2012; Jiang et al. 2019). To obtain this index, fit the runoff series for a certain period by using Gamma distribution firstly and then transform it into a standard normal distribution through an equal probability transformation. The specific calculation process is as follows:

1. Let represents the runoff data at the th month. Fifty-four years (1961–2014) of runoff data are available here, thus n ranges from 1 to 648. For a given time scale of k months (a time scale of 12-month is considered here), the cumulative runoff is calculated as:

   

   (1)
2. Fitting a two-parameter Gamma distribution, denoted as ∼ Gamma (μ, σ):

   

   (2)

   where μ and σ are the location and scale parameters in Gamma, μ > 0 and σ > 0. is the mathematical Gamma function.
3. The cumulative probability for a given time scale can be calculated as:

   

   (3)
4. Converting cumulative probability to a standard normal distribution function:

   

   (4)

When modeling time series over 30 years, we should pay attention to the nonstationarity of the time series which is potentially caused by climate change (Russo et al. 2013). The NSRI is then defined based on a nonstationary Gamma distribution with its location and scale parameters varying with covariates by using the GAMLSS framework (Rigby & Stasinopoulos 2005; Stasinopoulos & Rigby 2007).

The RS algorithm is used to fit the nonstationary Gamma model in the GAMLSS framework (Rigby & Stasinopoulos 2001, 2005). After fitting the nonstationary Gamma model, the remaining calculation process of NSRI is the same as SRI, converting the corresponding cumulative probability of nonstationary Gamma model to a standard normal distribution function. To prevent the model overfitting, Akaike information criterion (AIC) (Akaike 1974), together with the worm plot, is used to test its goodness of fit.

Since the NSRI proposed here is normalized in the same way as the traditional SRI, the same drought-level standards for both indices are suggested (Guttman 1999). As shown in Table 2, positive NSRI values indicate wet conditions, and negative values indicate dry conditions. The higher the NSRI value, the wetter it is, and the lower the value, the drier it is.

Since it is the cumulative runoff series fitted by a distribution, the stationarity of cumulative runoff series is to be tested. The augmented Dickey–Fuller (ADF) test, an augmented version of the original Dickey–Fuller test (Dickey & Fuller 1979), is then performed to the 12-month cumulative runoff series from 1961 to 2014 to verify their stationarities. The computed Dickey–Fuller statistic is compared with the corresponding critical value at a certain P-value. For example, when P = 0.05, the critical value of the ADF test result is −3.5 (with the sample size of 50). If the computed Dickey–Fuller statistic is greater than this critical value, it indicates that the null hypothesis of the ADF test is accepted at the 0.05 significance level, that is, the sequence to be analyzed is nonstationary. In this study, we used the package of ‘tseries’ in r to do the ADF test which can output the P-values and the Dickey–Fuller statistic values directly. The P-value of this test in this package is obtained through interpolating from the critical value table of Banerjee et al. (1993). If the computed P-value is larger than a given significance level, such as 0.05, the null hypothesis is accepted. Otherwise, it is rejected.

As shown in Table 3, the Dickey–Fuller statistics at the QL station are less than −3.5, and the corresponding P-values are less than 0.05, which means that, at the 0.05 significance level, the cumulative runoff series at the QL station exhibit a stationary state. While at ZMSK and YLX stations, the results are opposite. Therefore, we conclude that the cumulative runoff at these two stations can be considered as nonstationary series. To model such series, the NSRI needs to be developed.

The selection of a suitable covariate is of great importance in developing a NSRI. To test whether the runoff and the alternative variables over the study area are related to each other, a correlation test is carried out by using Spearman, Pearson and Kendall correlation tests. The test results are shown in Table 4, in which r indicates the correlation coefficient, sig. indicates the significance level. If sig. <0.05, there is a significant correlation between the two variables, and the higher the absolute value of r, the stronger the correlation. Among the six variables, five of them (precipitation, temperature, relative humidity, potential evapotranspiration and antecedent runoff) are strongly correlated with runoff with correlation coefficients greater than 0.5, and wind speed shows weak correlations. Therefore, except for wind speed, all variables are selected as covariates to develop the NSRI.

For better understating the interactions of different factors and finding the main influencing factors of runoff variation in different seasons, four different types of models are developed with adding different variables to each nonstationary cases, Model 0 (M0), Model 1 (M1), Model 2 (M2) and Model 3 (M3). M0 is the stationary case with constant Gamma parameters, and the remaining three are all nonstationary cases. M1, M2 and M3 take one, two and three variables as covariates, respectively. Different combinations of covariates are shown in Table 5 for details, with 16 different models developed in total. M11–M15 represent the nonstationary models with one covariate, M21–M27 represent the nonstationary models with two covariates and M31–M34 represent the nonstationary models with three covariates. The AIC values corresponding to each model are also displayed in the table. The model with the smallest AIC value is regarded as the best model.

As shown, the stationary model, M0 has larger AIC values, presenting the poor performance in the fitting. Comparatively, the nonstationary models, M1–M3, with assuming that location and scale parameters in Gamma change with the covariates, perform better in terms of the lower AIC values. This is also largely supported by the result of the stationary test. That is, the nonstationary time series show poor performance with stationary models and better performance with nonstationary models.

Comparing the results of different numbers of covariates, it can be found that the AIC values for M2 (with two covariates) show a little lower than those for M1 (with one covariate), also lower than those for M3 (with three covariates) (Table 5). Thus, M2 are selected for our analysis finally. Specifically, M21 (precipitation and temperature as covariates) and M22 (precipitation and relative humidity as covariates) are selected for ZMSK and YLX stations, respectively, according to the lowest AIC values.

To further evaluate the quality of fit, taking August runoff series at the YLX station as an example, the worm plots of M0 and the best nonstationary M22 are provided in Figure 2. The worm plot was first used by van Buuren (2001) to determine the model fit based on whether the residual point is within the acceptable area surrounded by the elliptic curve. According to Figure 2(a), it can be observed that although most residual points are within the acceptable area, minority points are at 95% confidence level boundary and even outside the acceptable area, which means that M0 does not meet the reliability requirements quite well. While, in Figure 2(b), the M22 residual points are all in the acceptable area, indicating that this model satisfies the reliability and fits well. The performance of the nonstationary model with meteorological elements as covariates was thus confirmed for August runoff series at the YLX station.

Figure 3, also taking the cumulative runoff series in August at the YLX station as an example, shows the five estimated percentile results (5th, 25th, 50th, 75th and 95th) from M0 and M22, as well as the observed runoff series. As shown in Figure 3(a), the quantile grayscale image is horizontally linear, and the dynamic change of runoff series under climate change conditions cannot be described well. Figure 3(b) reflects that the percentiles are constantly changing with the observed data, and M22 with two meteorological factors as covariates shows better performance in describing the nonlinear variation in the runoff data. It also implies that the influence of meteorological factors to runoff sequence cannot be neglected.

Overall, in the context of climate change, the nonstationary Gamma distribution with considering meteorological factors as covariates tends to more robust in describing the variation behavior exhibited in the runoff data in our study case.

Figure 4 shows the time series of drought index calculated from the observed runoff data at three hydrological stations. For the QL station, only stationary SRI sequences are calculated from 1961 to 2014. This study defines a drought event as a period of time when the index value is continuously negative and reaches a threshold of −1.0 or lower. In this respect, the dry periods for the QL station were from July 1970 to September 1971, from November 1972 to August 1974, from June 1978 to June 1980, from September 1982 to May 1983, and from May 1992 to August 1992. For ZMSK and YLX stations, both SRI and NSRI sequences are demonstrated in the figure. Overall, the SRI and NSRI series at both stations provide close drought index in most cases, while also defining several differences.

For the ZMSK station, the period of October 1990 to January 1991 is detected as normal state from the SRI, but moderate drought from the NSRI, the period of June to July 2011 is detected as wet from the SRI, while drought from the NSRI. For the YLX station, both SRI and NSRI recorded extreme drought events from August 1973 to July 1974. The results obtained are consistent with historical records and related literature (Yang et al. 2017). However, several significant differences between SRI and NSRI estimates can also be observed (Figure 4). For example, from August 1984 to May 1985, August 1999 to June 2000, and September 2010 to June 2011, SRI detects normal, while NSRI detects them as mild droughts.

The occurrence frequency of different drought grades is calculated from both the SRI and the NSRI. As shown in Figure 5, compared with the results of the SRI, moderate and severe droughts occur more frequently, while mild and extreme droughts are less frequent according to the NSRI for the ZMSK station. For the YLX station, the frequency of mild drought is higher, and the frequencies of moderate and severe droughts are slightly lower in terms of the NSRI.

It also can be found that, in Figure 5, the two stations present completely opposite results in the frequency of drought occurrence, especially at mild and moderate grades. For the ZMSK station, the frequency of mild drought decreased and that of moderate drought increased, while for the YLX station, the results are the opposite. This can be explained from the locations of the two stations. As displayed in Figure 1, the ZMSK station is located upstream of the YLX station. The moderate drought in the upstream (where the ZMSK station is located) may merely cause the weaker drought (mild drought) in the downstream (where the YLX station is located), since there are other tributaries flowing into the main stream and then alleviating the drought grade in the downstream to some extent.

The run-length theory can be used to extract drought charateristics, including the peak, duration and severity (Yevjevich 1967). Drought peak is defined as the minimum of the drought index during a drought event. Drought duration refers to the duration of drought from the start to the end of a drought event. Drought severity refers to the cumulative sum of the differences between the drought index and the threshold during the drought (Yevjevich 1967; Gu et al. 2019). A threshold of −1.0 was set to identify the occurrence of drought events in this study. Drought peak, duration and severity were then extracted from the drought events according to the run-length theory. As shown in Table 6, the QL station identifies 11 drought events, of which the drought peak reaches −2.92, the longest drought duration lasts for 25 months and the maximum drought severity is −59.83. The ZMSK station identifies more drought events (18 events) based on the NSRI, with higher drought peak (−2.42) and higher drought severity (−33.62), compared with the results of SRI (12 events, −2.29, −26.69). The YLX station also identifies more drought events (19 events) based on the NSRI, while the drought peak (−2.47) and the maximum drought severity (−37.66) are lower than those from the SRI (−2.63 and −43.73). The longest drought durations are quite close from both the NSRI and the SRI for this station.

Many meteorological factors show strong effects on droughts under the global warming in recent decades (e.g., Núñez et al. 2014; Zarch et al. 2015), so many of these factors are taken into account in the study of hydrological droughts. Here, six variables are considered, and three nonstationary models are developed. When only one variable is considered, M1, in which the variable of precipitation is taken as a covariate, shows the best for 9 of 24 series, and the series are concentrated basically in April, May, August, September and October (Table 5). It indicates that precipitation is the main factor affecting runoff changes during these months. It is understandable since these months are usually with more precipitation in the study area. M15, taking the variable of antecedent runoff as a covariate, shows the best for 13 of 24 series, and most of them are concentrated from November to about next April. This finding keeps consistent with the river recharge resources in different seasons. In rainy seasons (e.g., summer), precipitation contributes the runoff most, while in dry seasons (e.g., winter and spring), the river is mostly recharged by glacial and snow melting water, permafrost melting water and groundwater over the study area (Li et al. 2006; Wang et al. 2011). Pei et al. (2008) and Wang et al. (2011) concluded that about 75% of runoff are recharged by precipitation in rainy seasons and more than 75% are recharged by snow melting in dry seasons in the Heihe River basin.

In the case of taking two variables as covariates (M2 models), all the corresponding AIC values show slightly lower than those for M1 models (Table 5). It indicates that two variables are better in describing the nonstationary characteristics of the runoff data. In the case of taking three variables as covariates (M3 models), however, the AIC values have gone up instead. This may be the result of over-fitting caused by the over-complexity of the model. Generally speaking, when the number of covariates increases, the complexity of the model increases, the likelihood function also increases, it leads to AIC value decreases. Whereas when the number is too large, the growth rate of likelihood function slows down, which thus results in the increase of AIC (Akaike 1974). In addition, as stated by Agilan & Umamahesh (2018), the increase of covariates would bring greater uncertainty in the final assessment. Thus, M2, with two variables as covariates, are finally selected for our analysis.

Note that many studies also included the time as a covariate in a nonstationary analysis (Russo et al. 2013; Wang et al. 2015a; Javad & Somayeh 2018), although some studies have proved that time is not suitable in their cases due to its linear monotonic trends (Li et al. 2015). It is considered as an additional covariate for comparison to discuss the nonstationary runoff process in this study. As illustrated in Table 7, the corresponding AIC value is moderate, slightly higher than those from the selected models (M21 for the ZMSK station and M22 for the YLX station) and lower than those from the stationary models (M0). Besides, since using time as an explanatory variable has no physical meaning, it can only describe the general trend of a hydrological sequence over time (Vu & Mishra 2019). Thus, the time is not considered as a covariate in our nonstationary analysis.

To discuss the effects of different covariates to the drought index, we compare the drought index calculated from three types and 16 nonstationary models (five models in M1, seven in M2 and four in M3). It is observed that in Figure 6(a) and 6(c), in the case of taking the variable of precipitation as the covariate, no matter what kinds of nonstationary models are (M21, M22, M23 and M24 with two covariates or M31, M32 and M33 with three covariates), the derived drought index are quite close with those from M11 (in which only the variable of precipitation is selected as the covariate). That is, as long as the combinations of precipitation and other variables are taken as covariates for the nonstationary model, the derived drought index curves are quite close. It indicates that the hydrological drought in this study is mainly affected by precipitation, and the contribution of precipitation to hydrological droughts is great over the study area. This conclusion is also confirmed by Li et al. (2006) and Wang et al. (2011, 2015b). When other variables, rather than precipitation, are considered as covariates (M12, M13, M14, M15, M25, M26, M27 and M34), the derived drought index curves fluctuate a great deal, and large differences are found between them (Figure 6(b) and 6(d)). This indicates that the results of drought assessment differ greatly without considering precipitation. All these mean that, although other variables have certain effects to hydrological droughts, precipitation is still the major influence factor in hydrological drought in our case.

According to Figure 4, the SRI and the NSRI at ZMSK and YLX stations display the close drought curves in most cases, while there are still several differences. These differences are attempted to be explained by the meteorological variables. According to Jiang et al. (2014), the hydrological drought in the upper reach of Heihe River basin lag behind the meteorological drought by about 4 months. Therefore, we deduce that the hydrological drought situation of a certain period is possibly related with the precipitation that slides forward for 4 months.

For the ZMSK station, the period of October 1990 to January 1991 is detected as normal state from the SRI, while moderate drought from the NSRI. It is found that the precipitation sliding forward for 4 months (from June to September in 1990) was 328.7 mm, lower than the multi-year mean value of 333.7 mm for the same period. Less precipitation is more likely to cause drought, and thus the conclusion drawn from the NSRI (drought) seems to be more reliable. For the period of June to July 2011, it is found that the observed precipitation from February to March 2011 (6.4 mm) is less than the multi-year mean value (11.0 mm), and the mean temperature (9.04 °C) shows higher than the multi-year mean value (8.28 °C). Owing to the less precipitation and the higher temperature during this period, it is concluded that drought is more likely to occur, and thus the detection from the NSRI (drought) is more reasonable than that from the SRI (wet).

For the YLX station, the drought conditions for periods of August 1984 to May 1985, August 1999 to June 2000 and September 2010 to June 2011 are different from the SRI and the NSRI. It is found that the precipitation sliding forward for 4 months in these three periods (336.8, 327.8 and 333.6 mm) are less than their multi-year mean value for the same period (364.9, 358.3 and 345.5 mm), and the average relative humidity of these three periods (54.50, 53.48 and 51.76%) are also lower than the multi-year mean value (55.47, 54.27 and 52.76%). Therefore, the occurrence of drought (drawn from the NSRI), rather than the normal state (drawn from SRI), is regarded as more reasonable.

Another observation is that the same runoff value corresponds the same drought classification based on the SRI, which is because its calculation solely depends on runoff data (Shukla & Wood 2008), and the drought state detected from the SRI keeps stable. However, this conclusion is not applicable to the nonstationary model. The NSRI leads to differences in drought classification with the same runoff value. For example, a monthly runoff of 150 mm for the YLX station corresponds to a mild drought state in terms of the SRI, but corresponds to mild drought (May 1964), normal state (June 1975) and moderate drought (December 2004) in terms of the NSRI. A monthly runoff of 183 mm is classified as normal state from the SRI, while classified as moderate wet (February 1982), normal (August 1999) and mild drought (September 2010) states from the NSRI.

In a changing environment, hydrological drought can no longer be evaluated solely on runoff. Zhang et al. (2018) investigaed the impacts of both climate variations and human activities to hydrological droughts in the middle reach of Yangtze River in China. Wang et al. (2019b) incorporated the climate-driven and human-induced variables as covariates in the nonstationary model to indicate the nonstationary features of runoff series. In this study, the determination of NSRI depends not only on runoff data but also on the key factors affecting runoff. Variables including precipitation and temperature are considered as covariates in characterizing the nonstationary features in the runoff data, and thus drought state detected by the NSRI varies with the influence of covariates.

We also calculated the NSRI with time as a covariate (denoted as the TSRI here) for comparison (Figure 7). It is observed that TSRI differs a great deal with SRI especially for the first and last several years. This is because the location parameter is constant and close to its mean value in the SRI, while it is time-dependent in the TSRI. This condition leads to more differences between TSRI and SRI for the early and last years of the record period and less differences for the middle part of the period (Javad & Somayeh 2018). Furthermore, the linear monotonic trend of time also makes it not as appropriate as other meteorological variables in characterizing the nonstationary processes of runoff series (Li et al. 2015).

These results support that, compared with the traditional SRI, the new developed NSRI, with considering meteorological variables as covariates, is proved to be more robust and more applicable in this drought assessment.

Appropriate drought indices are often important tools for regional drought assessments. In the context of changing environment, the traditional SRI which is defined based on a stationary Gamma distribution is no longer suitable due to the nonstationary characteristics of runoff series (Villarini et al. 2009; Russo et al. 2013). Thus, the NSRI is developed in this study within the GAMLSS framework. The proposed new drought index considers climate change and incorporates six variables into calculation, accounting for the nonstationary features of runoff data. As a quantitative indicator, just like the SRI, the NSRI can assess drought events at different time scales (e.g., 3, 6, 9 and 12 months) and allow comparison of climatic conditions at different locations. The 12-month time-scale NSRI is applied to the upper reach of the Heihe River basin.

The conclusions can be summarized as follows. For the study area, (1) the fitting results of the nonstationary models to those nonstationary runoff data (ZMSK and YLX stations) are much better than the stationary models according to AIC values; (2) for the nonstationary models with one covariate, precipitation shows better in the fitting as a covariate in rainy seasons, and antecedent runoff shows better in dry seasons, which is highly related with the recharge sources in different seasons; the nonstationary models with two covariates perform better than those with one or three covariates; (3) the NSRI identifies more drought events than the SRI, and precipitation is the major factor affecting the hydrological drought conditions over the study area.

In general, the current study solves the problems of hydrological drought assessment based on nonstationary runoff data. It is proved that the proposed new drought index, NSRI, is capable of providing more reliable results for drought assessment over the study area. Revelant findings are also very helpful in developing strategies for coping with local droughts and water resource risk management. It should be noted that, although the findings here are in the context of NSRI at 12-month scale and limited stations, the framework and the methods in this study could be extended over further scales and regions.

This paper mainly considers the relation of several meteorological elements to hydrological droughts. In fact, hydrological droughts are also induced by factors like the underlying surface of the basin, vegetation coverage and human activities (such as reservoir opening and irrigation) (Shukla & Wood 2008; Fang et al. 2018b; Zhang et al. 2018; Wang et al. 2019a, 2019b, 2019c). The use of appropriate quantitative factors to characterize human activities will be considered in the future research for making the nonstationary drought index more reliable under the complex conditions. Additionally, to identify the hydrological drought propagation and its key influencing factors is also a hot topic and can be included in the future studies.

This study is supported by the Fundamental Research Funds for the Central Universities (No. 39435832015028).

The meteorological dataset is available from the China Meteorological Data Network (http://data.cma.cn/). The hydrological dataset is available from the Hydrological Yearbook of the People's Republic of China.