Assessment of drought evolution characteristics based on a nonparametric and trivariate integrated drought index

Journal Pre-proofs Research papers Assessment of drought evolution characteristics based on a nonparametric and trivariate integrated drought index Ying Zhang, Shengzhi Huang, Qiang Huang, Guoyong Leng, Hao Wang, Lu Wang PII: DOI: Reference:

S0022-1694(19)30965-5 https://doi.org/10.1016/j.jhydrol.2019.124230 HYDROL 124230

To appear in:

Journal of Hydrology

Received Date: Revised Date: Accepted Date:

2 July 2019 29 August 2019 11 October 2019

Please cite this article as: Zhang, Y., Huang, S., Huang, Q., Leng, G., Wang, H., Wang, L., Assessment of drought evolution characteristics based on a nonparametric and trivariate integrated drought index, Journal of Hydrology (2019), doi: https://doi.org/10.1016/j.jhydrol.2019.124230

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

Assessment of drought evolution characteristics based on a nonparametric and trivariate integrated drought index

Ying Zhanga, b, Shengzhi Huanga*1, Qiang Huanga, Guoyong Lengc, Hao Wangd, Lu Wanga

a State Key Laboratory of Eco-hydraulics in Northwest Arid Region of China, Xi’an University of Technology, Xi’an 710048, China b Yellow River Engn Consulting Co Ltd, Zhengzhou 450003, Peoples R China. c Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China d State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China

*

Corresponding author at: State Key Laboratory of Eco-hydraulics in Northwest Arid Region of China, Xi’an University of Technology, Xi’an 710048, China. Tel.: +86 29 82312801; fax: +86 29 82312797. E-mail Address: [email protected]. 1

Abstract Drought index is an important tool for drought research. The development of a reliable and simple integrated drought index is a critical step of current drought monitoring and risk assessments. The Gringorten point position formula was adopted to develop a trivariate standardized drought index (TSDI) based on precipitation (meteorology), runoff (hydrology), and soil moisture (agriculture) information in this study. This index can be used to characterize droughts from an integrated meteorological, hydrological, and agricultural point of view, and explore drought evolution characteristics. The possible driving forces of drought were explored by examining the associations between sunspot activities and large-scale atmospheric circulation anomalies. A typical Loess Plateau watershed, the Wei River Basin (WRB), was selected as a case study. Results showed that: 1) the TSDI is better at characterizing real drought conditions than other traditional single-variable indices; although it is similar to the parametric trivariate drought index based on copula function (Serinaldi and Grimaldi., 2007; Yang et al., 2018), it is much simpler to calculate and more suitable for characterizing drought of different regions; 2) the value of integrated drought index gradually increased in all areas, and its change point occurred in the early 1970s and 1990s mainly due to variations in runoff; 3) the evolution of the integrated drought appeared to be strongly influenced by sunspot activities and abnormal atmospheric circulation factors, where sunspots activities are dominant. This study provides a valuable reference for the development of an integrated drought index, which might be applied to local drought monitoring, preparedness and mitigation. 2

Keywords: drought; integrated drought index; nonparametric distribution; trivariate empirical frequency; the Gringorten point position formula

1 Introduction

Drought is a kind of natural disaster with extremely detrimental effects, which occurs frequently, affecting a wide range of areas (Huang et al., 2019; Fang et al., 2019a, b). It not only affects agricultural production, but can also hinder economic development, ultimately leading to social instability. China is drought-prone since ancient times. Recently, it has evolved in a big agricultural and developing country. In this context, drought detection and prevention are undoubtedly becoming priorities for the Chinese government (Huang et al., 2019). The 21st century environment is being threatened by global climate change (Zhao et al., 2019; Han et al., 2019). These alterations are accelerating the global water cycle (Milly et al., 2002) and increasing the frequency of extreme hydrological events (e.g., droughts and floods). Since the last century, air temperature has increased by 0.85 °C (IPCC, 2013) under the influence of global climate changes and human-induced greenhouse gas emission. Consequently, the world has been suffering severe droughts during the last decade, whose frequency is expected to increase in many regions (Huang et al., 2019). Drought events have seriously affected the daily life of local residents and impeded the growth of crops. In addition, they significantly hindered industrial production and development, eventually impacting national economic development and public health security (Sheffield and Wood, 2008; Mussa et al., 3

2015). Hence, it is urgent and necessary to conduct in-depth research on this topic. Drought indices have been widely applied by drought researches, being on the basis of numerous drought investigations and playing an important role in the assessment of drought parameters (Kao et al., 2010; Huang et al., 2016a, b). A large number of drought indices has been developed in order to characterize different types of droughts (Shukla and Wood, 2008; Hao and AghaKouchak, 2013), which provide reliable indices for the realization of drought risk assessments and the formulation of disaster mitigation strategies. Droughts are generally divided into four categories: meteorological droughts (resulting from a precipitation deficit), hydrologic droughts (resulting from insufficient runoff), agricultural droughts (resulting from a shortage in soil moisture), and socio-economic droughts (resulting from insufficient water supply with respect to water demand) (Guo et al., 2019a, b). The four types of droughts present some common characteristics (Huang et al., 2017a). As a matter of fact, one type of drought can evolve into another type (Huang et al., 2017a). For example, meteorological droughts may transform into hydrological or agricultural droughts (Wu et al., 2018a, b). A variety of indices have been developed for these four types of droughts based on different elements, and each with its own advantages and defects. Early drought indices were based on a single variable, like the standardized precipitation index (SPI; McKee et al., 1993), the standardized runoff index (SRI; Shukla and Wood, 2008), and the standardized soil moisture index (SSI; Hao and AghaKouchak, 2013), which have all been extensively used in drought monitoring. Among them, the SPI has been especially applied to identification of early drought 4

occurrences. Other drought indices, like the Palmer drought severity index (PDSI; Palmer, 1965), the self-calibrating PDSI (SC-PDSI; Yi et al. 2017), and the standardized precipitation evapotranspiration index (SPEI; Vicente-Serrano et al. 2010) have been also widely used. Nevertheless, investigations based on a single variable are increasingly considered unreliable, since droughts are influenced by multiple variables (e.g., precipitation, runoff, soil moisture) (Hao and AghaKouchak 2013; Kao and Govindaraju, 2010; Mishra and Singh 2010). In fact, drought events are attributed to multiple elements of water shortages over different periods of time (Yang et al., 2018).

