Most current drought indices rely on a representative parametric distribution function to fit data, which results in different tail behaviors. Additionally, a drought index based on a single variable may not be sufficient for monitoring drought conditions timely and reliably. Therefore, a nonparametric multivariate drought index (NMSDI) combined with the information of precipitation and streamflow was introduced in this study, without assuming representative parametric distributions. It was applied to characterize drought in the Yellow River Basin (YRB) on seasonal and annual scales. Results indicate that: (1) the variations of developed NMSDI is well consistent with those of 1-month SPI (standardized precipitation index) and SSI (standardized streamflow index), indicating the reliability and effectiveness of the newly proposed nonparametric based drought index; (2) a decreasing NMSDI trend was found over the period of 1952–2012 at seasonal and annual time scales, which would reverse in the future as suggested by the Hurst index; (3) no significant change points were detected for the annual NMSDI series over the YRB except one (i.e. the year 1991) in the middle streamflow sub-basin Wei River Basin (WRB) which was primarily caused by the combined effects of climate change and human activities.

## INTRODUCTION

Drought is a natural recurring climate extreme that usually occurs over long time spans and spreads across a large spatial scale (National Research Council 2010). Compared with other natural hazards like flood and hurricane, drought is one of the costliest and least understood hazards, exerting huge adverse influences on ecological and social systems (Wilhite 2000; Hao & Aghakouchak 2014). As the Federal Emergency Management Agency (Federal Emergency Management Agency 1995) reported that annual losses caused by droughts in the USA is $6–8 billion, and the 1987–1989 drought occurring across a large portion of the USA led to the losses of about $39.4 billion in direct and indirect ways, which is obviously larger than that induced by hurricanes and floods (WGA 2004). Global warming is expected to result in intensified water cycles, and extreme events have been becoming more frequent, severe and pronounced (e.g. drought and flood) worldwide (Kunkel 2003; Roy & Balling 2004). Therefore, understanding the characteristics of droughts is very important for hazard adaptation and mitigation in a changing environment (Hao & Aghakouchak 2013; Huang *et al.* 2014a, 2014b).

Drought indices are important tools for drought characterization and prediction due to their abilities to simplify the complex interactions among various climate and climate-related variables. They allow scholars to quantitatively assess climate anomalies in terms of their frequency, duration, and severity, as well as spatial extent (Wilhite 2000). During the past decades, various indicators have been proposed to characterize droughts, and each has its own advantages and disadvantages (Mishra & Singh 2010). Some of the frequently used indexes are the standardized precipitation index (SPI; McKee *et al.* 1993, 1995), Palmer drought severity index (PDSI; Palmer 1965), crop moisture index (CMI; Palmer 1968), soil moisture drought index (SMDI; Hollinger *et al.* 1993), vegetation condition index (VCI; Liu & Kogan 1996), etc. Amongst these drought indicators, SPI is one of the most frequently utilized indices for meteorological drought assessment (Shukla *et al.* 2011; AghaKouchak & Nakhjiri 2012) due to its simplicity, standardized feature, and flexibility across various time scales (e.g., 1-, 3-, 6-, 12-month) (Hao & Aghakouchak 2014; Farahmand & Aghakouchak 2015).

Nevertheless, SPI has a potential weakness in that an appropriate parametric probability distribution function (PDF) has to be assumed in fitting precipitation time series (Angelidis *et al.* 2012). In practice, many issues will arise from the assumption that samples should follow a specific distribution. Typically, in constructing the SPI, the precipitation records are fitted to a gamma PDF which, however, may not be the most appropriate selection of distribution (Guttman 1999; Quiring 2009). Quiring (2009) pointed out that the developed SPI values, especially in the tails, were extremely sensitive to the selection of parametric distribution function. In fact, there exists no worldwide-accepted parametric distribution for meteorological and hydrologic variables (Silverman 1986; Smakhtin 2001), and the use of a specific parametric distribution is expected to lead to a large deviation for the low or high quantiles (Sharma 2000). Thus, a nonparametric method, which is distribution-free, is a good alternative to parametric approaches in developing drought indicators.

