Characterizing Hydrological Drought and Water Scarcity Changes in the Future: A Case Study in the Jinghe River Basin of China

The assessment of future climate changes on drought and water scarcity is extremely important for water resources management. A modeling system is developed to study the potential status of hydrological drought and water scarcity in the future, and this modeling system is applied to the Jinghe River Basin (JRB) of China. Driven by high-resolution climate projections from the Regional Climate Modeling System (RegCM), the Variable Infiltration Capacity model is employed to produce future streamflow projections (2020–2099) under two Representative Concentration Pathway (RCP) scenarios. The copula-based method is applied to identify the correlation between drought variables (i.e., duration and severity), and to further quantify their joint risks. Based on a variety of hypothetical water use scenarios in the future, the water scarcity conditions including extreme cases are estimated through the Water Exploitation Index Plus (WEI+) indicator. The results indicate that the joint risks of drought variables at different return periods would decrease. In detail, the severity of future drought events would become less serious under different RCP scenarios when compared with that in the historical period. However, considering the increase in water consumption in the future, the water scarcity in JRB may not be alleviated in the future, and thus drought assessment alone may underestimate the severity of future water shortage. The results obtained from the modeling system can help policy makers to develop reasonable future water-saving planning schemes, as well as drought mitigation measures.


Introduction
Drought and water scarcity (WS) are inevitable research topics for water resource management in water-stressed areas [1][2][3]. Drought is a natural disaster that affects a large number of people and causes huge economic losses in the world [4,5]. WS is usually caused by the combined effects of natural factors and excessive water uses of human beings, and has already become a vital obstacle for socio-economic development in many parts of the world [6]. In the context of global climate change, many traditional arid regions may experience more severe drought events in the future, which have even been observed in the past few years [7][8][9]. Moreover, due to the increase in population and the improvement of living standards, water usage is significantly growing, which will result in massive WS in numerous regions in the future [2,10,11]. Consequently, it is desired to explore the severity of drought and further reveal the associated WS under climate change in order to manage water resources reasonably and meet the needs of sustainable development [6,11].
Drought can be classified from meteorology, agriculture, hydrology and socioeconomic aspects [12]. Among the four types of drought, hydrological drought (HD) can characterize the reduction of water

Runoff Simulation and Drought Identification
In this study, the climate projections (2020-2099) are employed to drive the Variable Infiltration Capacity (VIC) model [40,41] to obtain future runoff for two RCP scenarios. The VIC hydrological model is run at a 0.22-degree resolution (24-h time scale), and then the runoff and base flow of the simulated grid cells are processed by a river routing model [42]. Detailed information for the VIC