Several integrated drought indices have been proposed with the aim of

investigating drought events in a more integrated manner (Svoboda et al., 2002; Huang et al., 2015). For example, the linear weighted average was adopted to blend several drought indices. The US drought monitor (USDM) developed by Svoboda et al. (2002) can contain information on climate indices, numerical models, and inputs from regional and local experts. The aggregated drought index developed by Keyantash and Dracup (2004) can aggregate all the selected variables into a single time series, through a correlation-based principal component analysis (PCA). However, both of these aggregated indices have some weaknesses. First, pairs of drought variables are characterized by nonlinear relationships. Second, the assumptions of the PCA are not always met in real situations. Overall, the use of these two methods may lead to information distortion. Subsequently, the bivariate copula-based joint index proposed by Kao and Govindaraju (2010) has been applied for the quantification of drought characters 5

based on precipitation and streamflow information. This type of drought indices is multivariate and hence more effective than those previously mentioned. However, this joint distribution method becomes too complicated to allow a rapid and effective calculating process of dependence structures among variables when it involves more than three variables. Based on the concept of copulas, Hao and AghaKouchak (2013) developed the multivariate standardized drought index (MSDI). This index is calculated based on the joint cumulative probability of precipitation and soil moisture content. Thus, it is able to characterize drought conditions based on either of these two parameters. Unfortunately, this index includes only meteorological and agricultural drought information. Among the four main types of droughts, only three of them can be physically described. In order to improve their characterization, researchers should consider not only meteorological and agricultural droughts, but also hydrological drought. Recently, a trivariate drought index has been proposed by Yang et al. (2018) via employing copulas combining meteorological, hydrological, and agricultural variables, which is referred to as nonparametric multivariate standardized drought index (NMDI). Copulas strongly rely on the assumption that data samples follow a particular probability density function (pdf) (Huang et al., 2016a, b; Fang et al., 2019a, b). However, the given pdfs usually fail to reflect real drought conditions exactly (Huang et al., 2016a), leading to deviations from reality. In addition, the use of copulas involves many tedious steps and complex calculations. Therefore, it is necessary to explore simpler ways to obtain a reliable and effective integrated drought 6

index. The empirical frequency plotting position formula proposed by Gringorten (1963) proved to be suitable for the fitting of hydrological time series. Hence, it has been widely applied to hydrological frequency analyses (Hao and Aghakouchak, 2014; Huang et al., 2016a, b). Since it reflects a nonparametric empirical frequency distribution, it can skilfully avoid fitting optimized parametric distribution. This unbiased plotting position formula, which works without following any specific distributions, may solve the problems mentioned above. In brief, the main objectives of this study were to: (1) develop a simple but effective TSDI, verifying its reliability and superiority, (2) analyse the temporal-spatial evolution of the TSDI in the Wei River Basin (WRB), and (3) explore the driving forces of the evolution characteristics of droughts from an integrated perspective.

2 Study area and data

2.1 Study area The WRB, which represents a typical Loess Plateau watershed in an arid/semi-arid region of China, was selected as the study area for this study (Fig. 1). Located between 33.5–37.5°N and 103.5–110.5°E, the Wei River is the largest tributary of the Yellow River (Liu et al., 2019). The altitude of the WRB gradually increases from southeast to northwest. The spatial and temporal distribution of precipitation over this area is very uneven, which is one of the main triggers for the frequent droughts and waterlogging disasters in the region. The WRB is currently one of the most important 7

agricultural (e.g., grain, cotton, and food oil) and industrial production bases in China. Located near the border of western China, it plays an important role in the development of this region (Liu et al., 2018). Notably, under the joint influence of climate and human activities, precipitation and runoff have decreased in the WRB (Huang et al., 2014, 2017b; Liu et al., 2018). Therefore, it is of important significance to study the evolution characteristics of drought in an integrated manner in the WRB for its sustainable development of social economy.

2.2 Data source The digital data were collected from seven regions of the WRB, referring to the river system and administrative division: Beidao (BD), Linjiacun (LJC), Xianyang (XY), Zhangjiashan (ZJS), Lintong (LT), Huaxian (HX), Zhuangtou (ZT). The monthly precipitation data obtained from 21 meteorological stations in the WRB and its surroundings were utilized as follows (Fig. 1). The monthly runoff data relative to Linjiacun, Zhangjiashan, Huaxian, and other 4 hydrological stations in the WRB were obtained from the hydrologic manual. The Huaxian station is located in the lower reach of the basin, and its catchment area accounts for 97.16% of the entire basin. Hence, the total runoff registered at the Huaxian and Zhuangtou stations can be considered approximately equal to the runoff of the whole basin. In the study area, soil moisture data are very limited for in-situ observations. Due to difficulties during the data collection process, the monthly soil moisture data for depths about 10 cm had to be simulated using a variable infiltration capacity (VIC) model. Sure, it is well known that uncertainties exist in hydrological simulations. Nevertheless, the VIC 8

model performed well in China, with a high resolution of 0.25° × 0.25°, because it employed a high density of meteorological observations from the China Meteorological Administration (CMA) as model inputs and adopted streamflow observations from the Hydrology Bureau of China (HBC) for model calibration (Zhang et al., 2014). The accuracy of soil moisture in this study area simulated by VIC model can be guaranteed (Wu et al., 2011; Wang et al., 2012; Leng et al, 2015; Huang et al., 2017b). The considered meteorological, hydrological, and agricultural data covered the period between January, 1960–December, 2010. Moreover, each time series included 612 specimens. The agricultural drought induced yield losses data adopted in this paper were obtained from the Ministry of Agriculture and Rural Affairs of the People's Republic of China (http://zzys.agri.gov.cn/zaiqing.aspx), which include both the areas covered and affected by droughts (6.67 km2) between years 1960 and 2010.

3 Methodology