Furthermore, the existing available indices have some drawbacks including statistical incomparability and temporal inconsistency (Farahmand & Aghakouchak 2015). For instance, since SPI and PDSI have different scales, they cannot be directly compared with each other (Steinemann & Cavalcanti 2006). Additionally, droughts are categorized differently based on different types of shortages (e.g. the deficits of precipitation, streamflow, soil moisture, etc.) (Heim 2002). Univariate drought indicators can only focus on specific aspects of water shortages, leading to varying performances (e.g. in identification of the drought onset and termination). A meteorological drought (deficit in rainfall) may develop and end quickly, while the onset and termination of the corresponding hydrological drought (shortage in runoff) responds to the meteorological drought with some lag time (Pandey & Ramasastri 2001). Therefore, the current consensus amongst scholars is that no indicator based on a single variable is ideal for drought monitoring and prediction (AghaKouchak 2014). A holistic method for drought characterization needs an integrated indicator merging various sources of water shortages (precipitation, runoff, soil moisture, evapotranspiration, etc.). The advantage of standardized indicators is that they provide an opportunity to construct statistically consistent indicators on this basis of precipitation (SPI), runoff (SRI), soil moisture (SSI), etc. (Farahmand & Aghakouchak 2015). Therefore, characterizing droughts based on multiple climate variables using a nonparametric approach is very essential for developing spatially and temporally consistent drought indices, which is the major motivation of this study.

The Yellow River ranks the second largest river in China and the sixth largest river in the world in terms of its length (Shiau *et al.* 2007). In northern China, the Yellow River is a primary source of freshwater for nearly 107 million residents within the watershed, accounting for 8.7% of the total population in China (Wang *et al.* 2006). Meanwhile, the Yellow River Basin has been severely plagued by droughts from ancient times to the present (She & Xia 2013). Historically, one continuous drought hit the Yellow River Basin (YRB) in 1637–1643, and the extremely severe drought contributed to the collapse of the Ming Dynasty (Xie & Fu 2004). Another drought occurring between late 1920s and early 1930s had affected about 20 million people, and it is claimed that more than 3 million people died due to drought-related diseases and famine (Xie & Fu 2004). After the 1970s, the zero-flow phenomena which frequently occurred in the downstream of the Yellow River have gained wide attention. Highly frequent flow interruptions during the past 30 years have led to widespread unfavorable influences on local social and ecological development. Hence, it is of great importance to scientifically and reasonably monitor droughts based on a reliable drought indicator in the YRB.

Therefore, the primary objectives of this study are: (1) to develop a nonparametric multivariate standardized drought index (NMSDI) for comprehensively and effectively characterizing the drought conditions in the YRB; (2) to investigate the drought evolution characteristics including its trend and persistence in the future; (3) to identify possible change points of calculated NMSDI series, thus determining whether drought evolution in the YRB is non-stationary.

## STUDY AREA AND DATA

### The Yellow River Basin (YRB)

^{2}. The YRB lies between 95°E–119°E and 32°N–41°N, as presented in Figure 1. Because Chinese ancestors have lived in this region since prehistoric times, the Yellow River is always called the ‘Mother River of China’ (Fu

*et al.*2004). Note that the Chinese Loess Plateau where water loss and soil erosion is extremely severe is located in the middle of the river basin. Every year, nearly 1,060 million tons of sediment is delivered from the loess plateau to the Bohai Sea (Milliman & Meade 1983). The YRB is impacted by the continental monsoon (Wu

*et al.*2013). The annual precipitation in the YRB ranges from 123 to 1,021 mm and its annual pan evaporation is between 700 and 1,800 mm (Shao

*et al.*2006). The mean annual precipitation is about 378 mm, which increases from the northwest to the southeast of the basin (Wu

*et al.*2013). Because of the different topographies and large area, precipitation distribution shows a noticeable discrepancy. In order to systematically investigate the drought characteristics of the basin, the whole basin was partitioned into eight sub-basins on the basis of the secondary basin boundary in China, and their zone numbers are 26, 28, 29, 31, 33, 36, 40, and 41, respectively (Figure 1).