Runoff Simulation and Drought Identification
In this study, the climate projections (2020-2099) are employed to drive the Variable Infiltration Capacity (VIC) model [40,41] to obtain future runoff for two RCP scenarios. The VIC hydrological model is run at a 0.22-degree resolution (24-h time scale), and then the runoff and base flow of the simulated grid cells are processed by a river routing model [42]. Detailed information for the VIC model can be found in Zhang et al. [43]. The VIC model is calibrated for the period of 1981-1990 and validated for the period of 1998-2005 before being used to generate future runoff. The performance of the VIC model is assessed by the Nash-Sutcliffe efficiency coefficient (NSE, Equation (1)) [44,45] and determination coefficient (R 2 , Equation (2)).
It is worth noting that the simulation results of annual runoff volume are needed in the subsequent calculation process of WEI+. Although NSE and R2 can well reflect the fluctuation of monthly runoff time series, they cannot reflect the deviation degree of runoff volumes. Therefore, the deviation of the runoff volumes (DV, [46]) indicator is then calculated. The closer the DV value is to 1, the smaller the volumetric deviation is: where V obs,t and V sim,t represent the observed and simulated annual runoff volumes in year t, respectively; n is the number of the data.
The SDI is used to identify hydrological drought events and their corresponding characteristics (duration and severity). SDI can be calculated from monthly streamflows: where m is the time scale (can be 1, 3, 6, 12 et al.); Q i,j is the runoff in month j, year i; V i,j denotes the cumulative streamflow of the chosen time scale. The SDI values can be obtained by normalizing the V i,j values [13]: where V and s V represent the average and standard deviation of V i,j , respectively. The thresholds of HD according to SDI values are shown in Table 1 [13]. In this study, duration and severity of HD events are derived from monthly SDI values using the threshold approach [47][48][49]. The duration of a HD event is the time period from the beginning to the end of the event, and severity is the cumulative deviation of the SDI below the threshold. More details of HD identification process and HD duration and severity meanings are presented in Figure S1 and its explanations in the Supplementary Materials. In reality, duration and severity of HD events are correlated with each other. An individual variable (duration or severity) of the HD event can hardly provide the event a comprehensive description and thus may lead to misestimation of the HD risk. Therefore, the correlation between HD variables should be carefully taken into account in HD risk analysis. Developed by Sklar [50], copulas can be used to construct multivariate distributions in a flexible way in choosing marginal distributions and their dependence structures [51,52]. A two-dimensional copula function can be expressed as: where (F X , F Y ) represent the cumulative distribution functions (CDFs) of a two-dimensional random vector (X, Y).
Marginal distributions are fitted for the variables before establishing the bivariate distribution in Equation (6). Gamma, lognormal (LN), Weibull, generalized extreme value (GEV), Pearson type III (P-III), and Log Pearson type III (LP III) distributions are employed to establish marginal distributions for HD variables. The goodness-of-fit tests for marginal distributions are implemented by using the Akaike's information criteria (AIC) [53], root mean square error (RMSE) and Kolmogorov-Smirnov (KS) statistic test [54].
Because of the simplicity and good accuracy properties, Archimedean copulas have attracted wide attention in drought risk analyses [15,55]. Three Archimedean copulas: Frank, Clayton, and Gumbel-Hougaard are chosen as candidate models to derive bivariate distribution for the duration and severity of HD events. The parameters in the chosen copulas are obtained through the method-of-moments-like (MOM) estimator [56]. The Cramér-von Mises statistic (CvM) test [57], AIC and RMSE are used for goodness-of-fit tests of copulas.

Drought Risk Analysis
In this study, different kinds of return periods are adopted to quantify the risk of drought events. The return period is the mean interval time between two consecutive drought events. The univariate return periods of the two HD variables can be derived by estimating the exceedance probability from the fitted marginal distributions: where T D /T S denotes the return period of duration (D)/severity (S) larger than d/s; E(L) represents the average interarrival time of the HD variable; F D (d) and F S (s) are CDFs of duration and severity, respectively. Two types of joint return periods (JRPs) are considered in this study to quantitatively evaluate the joint risk of the two variables: where T AND DS denotes the JRP when D and S respectively exceed specific values of d and s; T OR DS denotes the JRP when D or S is higher than the certain thresholds (i.e., D ≥ d or S ≥ s). F DS (d,s) is the joint CDF of D and S.

Water Exploitation Index +
The WS indicator used in this study is the WEI+. It is applied to analyze the level of pressure exerted by human activities on natural water resources. WEI+ is calculated as the ratio of water withdrawal minus return water (i.e., net water consumption, in m 3 ) and renewable water resources (in m 3 ) [26,27]: where Abstractions refers to the total water resources extracted from the watershed (containing surface and ground water); Returns indicates the amount of water resources that have been exploited and re-entered into the hydrological cycle; RWR is the amount of renewable water resources. For basins influenced by human activities, two methods exist for renewable water calculation [26]: Method 2: RWR = Out f low + (Abstractions − Returns) − ∆S art (13) where ExIn (External Inflow) indicates the amount of surface water and groundwater inflow from outside the basin; P refers to the total precipitation; Eta (Evapotranspiration) denotes the total evapotranspiration; ∆S nat refers to the change of water resources reserves in natural environment; Outflow refers to the actual water resources flowing out of the drainage basin; ∆S art denotes changes in water reserves in artificially regulated lakes or artificial reservoirs. Thresholds of water stress according to the WEI+ are shown in Table 1 [27,58].