The selected precipitation (i.e., meteorological), runoff (i.e., hydrological), and soil moisture (i.e., agricultural) data, together with the empirical frequency cumulative probability distribution (estimated through the Gringorten plotting position formula (1963)), were considered to calculate the combined probability distribution. The developed TSDI could provide an integrated characterization of meteorological, hydrological, and agricultural droughts.

3.1 Trivariate Gringorten empirical frequency formula 9

The joint probability distribution function was employed to combine different drought elements. The precipitation, runoff, and soil moisture data at a specific time scale (e.g., monthly or annual) were represented by the random variables X, Y, and Z, respectively. The joint distribution of these three variables can be expressed as: P( X  x, Y  y, Z  z )  p

(1)

where p denotes the joint probability of the precipitation, runoff, and soil moisture. The TSDI can be defined based on the joint probability p, as follows (Hao and AghaKouchak, 2013):

TSDI   1 ( p)

(2)

where  is the standard normal distribution function. Similar to the MSDI, the TSDI is also derived from the (joint) probability of the selected variables, which can be used to obtain drought information. In Hao and AghaKouchak (2014), the joint distribution expressed in Eq. (1) was also constructed using the Gringorten plotting position formula, but considering only two variables. The alternative method proposed in this paper is based on the nonparametric joint distribution. However, we avoided making any assumptions on the distribution family and tried to alleviate the computational burden in fitting parametric distributions, in order to address a higher-dimension case. The trivariate empirical joint probability can be estimated with the plotting position formula of Gringorten (1963), as follows: H  xi , yi , zi   P  X  xi , Y  yi , Z  zi  

mi  0.44 N  0.12

(3)

where N is the number of the observations and mi is the number of occurrences of the pair (X, Y) for X  xi, Y  yi, and Z  zi ( 1  i  n ). The joint probability 10

derived from Eq. (3) is then applied to Eq. (1) to obtain p. Finally, the obtained p value should be applied to Eq. (2) substitute in order to obtain the TSDI.

3.2 Wavelet analysis and cross wavelet transform The cross wavelet transform (XWT) of the two time series xn and yn is defined as W XY  W X W Y * ,

where * denotes the complex conjugation. The cross wavelet power

can be further defined as W XY . The complex argument arg ( W xy ) can be interpreted as the local relative phase between xn and yn in the time-frequency space (Huang et al., 2016c). The theoretical distribution of the cross wavelet power of two time series with background power spectra PkX and PkY is given in Torrence and Compo (1998) as:  WnX  s WnY *  s   Zv  p  D  p  PkX PkY ,     v X Y  

(4)

Where Zv(p) is the confidence level associated with the probability of a pdf, defined as the square root of the product of two x2 distributions. In this study, the 5% significance level was calculated considering Z2(95%) = 3.999. Further details on the cross wavelet transform can be found in Torrence (1998) and Grinsted (2004).

4 Results

4.1 Comparison between the developed TSDI and other indices One of the targets of this study was to verify the reliability and superiority of the developed TSDI by comparing it with other drought indices, including traditional single-variable and other multivariate integrated drought indices. Meanwhile, the 11

TSDI is compared with a trivariate integrated drought index based on copula function, in order to further prove its convenience and reliability. The results of the comparisons are displayed in the following sections. 4.1.1 Reliability of the developed TSDI Fig. 2 presents some general variations based on different drought indices: the SPI (precipitation-based), SRI (runoff-based), SSI (soil moisture-based), and the TSDI (based on three previously mentioned variables). Generally, these indicators present similar temporal changing patterns. The correlation coefficients between SPI, SRI, SSI, and TSDI are shown in Table 1. The data indicate that drought characterization and monitoring performed using the TSDI can be considered reliable. Notably, the red line representing the TSDI is slightly lower than the other lines, suggesting that the drought conditions captured by the TSDI are a little more serious than those captured by single-element drought index. It can be seen for Table 1that all of them were > 0.65, and most of them were > 0.7, indicating that the developed TSDI is reliable compared with other drought indices. 4.1.2 Superiority of the developed TSDI The first advantage of the developed TSDI compared to other drought indices is its simpler calculation. The copula is a popular method for the formulation of integrated drought indices (e.g., the trivariate integrated drought index by Serinaldi and Grimaldi, 2007). However, the process of fitting the optimized distribution and choosing the appropriate copula function is complex and tedious. Would it be possible to obtain the same result in a simpler way? For comparison purposes, the Frank copula was used to 12

construct a trivariate integrated drought index (MSDI) based on precipitation, rainfall, and soil humidity. Compared with the single-variable drought index, the TSDI has a significant higher correlation with the MSDI (~ 0.95; Table 1). This indicates that the results of the empirical frequency method are not much different from those obtained from the copula function. Besides, Fig. 3 shows how the TSDI and MSDI have approximately the same trend and are highly overlapping. Hao and Aghakouchak (2013) have already demonstrated the advantages of the MSDI in capturing drought events. It can be inferred that the TSDI possesses similar advantages. The second advantage of the TSDI compared to other drought indices is its higher sensitivity to the occurrence of real droughts (Fig. 2). It can be seen from Fig. 2 that the TSDI can integrate the advantages of the SPI and SRI, accurately capturing the onset and the end of each drought. The low index values (with respect to the threshold) occur at different times and with different frequency between the TSDI to SPI, SRI, and SSI. The occurrence of a drought with a specific type does not indicate the occurrence of a drought in an integrated manner. Similarly, when a drought in an integrated manner occurs, a drought with a specific type may not occur. For example, a period characterized by a precipitation deficit can result in a meteorological drought. However, in case of a fast replenishment of the precipitation deficit following runoff surplus, the drought in an integrated manner will not occur, nor have adverse consequences on the human society. It should be noted that a long-term precipitation deficit would lead to 13