### Data

The observed gridded monthly precipitation and simulated monthly streamflow data based on the Variable Infiltration Capacity (VIC) model covering 1952–2012 in the YRB were utilized in this study. The VIC model (Liang *et al.* 1994) was applied to derive the land surface hydrologic fluxes/states. It is a semi-distributed macro-scale hydrologic model characterized by representing sub-grid variability of precipitation, soil moisture storage capacity, topography, vegetation classes, etc. (Liang *et al.* 1994; Nijssen *et al.* 1997). Meteorological forcing (e.g. precipitation, wind speed, temperature) for driving the VIC model was acquired from China Meteorological Administration (CMA), which was interpolated into 0.25 degrees. Land surface characteristics like vegetation, soil, and elevation were gained from Nijssen *et al.* (2001). Six parameters, i.e. the infiltration parameter *b*, the second and third soil layer depths (*d*_{2}, *d*_{3}), and the three parameters in baseflow scheme (*D _{m}*,

*D*,

_{s}*W*), were calibrated to match the long-term monthly streamflow observations. The dataset offers a current state-of-art estimate of land surface fluxes/states, and it is of great value in assessing the long-term water balance. Model calibration/validation was conducted for each river basin (including the Yellow River Basin) in China. In total, meteorological forcing from about 756 meteorological stations were interpolated into 0.25 degrees as inputs into the model, and streamflow records from 15 hydrological stations were used for model calibration/validations. An optimization algorithm of the multi-objective complex evolution of the University of Arizona (MOCOM-UA) (Yapo

_{s}*et al.*1998) was used to find the optimal parameter set. The Nash–Sutcliffe efficiency and the relative error between the simulated and observed mean annual streamflow were used to assess the model performance. Details for the model calibration/validation should be referred to Zhang

*et al.*(2014).

## METHODOLOGIES

### Nonparametric multivariate standardized drought index (NMSDI)

*et al.*2014b). The calculation of SSI is similar to SPI (Li

*et al.*2013). Since parametric methods have some limitations in developing drought indicators, the empirical probability which is distribution free can replace them to derive a nonparametric standardized index. The empirical Gringorten plotting position (Gringorten 1963) was adopted to calculate the marginal probability of precipitation (or other variables). where

*n*is the sample size,

*i*is the rank of non-zero precipitation series from the smallest, and

*P*(x

*) is the corresponding empirical probability. The outputs of Equation (1) can be turned into a Standardized Index (SI) as follows: where is the inverse of the standard normal distribution function, and*

_{i}*p*is the probability acquired from Equation (1).

*X*is precipitation and

*Y*is streamflow), their bivariate distribution is defined as: where

*p*is the joint probability of

_{j}*X*and

*Y*(e.g. precipitation and streamflow).

*et al.*(1999): where

*m*is the number of occurrences of the pair (

_{k}*x*,

_{i}*y*) satisfying the conditions and , and

_{i}*n*represents the sample size. Being similar to univariate drought indices, the joint probability can be standardized through Equation (2) to obtain a NMSDI value. It should be noted that the effects of natural and anthropogenic factors governing the changes in streamflow would propagate into the derived NMSDI values, thus complicating the detection of droughts. Since 1-month SPI and SSFI were used to develop NMSDI in this study, the time scale of the NMSDI obtained in this study is 1 month.

### The modified Mann-Kendall (MMK) trend test method