Hydrological Model Verification and Hydrological Drought Identification
The NSE and R 2 calculated from monthly flow discharges are 0.80 and 0.86, respectively, for the calibration period 1981-1990, and are 0.78 and 0.82, respectively, for the validation period of 1998-2005. These indicate that the runoff simulations and observations are in good agreement. Figure 2 shows the comparison between the monthly runoff simulated from climate model projections and from the gridded observations of 1987-2004 (with NSE and R 2 are 0.71 and 0.78, respectively). The simulated and observed annual runoff volumes produce an DV value of 1.04, which is very close to 1. All these indicate that the simulated runoff series driven by the outputs of RegCM perform well in capturing the runoff variations in observations, and thus the future runoff simulations can then be used for HD-and WS-related calculations. The monthly streamflow series of observations from 1960 to 2009 and projections from 2020 to 2099 are used to generate the SDI values for both history and two future scenarios. The SDI variations with corresponding monthly precipitation series under history and projected future scenarios are depicted in Figure 3, which helps understand the runoff deficit under different precipitation conditions. The average annual precipitation under RCP4.5 and RCP8.5 scenarios would increase by 108.1 mm and 135.9 mm respectively compared with that in the historical period. The SDIs over the historical period changes from −3.46 to 3.44, and SDIs under scenarios of RCP4.5 and RCP8.5 have ranges of [−1.75, 5.45] and [1.46, 4.06], respectively. The mean values of SDIs for RCP4.5 and RCP8.5 are, respectively, 0.29 and 0.33 higher than those in the historical period. We can preliminarily conclude that hydrological droughts would be relieved under projected future climate change scenarios.
HD events with their corresponding deficit characteristics (duration and severity) are identified from monthly SDI series. Table 2 provides the statistical information about the characteristics of HD occurrence, persistence and severity under history and two projected future scenarios. Mean interval time of HD events would be increased by 0.96 and 2.08 months in RCP4.5 and RCP8.5, respectively. Average HD duration and severity for projected future scenarios would slightly decrease when compared with those in the historical period. The maximum value of severity for historical HD events is 7.97 and it would be respectively reduced to 4.98 for RCP4.5 and 3.67 for RCP8.5. The Kendall correlation coefficient (τ) [19] between drought duration and severity under historical and future scenarios ranges from 0.60 to 0.67, which indicates that there is a strong correlation between the two variables. The monthly streamflow series of observations from 1960 to 2009 and projections from 2020 to 2099 are used to generate the SDI values for both history and two future scenarios. The SDI variations with corresponding monthly precipitation series under history and projected future scenarios are depicted in Figure 3, which helps understand the runoff deficit under different precipitation conditions. The average annual precipitation under RCP4.5 and RCP8.5 scenarios would increase by 108.1 mm and 135.9 mm respectively compared with that in the historical period. The SDIs over the historical period changes from −3.46 to 3.44, and SDIs under scenarios of RCP4.5 and RCP8.5 have ranges of [−1.75, 5.45] and [1.46, 4.06], respectively. The mean values of SDIs for RCP4.5 and RCP8.5 are, respectively, 0.29 and 0.33 higher than those in the historical period. We can preliminarily conclude that hydrological droughts would be relieved under projected future climate change scenarios.
HD events with their corresponding deficit characteristics (duration and severity) are identified from monthly SDI series. Table 2 provides the statistical information about the characteristics of HD occurrence, persistence and severity under history and two projected future scenarios. Mean interval time of HD events would be increased by 0.96 and 2.08 months in RCP4.5 and RCP8.5, respectively. Average HD duration and severity for projected future scenarios would slightly decrease when compared with those in the historical period. The maximum value of severity for historical HD events is 7.97 and it would be respectively reduced to 4.98 for RCP4.5 and 3.67 for RCP8.5. The Kendall correlation coefficient (τ) [19] between drought duration and severity under historical and future scenarios ranges from 0.60 to 0.67, which indicates that there is a strong correlation between the two variables.