decreasing runoff and soil moisture (see the first rectangle in Fig. 4). The SPI did not reach its threshold value, indicating that meteorological drought did not occur. Nevertheless, a drought in an integrated manner did occur under the effects of multiple factors. The second rectangle in Fig. 4 highlights the deficits of precipitation, runoff, and soil moisture, although none of their values reached the threshold of drought. However, a drought in an integrated manner had already occurred. Overall, it can be observed that the drought conditions reflected by the TSDI were closer to the actual conditions of the WRB than those characterized by traditional single-variable indices. Since the TSDI is reliable, easy to calculate, as well as effective, it can be considered a superior tool for effective representation of actual drought conditions in an integrated framework, which is closer to actual drought conditions. Usually, droughts with different types occur asynchronously in time and space. Using single variable drought index cannot characterize drought condition in an integrated manner. As shown in Fig. 5 (a), if only the drought index solely based on runoff or soil moisture is used, droughts caused by precipitation deficits tend to be ignored. Conversely, if the drought index solely based on precipitation, droughts caused by runoff deficits are expected to be ignored (Fig. 5 (b)). The developed TSDI index can capture drought in these two cases. Therefore, the comprehensive drought index integrating multiple variables is of important significance for local drought preparedness and mitigation, since it can be used as an effective reference to guide drought prevention, planning and management. 14

It should be noted that the TSDI is not intended to completely replace the traditional drought index. It exhibits strong correlations with traditional drought index such as SPI, SRI, and SSI, indicating a high reliability of the new developed integrated drought index in this study with multiple drought-related variables. Therefore, the proposed TSDI is convincing an operational drought index, particularly in watershed drought planning and management from an integrated perspective.

4.2 Evolution characteristics of droughts in an integrated manner A more accurate assessment of the drought risk in the WRB can support local planning and sustainable development of water resources as well as drought early warning and mitigation. Therefore, the TSDI was employed to analyse the temporal changing characteristics of drought in an integrated manner in the WRB. 4.2.1 The TSDI trends in the WRB The TSDI trends during 1960–2010 were analysed using the Mann-Kendall method (Yu, 2008). The correspondent results can be seen in Table 2. The TSDI showed a significant downward trend at the 99% confidence level in all the regions of the WRB: during 1960–2010, the drought in an integrated manner in the WRB gradually intensified and its impacts became more serious. In order to improve drought prevention, resistance, and governance, it is necessary to arouse the attention of relevant governmental departments. We further explored the factors influencing the trend changes by calculating the trends in precipitation, streamflow, and soil moisture over the whole basin. The characteristic value of precipitation is Z = -1.67 > -1.96, indicating a nonsignificant 15

decreasing trend at the 95% confidence level in precipitation in the WRB. The Mann-Kendall statistics regarding streamflow and soil moisture has values < -1.96, indicating a significant decreasing trend in streamflow and soil moisture in the WRB at the 95% confidence level. Therefore, the negative trending of precipitation, streamflow and soil moisture jointly contributes to the significant decreasing trend of TSDI in the WRB. 4.2.2 Change points in the annual TSDI series A heuristic segmentation algorithm (Feng, 2005) was utilized to ascertain the occurrence of change points in the annual TSDI, MSDI series of the WRB area. The selected threshold value (P0) and the minimum step (l0) were 0.95 and 25, respectively. The correspondent results are shown in Fig. 6. Change points were detected in both annual TSDI and MSDI series in the year of 1971 and 1993. The results indicate the occurrence of change points in the integrated drought of the WRB and. Hence, the stationarity of the integrated drought series was invalid. 4.3 Climate driving forces of the evolution characteristics of the TSDI Sunspot activity and climate oscillation change the current situation of global water circulation, leading to the spatial and temporal redistribution of water resources and determining the upper limit of available water resources in a region of interest (Sofia and Li, 2001). Solar radiation causes evaporation of water body, and the formation of warm and moist airflow changes the moisture content in the atmosphere (Sofia and Li, 2001). Thus, the distribution characteristics of precipitation and regional hydrological cycle process change (Sofia and Li, 2001; Zhai and Qian, 2017). As for the 16

phenomenon of climate oscillation (Nino3.4, AO, etc.), they alter the general circulation situation and affect regional hydrological elements through the movement of rain belt and the transport of water vapor (Reid, 1991; Fu et al., 2012). The study on climate driving forces of the evolution characteristics of the TSDI is of important scientific significance to the scientific management of regional water resources. Solar radiation can be represented by sunspots, and atmospheric circulation can be represented by atmospheric circulation anomaly factors. Hence, we mainly referred to climate change-related phenomena (i.e., sunspots and abnormal atmospheric circulation factors) as the driving forces behind the evolution of drought conditions in an integrated manner, which would be helpful for drought early warning. 4.3.1 Linkages between the TSDI and sunspots/atmospheric circulation anomaly factors After examining the trends and change points of the TSDI, it is necessary to explore the driving forces behind the observed drought changes in an integrated way. In this study, three driving factors are considered: sunspot activities, sea surface temperatures over the NINO3.4 region, and the Arctic oscillation (AO). Fig. 7 – Fig. 9 show the wavelet correlation maps between sunspots/atmospheric circulation anomaly factors and the TSDI series in the WRB for the period 1960–2010. The arrows indicate the phase differences: those pointing rightwards indicate that the time sequence phase changed accordingly, while those pointing leftwards indicate that the time sequence phase changed in an opposite way. Moreover, the arrows pointing upwards (at 90°) indicate that the driving factor changed 3 months before the TSDI, 17

while those pointing downwards (at 90°) indicate that the driving factor changed 3 months after the TSDI. Finally, higher wavelet transform coefficients correspond to stronger correlations between the factors and the TSDI (the highest correlations are shown in red). Solar radiation affects the quantity of clouds in the Earth’s atmosphere, regulating the cosmic ray flux and thus influences the climate (Sofia and Li, 2001). Fig. 7 shows the cross wavelet transform between sunspots and the TSDI. It can be observed that a negative correlation between the occurrence of sunspots and the TSDI over a period of 7–11 years (during 1967–1985) and a positive correlation over a period of 7–12 years (during 1985–2005) exist. El Niño Southern Oscillation (ENSO) corresponds to the strongest inter-annual change signal in the tropical air-sea coupling system (at low latitudes) (Reid, 1991). The abnormal changes in SST linked to the ENSO have profound effects on the global atmospheric circulation, as well as on weather and climate variability in China (especially over the NINO3.4 region) (Reid, 1991). Fig. 8 shows the cross wavelet transform between the NINO3.4 index and the TSDI. It can be seen that negative correlations between the ENSO and the TSDI exist over a period of 1–4 years (during 19621–1974) and 4–6 years (during 19841–1993). The AO represents the main mode of low-rate variation in the extratropical atmosphere of the northern hemisphere (Toreti et al., 2010). Climate anomalies in northern China are closely related to the AO (Wang et al., 2013). Fig. 9 shows the cross wavelet transforms between the AO and TSDI. It can be noticed a positive 18