The initial Mann–Kendall (MK) trend test method recommended by the World Meteorological Organization (Mitchell *et al.* 1966) is a nonparametric approach. However, the MK test results are expected to be impacted by the persistence of hydro-meteorological series. Hamed & Rao (1998) provided a modified Mann–Kendall (MMK) trend test via taking into account the lag-*i* autocorrelation to remove the persistence. Daufresne *et al.* (2009) pointed out the MMK test approach is more robust than the initial MK in capturing the trends of hydro-meteorological series. The detailed procedures of the MMK can refer to Huang *et al.* (2014a).

### The Rescaled Range (R/S) analysis

The R/S analysis was first developed by Hurst (1956) and applied to investigate the hydrological changes in the Nile River. Hereafter, R/S analysis has been extensively used in morphology, hydrology, cell reproduction, earthquake activity, etc. (Hosking 1984; Wang *et al.* 2013). The basic idea of the R/S analysis is to change the timescale of a specific series and explore its statistical features at different time scales (Oliver & Ballester 1996).

*C*is constant, and

*H*is the Hurst index.

Different Hurst indexes represent different types of time series. For *H* = 0.5, it means this series is totally independent; for , it indicates the future trend of the series may be opposite to the existing series, and a smaller *H* value exhibits a stronger opposite persistence (Oliver & Ballester 1996). For , it denotes the future trend of the series may be consistent with the existing series, and a larger *H* means a stronger persistence. Therefore, this method was adopted to investigate the persistence of drought in the future in the YRB.

### The heuristic segmentation method

Existing statistical test methods like the Mann–Kendell test, sliding *T* test, sliding *F* test, and rank sum test are frequently utilized to identify change points. All of them strongly rely on the assumption that input data should be stationary. However, hydrologic series present greatly nonlinear features due to the high variability of the hydrologic process. Thus, these approaches are difficult to capture real change points. The heuristic segmentation method proposed by Pedro *et al.* (2001) is based on the idea of sliding *T* test and is modified to divide a non-stationary series into several stationary series, thus overcoming the limitation of the above methods.

*et al.*2001). The averages of the left and right subsets of the pointer expressed by

*μ*

_{1}and

*μ*

_{2}, respectively, are calculated. For two Gaussian-distributed random series, the difference between

*μ*

_{1}and

*μ*

_{2}under the statistical significance is assessed by Student's

*t*-test statistic as follows: where is the pooled variance,

*s*

_{l}and

*s*

_{2}denote the standard deviations of the two series, respectively,

*N*

_{1}and

*N*

_{2}represent the number of points from the two series, respectively. Through moving the pointer along the given time series point-by-point, the statistic

*t*is calculated to evaluate the difference between the averages of the right side and left side time series. A larger

*t*means that the average values of the two time series tend to be more different. The largest

*t*value is seen as a good candidate for the cut point. Then, the statistical significance

*P*(

*t*

_{max}) is computed. Note that

*P*(

*t*

_{max}) is not the standard Student's

*t*-test because the series are not independent and cannot be obtained in a closed analytical form; therefore,

*P*(

*t*

_{max}) is approximately calculated as follows: where , and which was acquired from the Monte Carlo simulations,

*N*is the number of the time series to be cut,

*v*=

*N*–2, and

*I*(

_{x}*a*,

*b*) is the incomplete beta function. If the difference between the averages is not statistically significant (e.g.

*P*(

*t*

_{max}) is less than a threshold of 0.95), the time series will be not divided. Conversely, the time series is divided into two segments. If the time series is divided, the iteration of the same procedures on each new segment will continue until the acquired significant value is smaller than the threshold or the length of the acquired segments is less than the presupposed minimum segment length .

## RESULTS AND DISCUSSION

### The performance of the NMSDI

Moreover, in order to further verify the reliability of the calculated NMSDI, the historical severe droughts which occurred in the middle of the YRB were collected from census data to further confirm the performance of NMSDI. As shown in Figure 2(b), the occurrence time of droughts matched well with the low NMSDI values. Two typical drought events which occurred in 2005 and 2006 were picked out (YRCC 2006, 2007), which are perfectly reflected by NMSDI exhibited in the area covered by the large black rectangle in Figure 2(b). The consistence between the timing of low NMSDI values and the years when historical drought events occurred demonstrates the reliability of the developed NMSDI in characterizing drought.