Univariate and Bivariate Distributions
The univariate distributions for HD variables are achieved based on the identified HD records in both historical and future periods. Table 3 shows the parameters of AIC, RMSE and KS test results for the best fitted distributions for HD duration and severity. It can be seen that the Gamma and LN distributions would be employed for HD duration and severity in a historical period. The GEV distribution is suitable for both HD duration and severity under the RCP4.5 scenario. The HD duration and severity under the RCP8.5 scenario are fitted with GEV and Weibull distributions, respectively. The p-values of the KS test results for all fitted distributions are greater than 0.05, which indicate that these distributions can effectively describe the probability characteristics of the HD variables. Figure 4 illustrate the comparison of empirical and theoretical CDFs for HD duration and severity under historical period and projected future scenarios. It shows that the theoretical CDFs from best fitted distributions of the HD variables are very close to empirical distributions.  The dependence structure between HD duration and severity is then characterized by the copula functions. Table 4 shows the statistical test results for best fitted copulas. The Clayton, Gumbel and Frank copulas are, respectively, the best choices in modeling the joint distributions of HD duration and severity under history and two future scenarios. Figure 5 displays the bivariate CDFs based on the fitted copulas for HD duration and severity under history and two future scenarios. Table 4. Goodness-of-fit tests of best fitted copulas for HD duration and severity under historical The dependence structure between HD duration and severity is then characterized by the copula functions. Table 4 shows the statistical test results for best fitted copulas. The Clayton, Gumbel and Frank copulas are, respectively, the best choices in modeling the joint distributions of HD duration and severity under history and two future scenarios. Figure 5 displays the bivariate CDFs based on the fitted copulas for HD duration and severity under history and two future scenarios. Table 4. Goodness-of-fit tests of best fitted copulas for HD duration and severity under historical periods and projected future scenarios.

Return Periods of Hydrological Drought Events
Six return period levels from 3 to 100 (in Figure 6) for HD duration and severity are estimated according to Equations (7) and (8). Figure 6 shows that the values of the two variables for all return periods in both future scenarios would decrease significantly, especially in the RCP8.5 scenario. Under the RCP4.5 scenario, the 3-year and 100-year return periods of HD duration would reduce by 0.28 and 1.54 months, respectively, compared with those in the historical stage, and they would decrease by 0.4 months and 2.2 months, respectively, under the RCP8.5 scenario. The severity of drought with the same return periods would also decrease by 0.15 and 1.19 months under the RCP4.5 scenario, and by 0.22 and 2.07 months under the RCP8.5 scenario, respectively. These imply that HD events in the JRB would be significantly reduced under projected future scenarios.