correlation between the AO and the TSDI over a first period of 1–4 years (during 1960–1970) a second period of 7–8 years (during 1980–1987) and a third period of 4–5 years (during 1988–1992) exist. Overall, the variation occurred within the range of the rough black line outline. Meanwhile, the phase of the cross wavelet changed before and after the year of the change point of TSDI. These observations further prove the reliability of the captured change points. Moreover, they illustrate the relationship between the occurrence of variations in the TSDI sequence and that of environmental factors (i.e., sunspots). Some differences could be observed between the cross wavelet graph of sunspot occurrence and those of the NINO3.4 and AO indices. Sunspots have a stronger influence on the TSDI sequence, both in term of time and degree. Hence, sunspot activity was the dominant driving force on the drought in an integrated manner in the WRB to certain extent. Fig.7 and Fig. 9 suggest that the effects of different driving factors on the TSDI were different in the time and frequency domains. For the same driving factor, its effect is within different time and frequency domains, and the phase and retardation time of the driving factor with the TSDI are different, even in cases when their features are similar. As shown in Fig. 7, between the 1970s and the year of 2005, the sunspot and the TSDI phases changed from negative to positive. Thus, the influence mode of the sunspot on the TSDI possibly have changed. 4.3.2 Comparison between the TSDI and other single-variable indices related with sunspots 19

Fig.7 and Fig. 10 illustrate the cross wavelet maps between sunspots and the TSDI as well as the single-variable drought indices, respectively. Comparing these figures, we could notice how the TSDI effectively integrated precipitation, runoff, and soil moisture information, and it is more integrated than the univariate drought index. In Fig. 10 (b) and (c) the time domain 1987–1988 corresponded to a frequency domain of approximately 5 years. Both of them were included in the 95% confidence interval but could not be observed in Fig. 10 (a). Moreover, in Fig. 10 (b) the time domain 1998–2005 corresponded to a frequency domain of 8–12 years, while in Fig. 10 (c) the time domain 1974–1980 corresponded to a frequency domain of 6–7 years. These peculiar circumstances appeared all combined and offset in Fig. 7. The TSDI, being based on different drought-related variables, is expected to characterize drought conditions more accurately than single-variable drought indices. 4.3.3 Linkages between sunspots and atmospheric circulation anomaly factors Sunspots tend to have a very strong influence on both the AO and NINO3.4 indices (Fig. 11). Fig. 11 (a) shows the existence of a negative correlation between sunspots and the AO over a period of 8–12 years (during 1967–1980) and a positive correlation over a period of 7–14 years (during 1980–2007). Fig. 11 (b) shows the existence of a positive correlation between sunspots and the AO index over a period of 10–14 years (during 1967–1993). In addition, the arrows in the black coil in Fig. 11 are mostly downward, indicating that the AO and Nino3.4 indices present a certain degree lagged with sunspot changes. In other words, sunspots appeared to have a strong influence on these indices, thus indirectly influencing the TSDI. 20

Overall, sunspot activities had not only direct but also indirect impacts on the drought in an integrated manner in the WRB, mainly by influencing atmospheric circulation anomaly factors. Besides, sunspots seemed to have the strongest effect on the drought in an integrated way. Therefore, among all teleconnection factors, sunspot activities likely were the main driving force on the drought evolution in an integrated manner in the WRB.

5 Discussion

5.1 Comparison between parametric and nonparametric indices associated with yield losses Due to the lack of systematic and multifactor integrated drought data, integrated drought indices like the TSDI and MSDI cannot be directly compared with the measured data. In order to quantify the difference between the parametric and nonparametric drought indices, the SSI (closely related to agricultural drought) was compared with the drought induced yield losses data. In this study, the parametric SSI was referred to as SSIp, and the empirical SSI was referred to as SSIe. The SSIe was based on the Gringorten formula, while the SSIp was based on a generalized extreme distribution following a series of steps (e.g., selecting different distributions, estimating parameters, and performing a fitting optimization test) and standardized the result. In order to develop the SSIp, four distributions (i.e., Weibull, Gaussian, gamma, and generalized extreme distributions) were selected, and the root mean squared error (RMSE) was used to test the goodness 21

of fit (Table 4). The gamma distribution (recommended by the Chinese National Weather Service for hydrological frequency analysis) performed well, but not as much as the generalized extreme distribution. Table 5 shows the correlation between the SSI and the drought induced yield losses data (e.g., the SSIp based on the gamma distribution, the SSIp based on the generalized extreme distribution, and the SSIe). Table 5 shows the linkage between the SSIe and yield losses data is stronger than the SSIp based on the gamma distribution with good fitting effect. The results partly prove that a drought index based on an empirical frequency has a rather good performance. Hence, it is possible to avoid numerous tedious calculation steps, while maintaining a high level of reliability. 5.2 Analysis of the drought elements’ change points in the whole WRB To further analyse the causes of the observed change points in the TSDI series in the WRB, we performed a change point identification for the three drought-related variables.

Fig. 12 shows the results of the change point detection for precipitation,

streamflow, and soil moisture in the WRB with the period during 1960–2010. No change points were identified in the precipitation sequence, while some change points were detected in the annual streamflow series in corresponding year of 1969 and 1990, and in the soil moisture series in corresponding year of 1985. In the context of global warming, no significant declining trends in precipitation is observed. Variations in streamflow occurred in the early 1970s, when large-scale water conservancy projects and irrigation ditches were constructed. During the 1990s, the water consumption of the national economy greatly increased. An abrupt change 22