The comparison between NMSDI, SPI, and SSI in other sub-basins in the entire period also exhibits the same characteristics with sub-basin 26. For brevity, the related figures are omitted. The correlation coefficients between monthly NMSDI and SPI/SSFI spanning 1953–2012 in the eight subzones are shown in Table 1. It can be found that the correlation coefficients between monthly NMSDI and SPI/SSFI are very high, and all of them are significant at the 99% confident level. The results further verify the reliability and validity of the developed NMSDI. In view of the comprehensiveness and effectiveness of NMSDI in characterizing drought, it was applied to study the drought characteristics of the YRB in this study.

Subzones . | SPI . | SSFI . |
---|---|---|

26 | 0.96 | 0.91 |

28 | 0.93 | 0.93 |

29 | 0.91 | 0.92 |

31 | 0.89 | 0.88 |

33 | 0.92 | 0.92 |

36 | 0.95 | 0.92 |

40 | 0.88 | 0.93 |

41 | 0.95 | 0.91 |

Subzones . | SPI . | SSFI . |
---|---|---|

26 | 0.96 | 0.91 |

28 | 0.93 | 0.93 |

29 | 0.91 | 0.92 |

31 | 0.89 | 0.88 |

33 | 0.92 | 0.92 |

36 | 0.95 | 0.92 |

40 | 0.88 | 0.93 |

41 | 0.95 | 0.91 |

### The spatial characteristic of NMSDI in the YRB

*et al.*2013), which further verifies the validity of the NMSDI. Generally, the annual NMSDI in the basin has an obvious difference. The upstream and northern basins are drier than the other portions of the basin, which is consistent with She & Xia (2013).

### The temporal change in NMSDI

### The identification of change points of annual NMSDI

*P*

_{0}was chosen as 0.95 and was selected as 25, and the identification results of change points in sub-basin 26 and 31 are exhibited in Figure 8. Figure 8(a) indicates no change point was detected in sub-basin 26 due to the probability of its largest

*T*smaller than the threshold (

*P*

_{0}) during its first iteration, whilst Figure 8(b) indicates that one change point was found in sub-basin 31 because of the probability of its largest

*T*larger than the threshold (

*P*

_{0}) during its first iteration. The identification of change points in other sub-basins was also performed. The identification results show that only the WRB (sub-basin 31) had a change point (1991), indicating that the annual NMSDI series in the YRB except for the WRB remained stable despite the variations caused by the joint impacts of climate change and human activities.

The average annual precipitation in the WRB decreased by 11.2% from 625 mm (before 1991) to 555 mm (after 1991). The decrease of precipitation availability has a certain contribution to the change point in the WRB. In addition, Huang *et al.* (2014c) pointed out that the increasing air temperature and the initially decreasing then increasing potential evaporation had a certain effect on the change points of the relationship between rainfall and runoff in the WRB. Since the NMSDI is calculated based on precipitation and runoff, the variations of temperature and potential evaporation may partly lead to the change point of annual NMSDI series in the WRB. However, the quantitative impact of temperature and potential evaporation on the change point of annual NMSDI series need to be further investigated in future studies.

In addition to climate change, anthropogenic interferences may play an important role in leading to the change point. Soil loss is greatly severe in the WRB. Hence, the soil conservation practices including grass-planting, afforestation, building check dams, and creation of level terraces have been performed since 1960. The detailed areas of soil conservation practices are shown in Table 2.

. | Area of soil conservation practices (km^{2}). | ||||
---|---|---|---|---|---|

Time . | Level terrance . | Afforestation . | Grass-planting . | Check dam . | Total . |

1960 | 172 | 327 | 47 | 6 | 552 |

1970 | 974 | 1,309 | 193 | 24 | 2,500 |

1980 | 2,918 | 3,886 | 446 | 64 | 7,314 |