Return Periods of Hydrological Drought Events
Six return period levels from 3 to 100 (in Figure 6) for HD duration and severity are estimated according to Equations (7) and (8). Figure 6 shows that the values of the two variables for all return periods in both future scenarios would decrease significantly, especially in the RCP8.5 scenario. Under the RCP4.5 scenario, the 3-year and 100-year return periods of HD duration would reduce by 0.28 and 1.54 months, respectively, compared with those in the historical stage, and they would decrease by 0.4 months and 2.2 months, respectively, under the RCP8.5 scenario. The severity of drought with the same return periods would also decrease by 0.15 and 1.19 months under the RCP4.5 scenario, and by 0.22 and 2.07 months under the RCP8.5 scenario, respectively. These imply that HD events in the JRB would be significantly reduced under projected future scenarios. For the reason that different combinations of HD variables can lead to the same JRP, the JRPs for two variables in the historical period and future scenarios calculated by Equations (9) and (10) are shown in the contour plots in Figure 7. The JRPs at six levels (3, 5, 10, 20, 50, 100 years) for historical period and projected future scenarios are quite different. According to the contour plots, once certain values of HD variables are given (Table 5), the corresponding JRPs are much larger in future scenarios than that in the historical period. For example, when s is 2.8 and d is 5.5 in Figure 7a, the Tand for the historical period and projected RCP4.5 scenario are 10 and 20 years, respectively. This means that, with the same JRP, the future HD events would become less serious than those in historical situations. Table 5. Joint return periods for hydrological drought events at given duration and severity in Figure  6.  For the reason that different combinations of HD variables can lead to the same JRP, the JRPs for two variables in the historical period and future scenarios calculated by Equations (9) and (10) are shown in the contour plots in Figure 7. The JRPs at six levels (3, 5, 10, 20, 50, 100 years) for historical period and projected future scenarios are quite different. According to the contour plots, once certain values of HD variables are given (Table 5), the corresponding JRPs are much larger in future scenarios than that in the historical period. For example, when s is 2.8 and d is 5.5 in Figure 7a, the Tand for the historical period and projected RCP4.5 scenario are 10 and 20 years, respectively. This means that, with the same JRP, the future HD events would become less serious than those in historical situations. Table 5. Joint return periods for hydrological drought events at given duration and severity in Figure 6.

Water Scarcity in the Jinghe River Basin
The Jinghe River flows across three provinces of Ningxia, Gansu, and Shanxi. The annual renewable water resources produced by the JRB and the changes in water availability in major artificial reservoirs (for 2000-2017, in m 3 ) are obtained from the Water Resources Bulletins of the three provinces. The annual net water consumption can be evaluated according to Equation (13), and the annual WEI+ of the watershed can then be obtained based on Equation (8). The WEI+ of the JRB from 2000 to 2017 are shown in Table 6. The average annual net water consumption in 2000-2008 is 881 million m 3 and the average WEI+ is 0.59. The average annual net water consumption in 2009-2017 is 1093 million m 3 , and the average WEI+ is 0.62. Water consumption has increased significantly in these 9 years, while the increase in WEI+ is not visible, which may be due to precipitation increases during these years in the JRB.

Water Scarcity in the Jinghe River Basin
The Jinghe River flows across three provinces of Ningxia, Gansu, and Shanxi. The annual renewable water resources produced by the JRB and the changes in water availability in major artificial reservoirs (for 2000-2017, in m 3 ) are obtained from the Water Resources Bulletins of the three provinces. The annual net water consumption can be evaluated according to Equation (13), and the annual WEI+ of the watershed can then be obtained based on Equation (8). The WEI+ of the JRB from 2000 to 2017 are shown in Table 6. The average annual net water consumption in 2000-2008 is 881 million m 3 and the average WEI+ is 0.59. The average annual net water consumption in 2009-2017 is 1093 million m 3 , and the average WEI+ is 0.62. Water consumption has increased significantly in these 9 years, while the increase in WEI+ is not visible, which may be due to precipitation increases during these years in the JRB. With the development of our socio-economic society, water consumption is expected to increase continuously in the future [2]. To investigate possible conditions of WS in the future and explore the ultimate carrying capacity of water resources in the JRB, five water consumption scenarios for the future are set: the future water consumption is equivalent to the averaged water consumption in the past nine years (2009-2017), denoted as S0; the water consumption at the end of the century is 10%, 20%, 30% and 40% higher than that in the past 9 years and has a constant increasing rate for every ten years, which are respectively denoted as S1, S2, S3 and S4. The fluctuation of water reserves in the artificial reservoir is set to be zero since it has very little impact on the total water resources. The WEI+ values for the five water consumption scenarios under climate change projections of RCP4.5 and RCP8.5 are then calculated. The frequencies of WEI+ located in different intervals for each scenario are shown in Table 7.  Table 7 shows that the WEI+ would be most frequently located in the interval 0.7 to 0.8 under the five water consumption scenarios, with an average of 26.4 times under RCP4.5 and 16.6 times under RCP8.5. The values of WEI+ from S0 to S5 under the RCP4.5 scenario are mainly concentrated (occupied 73.5%) in the interval of [0.6, 0.9]. While the values of WEI+ are more dispersedly distributed between 0.4 and 0.9 under RCP8.5, indicating that the water pressure will vary significantly under this scenario. For the extreme scenario of S4, the WEI+ greater than 1 would be observed on six occasions under both climate change scenarios, which means that non-renewable water is available or water resources may need to be imported from outside to meet the demand in some years. The 10-year average WEI+ is applied to reflect the water pressure fluctuations at different time periods under different scenarios in the future ( Table 8). The values of WEI + at the end of the century under S3 and S4 scenario for RCP4.5 and S4 scenario for RCP8.5 are all greater than 0.8, which indicate that serious WS would occur at the end of the century. The WEI+ averages of S0-S4 in each time period show that the most severe WS under the RCP4.5 scenario is from 2059 to 2069, and the corresponding mean value of WEI+ is 0.81. According to the classification of water resources pressure shown in Table 1, severe WS occurs when WEI+ exceeds 0.4, and some researchers believe that the freshwater ecosystems of a river basin cannot be healthily maintained in this state [26,27]. However, some studies argue that the utilization rate of water resources can be greatly improved, and thus a threshold of 0.6 would be a more appropriate value [27]. Although it requires further in-depth investigation to set the threshold, according to the comprehensive results of this study, the JRB will still be in a state of extreme WS for a long time under the projected future scenarios, even with a significant increase in precipitation.