in soil moisture occurred in 1985, which was perhaps related to the application of cultivated land policies and measures for the protection of forest lands (Huang et al., 2015). The three drought-related variables changed differently. Therefore, the extrapolation of aberrances in drought from either of these conditions alone would not provide reliable results. In this case, the change points identified in the TSDI trend should provide a more realistic representation of the characteristics of drought variation than those identified in the MSDI series. In the 1990s, the water consumption of residents, national economy, soil and water conservation, and non-use caused by evaporation (linked to warming) all increased (Huang et al., 2015). These above factors could have all contributed to drought variability. The hydrological sequence in the WRB became non-stationary. Therefore, it can be reasonably concluded that the variations in runoff and soil moisture were caused by both climate change and human activities (Huang et al., 2015). In the case of the TSDI combining three factors, the impact of human activities might account for a larger proportion. The described quantitative methods can be used in the future study to further analyse the contributions to the variations in the TSDI series. 5.3 The period of the TSDI The period of the TSDI has been obtained based on wavelet analysis (Rao, 1996). The time-frequency distribution of the annual TSDI series is shown in Fig. 13. In Fig. 13, different colours correspond to different phases. Red and blue indicate positive and negative phases, respectively. The interval between the same phases represents the period of the TSDI series. In the wavelet analysis diagram, the horizontal 23

coordinate represents the period. A primary period of up to 25 years and a secondary period of nearly 7 years over the whole WRB during 1960–2010 are identified. Similar results have been obtained by Huang et al. (2016a). A careful observation of the 7-years period in Fig. 13 (a) shows that the phase changes significantly from one to the other side of the red vertical dotted line. The two dotted lines corresponding to the year of 1971 and 1993, which are the two change points in the TSDI series (Fig. 6). These observations confirm our previous analysis on the variation in TSDI. Additionally, the periods of the TSDI tended to increase after the year of 1990, independently from the primary and secondary periods.

6 Conclusions

A trivariate standardized integrated drought index based on the empirical frequency cumulative probability distribution was constructed in this study, which combines the information of precipitation, streamflow, and soil moisture data. This integrated drought index is developed on a multivariate joint distribution function whose marginal distributions are different, thus avoiding subjective errors caused by subjective empowerment and the non-linear distortions linked to the PCA. The use of copulas is avoided in the development of the drought index. The calculation of TSDI is simple and reliable at the same time. Overall, it is demonstrated that the developed TSDI is able to characterize meteorological, hydrological and agricultural droughts in an integrated manner. Although the TSDI was just applied to the WRB in this study, its use can be extended to any other regions without limitation. Our main conclusions 24

are as follows: (1) The drought conditions reflected by the TSDI are closer to reality than those reflected by traditional single-variable indices. The results of the TSDI were similar to those of a trivariate drought index based on copula function, which but was simpler to calculate and better applicable to different regions. (2) A significant downward trend with strong persistence was observed for the TSDI in all regions of the WRB, which indicated a gradual intensification of the drought in an integrated manner in the WRB during 1960–2010. (3) Two change points in the TSDI series were identified in the year of 1971 and 1993, indicating that the drought in an integrated manner in the WRB is non-stationary. This characteristic is probably linked to a combination of changing climate (section 4.3.1) and human activities (section 5.2), especially to human activities. (4) A primary period of up to 25 years and a secondary period of approximately 7 years for the TSDI series are identified over the whole WRB during 1960–2010. Moreover, its period changed before and after the change points of the TSDI. (5) The evolution of the drought in an integrated manner in the WRB was closely related to sunspot activities and atmospheric circulation anomaly factors. Among them, sunspot activities exert the dominant influence. Sunspot activities show both direct and indirect impacts on the drought in an integrated framework. Such conclusions were confirmed by the results of the change point identification, in which the relationships between sunspots, atmospheric circulation anomaly factors, and the TSDI before and after the year of variation changed considerably. 25

Overall, this study proved the simplicity and reliability of a new trivariate integrated drought index, named TSDI, which provided a valuable reference for drought index development. Acknowledgements This research was jointly funded by the National Natural Science Foundation of China (grant number 51709221), the National Key Research and Development Program of China (grant number 2017YFC0405900), the Planning Project of Science and Technology of Water Resources of Shaanxi (grant numbers 2015slkj-27 and 2017slkj-19), China Postdoctoral Science Foundation Grant (2018M640155), the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China Institute of Water Resources and Hydropower Research, grant number IWHR-SKL-KF201803), the Belt and Road Special Foundation of the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (2018490711), and the Doctorate Innovation Funding of Xi'an University of Technology (grant number 310-252071712). Authors would like to extend our sincere appreciation to the editor and two anonymous reviewers for their constructive comments, which help to improve the quality of the manuscript substantially.

References Banerjee, S., Mitra, M., 2013. Application of cross wavelet transform and Wavelet Coherence for classification of ECG patterns. International Conference on Emerging Applications of Information Technology (pp.89-92). IEEE. Fu, C., James, A. L., Wachowiak, M.P., 2012. Analyzing the combined influence of solar activity and el niño on streamflows across southern canada. Water Resources Research, 48(5). 26