1990 | 4,758 | 8,468 | 2,320 | 78 | 15,624 |

2000 | 9,088 | 14,924 | 3,648 | 134 | 27,794 |

2006 | 11,779 | 17,157 | 4,265 | 143 | 33,344 |

. | Area of soil conservation practices (km^{2}). | ||||
---|---|---|---|---|---|

Time . | Level terrance . | Afforestation . | Grass-planting . | Check dam . | Total . |

1960 | 172 | 327 | 47 | 6 | 552 |

1970 | 974 | 1,309 | 193 | 24 | 2,500 |

1980 | 2,918 | 3,886 | 446 | 64 | 7,314 |

1990 | 4,758 | 8,468 | 2,320 | 78 | 15,624 |

2000 | 9,088 | 14,924 | 3,648 | 134 | 27,794 |

2006 | 11,779 | 17,157 | 4,265 | 143 | 33,344 |

The area of soil conservation is increasing, and the growth rate has rapidly increased after 1990, indicating that local anthropogenic interferences on underlying surface have greatly intensified since 1990. The soil conservation practices alter local microtopography, intercept precipitation, facilitate infiltration rate of water flow, and retain or slow down streamflow, thus reducing streamflow (Chang *et al.* 2015). Furthermore, with the deepening of China's reform and opening up, a rapid economic development has been achieved since 1990. The mean economic growth speed in the WRB is about 10.5%. The high-speed economy requires numerous water resources to support its sustainable development. The average industrial and domestic water consumption in the WRB after 1990 is about 4.3 billion m^{3}/year, which increased by 52.6% compared with that before 1990 (Huang & Fan 2013), directly resulting in the decrease of runoff in the WRB. Therefore, the coupled impacts of climate change and human activities contribute to the change point in the WRB.

## DISCUSSION

It is well demonstrated in the literature that future drought will become more frequent, severe, and pronounced in most places with global warming. However, the declining drought tendency as shown in this study are based on the Hurst index of the R/S analysis which is conducted by changing the timescale of a specific series and exploring its statistical features at different time scales (Oliver & Ballester 1996). Therefore, the calculated persistence is a projection of the future possible conditions in a probabilistic way in nature. Although the R/S analysis has been extensively applied in morphology, hydrology, cell reproduction, earthquake activity, etc. (Hosking 1984; Wang *et al.* 2013), the persistence based on this method cannot be used directly for future drought projections. Instead, to better project the future droughts, more physically based approaches such as the global climate model (GCM) should be adopted.

## CONCLUSIONS

A reliable drought indicator is important for drought mitigation and water resources management. This study introduces a nonparametric multivariate drought index (NMSDI) coupled with the information of precipitation and streamflow, which is distribution free and can eliminate the limitation of parametric approaches. The integrated drought index was used to characterize drought in the YRB on seasonal and annual scales. The primary results are as follows:

The variations of 1-month NMSDI is well consistent with that of 1-month SPI and SSI, indicating the reliability and effectiveness of the nonparametric method in developing an integrated drought index. Generally, NMSDI is more sensitive to identify drought onset, persistence and termination.

The upstream and northern portion of the basin are drier than other areas. The seasonal NMSDI pattern is similar to the annual NMSDI pattern to a certain degree. Overall, the drought risk in spring and summer is larger than autumn and winter.

The YRB is dominated by a decreasing trend of NMSDI at seasonal and annual scales, indicating an increasing risk of drought during the past several decades. However, the decreasing trend of NMSDI may be the opposite in the future, which suggests that drought risk may reduce in the basin.

Annual NMSDI series in the YRB is relatively stable, except for the WRB in which a change point (1991) was identified mainly due to the combined impacts of climate change and human activities.

## ACKNOWLEDGEMENTS

This research was supported by the National Major Fundamental Research Program (2011CB403306-2), the National Natural Fund Major Research Plan (51190093) and the Natural Science Foundation of China (51179148, 51179149, and 51309188).