Conclusions and Discussions
The presented modeling framework in this study allows us to estimate the univariate and joint risks of hydrological drought variables, and also provide quantitative analyses for water scarcity under climate change. The WEI+ index was innovatively applied to evaluate water scarcity in an arid basin (JRB) of China under climate change. Lower frequency with lower mean duration and severity of HD events were identified under RCP8.5 scenario than those characterized under the RCP4.5 scenario. HD duration and severity under both RCP scenarios are projected to decrease remarkably when compared to those in the past. The copula-based bivariate HD risk assessment model also reveals that climate change under projected future scenarios would alleviate drought situation in the JRB. However, the WS assessment results show that the shortage of water resources in JRB is still severe in the future.
Drought indicators can flexibly evaluate the degree of water resources deviating from normal at different time scales, such as the monthly SDI estimated in this study. The developed modeling framework has important implications for HD risk management. However, drought assessment can only reflect the relative changes of water resources in a natural state, without taking into account water consumption changes in the future. Therefore, simply assessing the severity of future droughts cannot accurately reflect the severity of water shortages, and decision-makers cannot provide valuable water resources management policies. Considering the increase in human water consumption in the future, the degree of water shortage may not be alleviated by the increase in precipitation and water availability. This makes it necessary to introduce quantitative assessment methods for future water scarcity evaluation.
Although setting thresholds for the WEI+ is still debatable, the authors suggest that the WEI+ thresholds for the basin should be set based on different environmental protection objectives and economic development goals. Moreover, an obvious disadvantage of the WEI+ indicator is the lack of consideration of water quality, which is closely related to the quality of "return water" (exploited and re-entered into the hydrological cycle). The WEI+ may be improved or different WS assessment methods can be used to comprehensively assess water stress and water quality in the future. However, a river basin often covers several provinces, and the conditions of water use are complex, which requires the relevant departments to do a more collaborative job on basic data statistics. The results of this study are of great significance for policy makers and local stakeholders to make a prospective long-term plan. Based on the assessment of future drought and water shortage status, an appropriate plan can be made to support the sustainable development of social economy. This research is repeatable and can be used in other basins in traditional arid or semi-arid areas. In addition, due to the increase in population and the acceleration of urbanization, the water shortages in large cities and their surrounding areas in some developed countries are becoming increasingly severe, and some water-rich countries have also experienced WS events in recent years. The presented modelling system can also be applicable in such areas.