Feng, G. L., Gong, Z. Q., Dong, W. J., et al., 2005. Abrupt climate change detection based on heuristic segmentation algorithm. Acta Physica Sinica, 54(11), 5494-5499. Fang W., Huang S., Huang G., et al., 2019a. Copulas-based risk analysis for inter-seasonal combinations of wet and dry conditions under a changing climate. International Journal of Climatology, 39(4), 2005-2021. Fang, W., Huang, S.Z., Huang, Q., et al., 2019b. Probabilistic assessment of remote sensing-based terrestrial vegetation vulnerability to drought stress of the Loess Plateau in China. Remote Sensing of Environment, 232, 111290. Grinsted, A., Moore, J. C., Jevrejeva, S., 2004. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics, 11(5/6), 561-566. Gringorten, I. I., 1963. A plotting rule for extreme probability paper.J. Geophys. Res., 68, 813–814, doi:10.1029/JZ068i003p00813. Guo, Y., Huang, S., Huang Q., et al., 2019a. Assessing socioeconomic drought based on an improved Multivariate Standardized Reliability and Resilience Index. Journal of Hydrology, 568, 904-918. Guo, Y., Huang, S., Huang, Q., et al., 2019b. Copulas-based bivariate socioeconomic drought dynamic risk assessment in a changing environment. Journal of Hydrology, 575, 1052-1064. Hao, Z., Aghakouchak, A., 2013. Multivariate standardized drought index: a parametric multi-index model. Advances in Water Resources, 57(9), 12-18. Hao, Z., Aghakouchak, A., 2014. A nonparametric multivariate multi-index drought monitoring framework. Journal of Hydrometeorology, 15(1), 89-101. Huang S., Chang J., Huang Q., et al., 2014. Spatio-temporal changes and frequency analysis of drought in the Wei River Basin, China. Water Resources Management 28(10), 3095-3110. Huang, S., Chang, J., Huang, Q., et al., 2015. Identification of abrupt changes of the relationship between rainfall and runoff in the Wei River basin, China. Theoretical and Applied Climatology 120(1-2), 299-310. Huang, S., Huang, Q., Leng, G., et al., 2016a. A hybrid index for characterizing drought based on a nonparametric kernel estimator. Journal of Applied Meteorology & Climatology, 55(6). Huang, S., Huang, Q., Leng, G., et al., 2016b. A nonparametric multivariate standardized drought index for characterizing socioeconomic drought: A case study in the Heihe River Basin. 27

Journal of Hydrology, 542 875-883. Huang, S., Huang, Q., Chang, J., et al., 2016c. Linkages between hydrological drought, climate indices and human activities: a case study in the Columbia River basin. International Journal of climatology 36(1), 280-290. Huang, S., Li, P., Huang, Q., et al., 2017a. The propagation from meteorological to hydrological drought and its potential influence factors. Journal of Hydrology, 547, 184–195 Huang, S., Huang, Q., Leng, G., et al., 2017b. Variations in annual water-energy balance and their correlations with vegetation and soil moisture dynamics: a case study in the wei river basin, china. Journal of Hydrology, 546, 515-525. Huang, S., Wang, L., Wang, H., et al., 2019. Spatio-temporal characteristics of drought structure across China using an integrated drought index. Agricultural Water Management 218, 182-192. Han, Z., Huang, S., Huang, Q., , et al., 2019. Assessing GRACE-based terrestrial water storage anomalies dynamics at multi-timescales and their correlations with teleconnection factors in Yunnan Province, China. Journal of Hydrology, 574, 836-850. Keyantash, J. A., Dracup, J. A., 2004. An aggregate drought index: assessing drought severity based on fluctuations in the hydrologic cycle and surface water storage. Water Resources Research, 40(9), 333-341. Kao, S. C. , Govindaraju, R. S.,2010. A copula-based joint deficit index for droughts. Journal of Hydrology (Amsterdam), 380(1-2), 121-134. IPCC (Intergovernmental Panel on Climate Change), 2013, Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment. Leng, G., Tang, Q., Rayburg, S., 2015. Climate change impacts on meteorological, agricultural and hydrological droughts in China. Global and Planetary Change, 126, 23-34. Liu, S., Huang, S., Xie, Y., et al., 2018. Spatial-temporal changes of rainfall erosivity in the loess plateau, China: Changing patterns, causes and implications. Catena, 166, 279-289 Liu, S., Huang, S., Xie, Y., et al., 2019. Spatial-temporal changes in vegetation cover in a typical semi-humid and semi-arid region in China: Changing patterns, causes and implications. Ecological indicator, 98, 462-475. Liu, Y., Zhu, Y., Ren, L., et al., 2017. A multiscalar palmer drought severity index. Geophysical Research Letters, 44(13), 6850-6858. 28

McKee, T. B., N. J. Doesken, J. Kleist, 1993. The relationship of drought frequency and duration to time scales. Preprints, Eighth Conf. on Applied Climatology, Anaheim, CA, Amer. Meteor. Soc., 179–184. Mishra, A. K., Singh, V. P., 2010. A review of drought concepts. Journal of Hydrology, 391(1), 202-216. Milly, P. C. D., Wetherald, R. T., Dunne, K. A. , et al., 2002. Increasing risk of great floods in a changing climate. Nature, 415, 514–517. Mussá, F. E. F., Zhou, Y., Maskey, S., et al., 2015. Groundwater as an emergency source for drought mitigation in the Crocodile River catchment, South Africa. Hydrology and Earth System Sciences, 19, 1093–1106. Palmer, W., 1965. Meteorological drought. Weather Bureau Research Paper 45. U.S. Weather Bureau, Washington, DC, 58 pp. Reid, G. C., 1991. Solar total irradiance variations and the global sea surface temperature record. Journal of Geophysical Research: Atmospheres, 96. Rao, T. S., Indukumar, K. C., 1996. Spectral and wavelet methods for the analysis of nonlinear and nonstationary time series. Journal of the Franklin Institute, 333(3), 425-452. Serinaldi, F., Grimaldi, S., 2007. Fully nested 3copula: procedure and application on hydrological data. Journal of Hydrologic Engineering, 12(4), 420-430. Shukla, S., A. W. Wood, 2008. Use of a standardized runoff index for characterizing hydrologic drought. Geophys. Res. Lett., 35, L02405, doi:10.1029/2007GL032487. Serinaldi, F., Grimaldi, S., 2007. Fully nested 3-copula: procedure and application on hydrological data. Journal of Hydrologic Engineering, 12(4), 420-430. Sofia, S., Li, L., 2001. Solar variability and climate. Journal of Geophysical Research, 106(A7), 12969. Svoboda, M., Lecomte, D., Hayes, M., et al., 2002. The drought monitor. Bulletin of the American Meteorological Society, 83(8), 1181-1190. Sheffield, J., Wood, E. F., 2008. Global trends and variability in soil moisture and drought characteristics. Journal of Climate, 21, 432–458. Toreti A, Xoplaki E, Maraun D, et al., 2010. Characterisation of extreme winter precipitation in Mediterranean coastal sites and associated anomalous atmospheric circulation patterns. Natural Hazards and Earth System Science 10, 1037-1050. 29

Torrence C, Compo G P.,1998. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society. 79(1), 61-78. Vicente-Serrano, S. M. , Beguería, Santiago, et al., 2010. A multiscalar drought index sensitive to global warming: the standardized precipitation evapotranspiration index. Journal of Climate, 23(7), 1696-1718. Wang, A., Lettenmaier, D. P., Sheffield, J., 2011. Soil moisture drought in China, 1950–2006. Journal of Climate, 24(13), 3257-3271. Wang WG, Shao QX, Yang T, et al., 2013. Changes in daily temperature and precipitation extremes in the Yellow River Basin, China. Stochastic Environmental Research and Risk Assessment 27(2), 401-421. Wu, J., Chen, X., Yao, H., et al., 2018b. Hydrological drought instantaneous propagation speed based on the variable motion relationship of speed‐ time process. Water Resources Research, 54(11), 9549–9565. Wu J., Liu Z., Yao H., et al., 2018a. Impacts of reservoir operations on multi-scale correlations between hydrological drought and meteorological drought. Journal of Hydrology, S0022169418304761-. Wu, Z., Lu, G., Wen, L., et al., 2011. Reconstructing and analyzing China's fifty-nine year (1951–2009) drought history using hydrological model simulation. Hydrology and Earth System Sciences, 15(9), 2881-2894. Yang, J., Chang, J., Wang, Y., et al., 2018. Comprehensive drought characteristics analysis based on a nonlinear multivariate drought index. Journal of Hydrology, 557, 651-667. Yan-Sheng, Y. U. , Xing-Wei, C., 2008. Analysis of future trend characteristics of hydrological time series based on r/s and mann-kendall methods. Journal of Water Resources & Water Engineering. Zhai, Qian., 2017. Influence of solar activity on the precipitation in the north-central china. New Astronomy, 51, 161-168. Zhang, X., Tang, Q., Pan, M., 2014. A long-term land surface hydrologic fluxes and states dataset for China. J Hydrometeor. 15, 957-972.

30

Fig. 1. Location of the Wei River Basin.

Fig. 2. Comparison between monthly TSDI and three single variable drought indexes from 1960 to 2010 in the WRB.

31

Fig. 3. Comparison for monthly TSDI and MSDI from 1960 to 2010 in the WRB.

Fig. 4. Details of comparison for monthly TSDI and 3 single variable drought indexes from 1960 to 2010 in the WRB.

32

(a)

(b) Fig. 5 Schematic diagram of drought in different situations

33

Fig. 6. The result of change point identification of the TSDI & MSDI series in the WRB.

Fig. 7. The wavelet analysis of annual TSDI series and sunspot covering 1960-2010 in the WRB. 34

Fig. 8. The wavelet analysis of annual TSDI series and Nino 3.4 covering 1960-2010 in the WRB.

Fig. 9. The wavelet analysis of annual TSDI series and AO covering 1960-2010 in the WRB. 35

（b）

（a）

（c）

Fig. 10. The wavelet analysis of SPI, SRI, SSI series and sunspot covering 1960-2010 in the WRB.

36

(a)

(b)

Fig. 11. The wavelet analysis of AO, Nino3.4 and sunspot covering 1960-2010 in the WRB.

37

Fig. 12. The result of change point identification of precipitation, streamflow & soil moisture series in the WRB.

38

(a)

(b)

Fig. 13. The time-frequency distribution of annual TSDI series covering 1960-2010 in the WRB (a); the wavelet analysis of annual TSDI series covering 1960-2010 in the WRB (b).

39

Table 1 The Pearson correlation coefficient r between the TSDI and different drought indexes in all partitions. Index

BD

LJC

XY

ZJS

LT

HX

ZT

The whole basin

SPI SRI SSI MSDI

0.72 0.71 0.76 0.95

0.71 0.73 0.77 0.95

0.76 0.78 0.80 0.96

0.69 0.75 0.77 0.95

0.75 0.77 0.78 0.96

0.75 0.77 0.78 0.94

0.66 0.73 0.77 0. 95

0.79 0.78 0.81 0.96

Note: r represents the Pearson correlation coefficient result, if 0.5  r  0.8 , indicating a significant correlation between two series; if 0.8  r  1 , indicating a high correlation between two series.

Table 2 The Z values of MK statistics of TSDI in each region of WRB during 1960-2010.

Z value

BD

LJC

XY

ZJS

LT

HX

ZT

The whole basin

-4.4**

-4.21**

-2.96**

-4.08**

-2.97**

-3.29**

-4.54**

-3.04**

Note: Z value represents the Mann-Kendall test result, if Z  1.64 , indicating a failure to reject hypothesis at 95% confident level; if Z  2.32 , indicating a failure to reject hypothesis at 99% confident level. “*” and “**” represent significant at the 95% and 99% confidence level, respectively.

40

Table 3 The statistics of the MK of precipitation, streamflow and soil moisture in the WRB between 1960-2010.

Z value

precipitation

streamflow

soil moisture

-1.67*

-3.8**

-3.10**

Note: “*” and “**” represent significant at 95% and 99% confidence level, respectively.

Table 4 The Root mean squared error (RMSE) of 4 distributions.

Distributions

Weibull distribution

Gaussian distribution

Gamma distribution

Generalized extreme distribution

RMSE

0.102

0.259

0.058

0.012

Table 5 The Correlation between SSI based on different distributions and actual agricultural drought damage data.

Drought attacked areas Drought inundated areas

Gamma distribution

Generalized extreme distribution

Empirical frequency distribution

-0.596

-0.614

-0.618

-0.596

-0.632

-0.626

41

Highlights  An alternative nonparametric trivariate integrated drought index was proposed.  It is easier to calculate but more reliable in capturing droughts compared with univariate and parametric ones.  The drought in an integrated manner gradually intensified in the WRB.  The non-stationarity of integrated drought series is mainly caused by human activities and sunspots.  Sunspots are the dominant teleconnection factors exerting both direct and indirect effects on drought evolution.

42