Climate Change Impacts on Extreme Flows Under IPCC RCP Scenarios in the Mountainous Kaidu Watershed, Tarim River Basin

Abstract

In the 21st century, heavier rainfall events and warmer temperatures in mountainous regions have significant impacts on hydrological processes and the occurrence of flood/drought extremes. Long-term modeling and peak flow detection of streamflow series are crucial in understanding the behavior of flood and drought. This study was conducted to analyze the impacts of future climate change on extreme flows in the Kaidu River Basin, northwestern China. The soil water assessment tool (SWAT) was used for hydrological modeling. The projected future precipitation and temperature under Intergovernmental Panel on Climate Change (IPCC) representative concentration pathway (RCP) scenarios were downscaled and used to drive the validated SWAT model. A generalized extreme value (GEV) distribution was employed to assess the probability distribution of flood events. The modeling results showed that the simulated discharge well matched the observed ones both in the calibration and validation periods. Comparing with the historical period, the ensemble with 15 general circulation models (GCMs) showed that the annual precipitation will increase by 7.9–16.1% in the future, and extreme precipitation events will increase in winter months. Future temperature will increase from 0.42 °C/10 a to 0.70 °C/10 a. However, with respect to the hydrological response to climate change, annual mean runoff will decrease by 21.5–40.0% under the mean conditions of the four RCP scenarios. A reduction in streamflow will occur in winter, while significantly increased discharge will occur from April to May. In addition, designed floods for return periods of five, 10 and 20 years in the future, as predicted by the GEV distribution, will decrease by 3–20% over the entire Kaidu watershed compared to those in the historical period. The results will be used to help local water resource management with hazard warning and flood control.

1. Introduction

The increased frequency and magnitude of extreme climate events under climate change has been one of the biggest challenges to human societies [1]. Many researchers have indicated that climate change will change hydrological cycles and have a severe impact on the availability of water resources and the occurrence of hydrological extremes [2]. Particularly in northwestern China, many rivers are recharged by snow melt [3]. Previous studies have shown that, hydrological processes of several main basins in Tianshan and Kunlun Mountains have been significantly influenced by climate change [4]. Extreme precipitation events have the potential to aggravate floods and droughts [5,6,7], while the continuous rise in air/surface temperature may lead to the increasing snow melting and evapotranspiration (ET) in such mountainous river basins [8]. Understanding variations in historical discharge extremes and evaluating future changes will greatly contribute to flood hazard control and water resource management in inland river basins [9,10].

Hydrological models play an important role in modern water resource management [11]. Climate change impacts on the hydrological cycle can be investigated by integrating projected future scenarios based on downscaled general circulation models (GCM) data in hydrological models [12]. Considerable uncertainties exist in climate change impacts on hydrological systems. Moreover, compared to the effects of climate change on average conditions, changes in extremes are much more uncertain [13]. The uncertainties reside primarily in GCM, GCM initial conditions, future emission scenarios, downscaling techniques, and the structure of hydrological models and their parameters. However, it is difficult to address these sources of uncertainty effectively, and it is particularly difficult to handle all of them at the same time through effective methods [14]. Kay et al. [15] found that the GCM structure is the most significant uncertainty source when compared to the uncertainties originating from emission scenarios, GCM initial conditions, downscaling techniques, hydrological model structures, and hydrological model parameters. Chen et al. [16] also showed that the choice of GCM is the major uncertainty contributor. Najafi et al. [17] indicated that the uncertainties associated hydrologic model was larger than GCM. Dams et al. [13] indicated that the uncertainty due to climate change scenarios is larger than the uncertainty introduced by the structure of hydrological models. To reduce the uncertainty in climate predictions [2], simplifications in GCM and downscaling methods [1,18,19] as well as ensemble methods of climate change scenarios have been commonly used to simplify the GCM outputs [20]. Several types of GCM ensemble have been used in recent studies [12], such as the arithmetic ensemble means (AEMs), the Bayesian model average (BMA), and the reliability ensemble average (REA) methods. Ntegeka et al. [20] described and used tailored climate change scenarios to limit the number of future scenarios. Hosseinzadehtalaei et al. [21] calculated the response of extreme precipitation to future climate change through sensitivity and uncertainty analysis on GCM, initial conditions of the GCM and representative concentration pathways (RCPs). However, the source and contribution of uncertainties in the context of climate change and the inherent impact on hydrology are regionally dependent, and more in-depth research is needed [22].

In China, most studies analyzing the climate change impact on extreme floods have focused on several major rivers in the southern and southeastern parts of the country [23]. Recently, the annual and seasonal changes of extreme hydrological events in inland rivers of northwestern China have attracted increasing attention as well. However, few studies have been performed to examine the impact of climate change in this region [24]. The Tarim River is located in northwestern China and is the longest inland river in the country [25]. Under the dual influence of climate change and human activities, the hydrological cycle in the headwater regions has been severely changed. Indeed, some upstream rivers have completely ceased flowing [26]. Currently, only four rivers, namely the Kaidu, Aksu, Hotan, and Yarkand Rivers, flow to the main stream of the Tarim River [27]. Among these source streams, water resources from the Kaidu River passing through Bosten Lake are the most important supplement of the lower Tarim River [11]. Therefore, the analysis and prediction of the runoff of the Kaidu River are essential for water resource management and disaster warnings [28].

Previous studies of the hydrological cycle in the Tarim River have mainly focused on fractal patterns of runoff, trends and change point detection based on annual or monthly streamflow data [29,30,31,32]. Few studies have been conducted on extreme hydrological events due to the lack of essential daily data sets [33]. Hydrological models can simulate the variability in runoff at different space-time scales and fill the data gaps of conventional observation stations, particularly at the daily or hourly scales [8]. Recent studies have evaluated the climate change impacts on the hydrology of the Tarim River headstreams using hydrological models. Liu et al. [34] investigated the impacts of climate change on the hydrology of three upstream catchments of the Tarim River using the variable infiltration capacity (VIC) model. Liu et al. [35] used both lumped and distributed hydrological models to assess the climate change impact on hydrology under emission scenarios (SERS) of the IPCC through 2050 s. Ma et al. [36] analyzed the response of snowmelt runoff to future climate in the Kaidu watershed by coupling GCM outputs with the snowmelt runoff model (SRM). Meng et al. [37] found that climatologically parameters are contributing to 92% of the runoff variability, and human activity was responsible for 7.72% of the runoff change in the Aksu watershed. Liu et al. [38,39] coupled an average ensemble of 18 available GCMs with a well-calibrated Mike System Hydrological European (MIKE SHE) model and simulated hydrological processes under future climate change scenarios. Their results indicated that snow storage at high altitudes will decrease due to increasing evaporation, which will therefore decrease the available water in the downstream region.

Through previous studies, the mechanism of climate change impacts on the hydrology of the Tarim headstreams has been gradually recognized [34,35,38,39,40]. However, research on extreme runoff prediction is still lacking [41]. Variations and tendencies of extreme flow events under historical and future climate change scenarios have not been investigated. Previous studies have shown that the frequency and intensity of extreme floods in the Kaidu River have increased over the past decades [28,42]. Therefore, the objective of this research is to assess the potential effects of future climate change on extreme floods in the upper Kaidu River and to use the advances noted above to enhance our understanding of future climate change and local hydrology variations. The SWAT model is used to simulate river runoff. The uncertainty in runoff projections related to GCMs and RCP scenarios are also investigated. Finally, the variation in flow extremes is revealed by a frequency analysis model. The obtained results will be used to help local water authorities establish effective flood control policies and thus improve water resource security in the region.

2. Study Area

The Kaidu watershed is located in the northeast Tarim River Basin and has an area of 1.90 × 10⁴ km² (Figure 1). The elevation of the Kaidu watershed ranges from 1400 m to 4794 m. The land cover is dominated by alpine meadow and grassland (81%), followed by surface water bodies (8%), bare rock (9.5%), and forest (1.5%). The soils in this region include fine and coarse sand, loamy sand, slit clay loam, sandy loam, loam, and rock [11].

In the study area, meteorological and hydrological observations began in the mid-1950s. There are two meteorological stations. One is the Bayblk (BYBLK) station, located in the upper reaches of the Kaidu River with an elevation of 2500 m, and the other is the Dashankou (DSK) station with an elevation of 1500 m, located at the outlet of the Kaidu watershed. The climate is typical of mountainous environments in arid regions, with sparse precipitation, low temperature and high evaporation. Average annual precipitation in the study area is approxmately 384 mm/year, and more than 80% of the total precipitation is distributed from May to September. The annual average temperature is approximately −4.16 °C, and the mean temperature in Summer is 9.8 °C. Pan evaporation is approximately 1157 mm/year, which is much higher than the precipitation in the region. And regional average actual ET is approximately 193 mm/year [43]. Snowmelt is the most important source of water in spring. Several studies have indicated that precipitation and temperature have increased significantly during the past 50 years, especially since 1990 [35].

Hydrological observation also conducted at the BYBLK and DSK stations. Drainage area controlled by the BYBLK station is 6.65 × 10³ km², and the drainage area controlled by the DSK station is 19.01 × 10³ km². There are two flood seasons in spring and summer in Kaidu River. Spring flow peaks occur in April to May due to the melting of snow. Summer flooding results from the combined output of rainfall and snowmelt in high-altitude regions [44]. The river flood events not only represent destructive natural hazards in the mountainous region but also affect the lower reaches of the Tarim River [28]. Therefore, runoff prediction and flood risk assessment are essential for water resource management and flood control in the study area.

3. Data and Methodology

3.1. Data Availability

A digital elevation model (DEM) with 90 × 90 m resolution, land use/cover, and soil types were used to set up the hydrological model. Land use and land cover map in 2010 based on Landsat Thematic Mapper (TM) images were processed by the Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences (XIGE) [26]. Four land cover types (i.e., meadow, forest, rock, and surface water body) were detected in the Kaidu watershed and used for hydrological simulations. The soil in the study area was classified into seven types. Related physical parameters for each soil type were defined according to the FAO soil property reference [39]. The initial soil physical parameters (e.g., hydraulic conductivity, available water capacity) for each soil type were obtained by using the SPAW software developed by USDA (Saxton and Rawls, 2006).

Daily discharges at the DSK and BYBLK hydrological stations were collected from 1996 to 2011. Daily air temperatures, precipitation, humidity, and wind speed at the BYBLK and DSK stations from 1961 to 2011 were collected from the Tarim River Management Bureau (TRMB).

The daily precipitation and temperature outputs of 32 GCMs were downloaded from the dataset of Phase 5 of the Coupled Model Inter-comparison Project (CMIP5) [35,39].

3.2. Methodology

3.2.1. Hydrological Modeling

The SWAT model is a physically distributed hydrological model that has been widely used to assess the impacts of climate change and human activities on water [37,45]. To simplify the watershed and facilitate calculation, the study area was divided into 22 subbasins with 271 hydrologic response units (HRUs) based on land use, soil type, and slope.

$${SC}_{ti}={SC}_0+{\displaystyle \sum }_{i=1}^t(P_{da{y}_i}-Q_{sur{f}_i}-E_{a_i}-V_{see{p}_i}-Q_{g{w}_i})$$

(1)

where $S{C}_{ti}$ is the soil water content (mm) at time t, $S{C}_0$ is the initial water content (mm), t is the simulation period (days), $P_{da{y}_i}$ is the precipitation on the ith day (mm), $Q_{sur{f}_i}$ is the surface runoff on the ith day (mm), $E_{a_i}$ is ET on the ith day (mm), $V_{see{p}_i}$ is the water entering the vadose zone from the soil profile on the ith day (mm), and $Q_{g{w}_i}$ is the base flow on the ith day (mm).

Surface runoff was simulated using the soil conservation service (SCS) curve numbers (USDA-SCS, 1972) in SWAT. Redistribution of water between soil layers was calculated using a kinematic storage model. The underground flow was calculated as a function of the saturated hydraulic conductivity. The Muskingum routing method was applied to calculate the channel flow routing, and the model used a degree-day approach to estimate snow accumulation and melting.

The FAO Penman–Monteith method and Hargreaves method were used to estimate ET. Daily precipitation, temperature, wind speed, and relative humidity data were obtained from the BYBLK and DSK stations during the historical period (i.e., 1960–2010), and the solar radiation inputs were generated using the built-in SWAT stochastic weather generator based on data from nearby weather gauges. PET and ET, calculated according to the Penman–Monteith method, were close to those reported in the literature [43]. However, only temperature can be used to calculate PET in future periods because of the higher uncertainties in the GCM outputs for the additional variables. Therefore, we used PET calculated via the Penman–Monteith method as a standard value and calibrated the parameters of the Hargreaves method, which is only based on temperature and geographical location [46]. The calibrated Hargreaves method also showed reasonable PET results in the study area. Therefore, the PET data were calculated according to the calibrated Hargreaves method for both historical and future periods.

The sequential uncertainty fitting (SUFI-2) program version 2 embedded in the SWATCUP software [47] was used for sensitivity analysis and model calibration. The T-states and p-values of the parameters were used before and after improving the statistical model to obtain parameter sensitivity [48,49]. The period from 1996 to1997 was chosen as the spin-up period, while calibration and validation periods were from 1998 to 2007 and 2008 to 2010, respectively. Daily discharge data from the BYBLK and DSK stations were applied for model validation. The Nash–Sutcliffe coefficient (NS), water balance bias (WB) and correlation coefficient (R) were applied to evaluate the simulation accuracy [50].

$$NS=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(q_{o,i}-q_{s,i}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(q_{o,i}-{\overline{q}}_o\right)}^2}$$

(2)

$$WB=100\%\times (1-\frac{{{\displaystyle \sum }}_{i=1}^n\left|q_{o,i}-q_{s,i}\right|}{{{\displaystyle \sum }}_{i=1}^n{q}_{o,i}}$$

(3)

$$R=\frac{{{\displaystyle \sum }}_{i=1}^n\left(q_{o,i}-{\overline{q}}_o\right)\left(q_{s,i}-{\overline{q}}_s\right)}{{{\displaystyle \sum }}_{i=1}^n{\left[{\left(q_{o,i}-{\overline{q}}_o\right)}^2\right]}^{\frac{1}{2}}{{\displaystyle \sum }}_{i=1}^n{\left[{\left(q_{s,i}-{\overline{q}}_s\right)}^2\right]}^{\frac{1}{2}}}$$

(4)

where $q_{o,i}$ is the observed discharge at the time step $i$, $q_{s,i}$ is modeled discharge at the time step $i$, ${\overline{q}}_o$ is the mean observed discharge, and n is the total number of time steps. With the calibrated parameters, validation was performed using daily discharge data from the BYBLK and DSK stations.

3.2.2. Downscaling of Future Climate Scenarios

The global climate is a complex system; therefore, when global climate and regional models (e.g., GCMs and RCMs) are used to model the climatic system, considerable simplification must be conducted in the description of physical processes [51]. To incorporate the uncertainty of social actions, a range of future greenhouse gases (GHG) scenarios is taken into consideration [12]. Multimodel ensemble simulations using RCMs and GCMs have been shown to outperform individual models and can be used to more accurately model a given property through the provision of a large sample size, e.g., of a climatological mean of the frequency of a rare event [1].

Monthly precipitation totals and temperature were selected for intercomparison and assessment of the performance of GCMs. For average annual precipitation, monthly precipitation, and mean monthly temperature, the similarities between the GCM and observational data are quantified in terms of the root-mean-square relative error (RMSRE) and the relative bias (BIAS) which detail described in Liu et al. [35]. Further, they are calculated by using the observed climate data at BYBLK station as the referenced baseline. The period ranging from 1961 to 2000 was considered for the GCM control runs (CN) as well as for the observation data at the BYBLK station. In general, the GCMs present acceptable results; the mean annual temperature BIAS is 12.2%, and the RMSREs are less than 6%. For precipitation, BIAS ranges from −36% to 98% at the yearly scale, and the RMSREs of annual precipitation exceeds 18%. The climate models are less reliable and consistent in predicting precipitation variations than in describing temperature variations at the monthly scale.

In addition, extreme precipitation events may be the main cause of extreme flow, especially in Summer. Therefore, we performed a comparison of the probability distributions to assess the performance of the GCMs in simulating precipitation events [36,37]. For each GCM runs, extreme precipitation sequence were extracted based daily precipitation series, and when the precipitation bigger than 0.1 mm/day were selected as a wet day [35]. The cumulative distributed functions of extreme precipitations are shown in Figure 2. The black line is calculated from observed data from 1961 to 2000 at BYBLK station, and the grey lines are calculated from the selected GCM extreme precipitation series. After examining the consistency between the GCM monthly total precipitation and observation data at the BYBLK station, and also a comparison of the probability distributions of extreme precipitation events, 15 GCMs with a total of 31 runs (i.e., CNRM-CM5_r1i1p1, EC-EARTH_12i1p1, EC-EARTH_r8i1p1, MIROC-ESM_r1i1p1, MIROC-ESM-CHEM_r1i1p1, MPI-ESM-LR_r1i1p1, MPI-ESM-LR_r2i1p1, MPI-ESM-MR_r1i1p1, MRI-CGCM3_r1i1p1, bcc-csm1-1_r1i1p1, bcc-csm1-1-m_r1i1p1, BNU-ESM_r1i1p1, CanESM2_r1i1p1, CanESM2_r2i1p1, CanESM2_r3i1p1, CanESM2_r4i1p1, CanESM2_r5i1p1, CCSM4_r1i1p1, CCSM4_r2i1p1, CESM1-CAM5_r1i1p1, FGOALS-g2_r1i1p1, GFDL-CM3_r1i1p1, GFDL-ESM2G_r1i1p1, IPSL-CM5A-LR_r1i1p1, IPSL-CM5A-LR_r2i1p1, IPSL-CM5A-LR_r3i1p1, IPSL-CM5A-LR_r4i1p1, IPSL-CM5A-MR_r1i1p1, MIROC5_r1i1p1, MIROC5_r2i1p1), were acceptable. In general, it is an ensemble of the prediction results from different GCM (and the results from the same GCM with different initial conditions, initialization method and perturbed physics [52].

Changes in precipitation and temperature under the future climate were extracted from the simulation results of the GCMs and applied to the observed records over the historical period at the BYBLK station. For temperature, a simpler delta approach was used, and the change in temperature was described as the absolute difference between the future and control periods [38]. The daily temperature in the future period can be calculated as follows:

$$T_{DO\left(k,i\right)}={\mathrm{T}}_{O\left(k,i\right)}+\left(T_{F\left(k\right)}-T_{H\left(k\right)}\right)$$

(5)

where $T_{DO\left(k,i\right)}$ is the future temperature on the ith day in the kth month, $T_{O\left(k,i\right)}$ is the observed temperature on the ith day in the kth month, and $T_{F\left(k\right)}$ and $T_{H\left(k\right)}$ are the mean temperatures predicted by the GCMs in the kth month for the periods from 2041 to 2080 and from 1961 to 2000, respectively.

A large deviation was noted in precipitation predicted by the GCMs. A quantile perturbation method (QPM) which proposed by Willems et al. [53] and improved by Liu et al. [38] was applied to generate the daily precipitation series under future climate change scenarios in the period of 2041 to 2080. Detail description of the modified QPM can be found in Liu et al. [38].

The uncertainties of the projected precipitation and temperatures have been accounted for in this study by making use of an ensemble modeling approach [20]. In relation to the baseline (i.e., 1961–2000), all change signals of precipitation and temperature in future period by each GCM runs were extracted statistically. Changes under each emission scenario and for all considered runs that passed the performance evaluation were clustered into low, mean, and high levels before input into the SWAT model. The low and high levels correspond to the minimum and maximum values of the changes, respectively, and expresses the ranges of climate change uncertainty. The mean level corresponds to the mean values of the changes. For hydrological impact simulation, the average tendency of multiple GCMs was recommended [39]. Modeling results for daily maximum and minimum temperature and mean precipitation from a multi-model-ensemble of 15 GCM combinations for four GHG concentration pathways were employed to provide daily time series for the SWAT model.

3.2.3. Frequency Analysis of Extreme Floods

Fisher and Tipett [54] presented three extreme value distributions: the Gumbel distribution, Fréchet distribution and Weibull distribution. Jenkinson [55] and Coles [56] improved the three extreme value distributions as a three-parameter extreme value distribution, called the generalized extreme value (GEV) distribution. The probability density function (PDF) was calculated as follows:

$$F\left(x\right)=\exp \left\{-{\left[1+\epsilon \left(\frac{x-\mu }{\sigma }\right)\right]}^{-\frac{1}{\epsilon }}\right\}$$

(6)

$$\left[1+\epsilon \left(\frac{x-\mu }{\sigma }\right)\right]>0$$

(7)

where $\sigma$ is the scale factor, $\mu$ is the location factor, and $\epsilon$ is the shape factor. The GEV distribution includes three distributions. When $\epsilon$ = 0, $\epsilon$ > 0, and $\epsilon$ < 0, the GEV distribution could be described as a Gumbel distribution, Frèchet distribution, and Weibull distribution, respectively. The advance of the GEV distribution is that the unification of the three distributions can avoid shortcomings of a single distribution.

For parameter estimation of the extreme distributions, maximum likelihood estimation (MLE) was used. Compared with other parameter estimation methods, MLE can be applied for every population, has good asymptotic behavior in the case of large sample, and obtains consistent and effective parameters [57]. Suppose that {x1, x2, …, xn} is independent and identically distributed with the probability distribution function F(x), the parameter estimation of the GEV distribution using the MLE method can be obtained from the following log-likelihood function:

$$L\left(\theta \right)=L\left(\mu ,\sigma ,\epsilon \right)=-n\ln \sigma -{\displaystyle \sum }_{i=1}^n{\left[1+\epsilon \left(\frac{x_i-\mu }{\sigma }\right)\right]}^{-\frac{1}{\epsilon }}-\left(1+\frac{1}{\epsilon }\right){\displaystyle \sum }_{i=1}^n\ln \left[1+\epsilon \left(\frac{x_i-\mu }{\sigma }\right)\right]$$

(8)

where $\mathsf{\theta }=\left(\mu ,\sigma ,\epsilon \right)$; point $\left(\mu ,\sigma ,\epsilon \right)$ reaches the MLE when the function reaches the maximum point. In this study, numerical methods were used to solve the function.

The return period refers to the number of years between two floods, which is a key parameter for studying flood problems. It is a measure of safety (the ‘inverse risk’). In combination with the method of periodic maxima, the return period (PR) is calculated as the inverse of the population survival distribution of the annual maxima:

$$RP=\frac{1}{1-F\left(x\right)}$$

(9)

Figure 3 shows the comprehensive method by coupling the distributed hydrological model (SWAT) and GCMs with a statistical frequency analysis.

4. Results and Discussion

4.1. Calibration and Validation of the Hydrological Model

Totally 27 parameters were selected to participate in the calibration of the model. In order to obtain the sensitivity information of each parameter in the model, the model was simulated 1000 times to obtain the sensitivity information of the parameters. The results indicate that there were several parameters that show high sensitivity, including the effective hydraulic conductivity in main channel alluvium (CH_K2), groundwater delay (days) (GW_DELAY), snow pack temperature lag factor (TIMP), saturated hydraulic conductivity (SOL_K) and temperature lapse rate (TLAPS). The T-states of the CH_K2, GW_DELAY, TIMP, SOL_K, and TLAPS are 26.6, 13.7, 6.66, 5.79 and 2.35, respectively. Moreover, the p value of the sensitive parameters are below 0.1 with the great significance. The modeled results match the observations well (Figure 4). The evaluation statistics for model performance are listed in

Table 1. SWAT model performance for daily and monthly runoff.

		Calibration (1998–2006)			Validation (2007–2011)
		EF	R	RE (%)	EF	R	RE (%)
Monthly	DSK	0.86	0.95	3.00	0.85	0.92	3.7
	BYBLK	0.74	0.89	−6.32	0.71	0.86	5.88
Daily	DSK	0.76	0.91	3.44	0.72	0.90	3.07
	BYBLK	0.62	0.83	−5.12	0.61	0.77	4.66

. During the calibration period, the NS values are 0.76 and 0.62 at the DSK and BYBLK hydrological stations, respectively, the R values at the DSK and BYBLK stations reach 0.91 and 0.83, respectively, and the WB varies between 3.44% and −5.12% at the DSK and BYBLK stations, respectively. For the monthly scale, the model performances increase at both stations. The timing and volume of the simulated flow peaks are slightly shifting and higher than the observed values. Following the guidelines by Moriasi et al. [58], these simulation results can be considered acceptable. However, the NS value for daily discharge simulation at the BYBLK station is a litter lower than the references value (>0.65) according to Ritter and Muñoz-Carpena [59].

4.2. Predicted Changes in Temperature and Precipitation

Under future climate change scenarios, both the magnitude and frequency of the precipitation would change. The perturbation factor is the changes of precipitation or temperature time series between the future and control period. For temperature, the perturbation factors defined as absolute differences between future and control periods; for precipitation, they are defined as relative changes (ratios of the values during future versus control periods). Figure 5 shows the mean monthly frequency changes determined by the GCMs for different RCPs. The results indicate that, under the mean level, the variation of rainy days in each month is the largest under the RCP2.6 scenario, ranging from −8% to 6%; and the most moderate under RCP 6.0 scenario, ranging from −1% to 2%. Under the RCP4.5 and RCP8.5 scenarios, the change of rainfall days was found similar. The number of rainy days would increase under RCP4.5 and RCP8.5 from September to the following January, and also would increase from March to June; in February, July, and August, the number of rainy days would exhibit a moderate declining of approximately 3%. Under the low and high levels, the changes of rainy day numbers would be from −20% to −3% and 6 to 22%, respectively.

Figure 6 shows the precipitation data for the 0.01, 0.05, 0.1, 0.2, and 0.5 quantiles, which were selected to illustrate the monthly quantile perturbations of precipitation intensity in the future period relative to those of the history period. The results indicate that the rain intensities for different quantiles are not obviously changed in summer; the changes in winter are significant. The intensity of extreme precipitation (at the 0.01 quantile) in January under RCP4.5 is most significant, with an average increase of 26%.

The average monthly precipitations in the historical period and future climate change scenarios (under low, mean, and high levels) show notable variations (Figure 7). For the mean level, the future precipitation in the Kaidu watershed exhibits an increasing trend under all RCP scenarios, while monthly precipitation presents a decreasing trend in spring (i.e., March, April, and May), autumn (i.e., September, October and November) and winter (i.e., December, January, and February) and is almost unchanged in summer (i.e., June, July, and August). When comparing the changes among the RCP scenarios, the precipitation shows the largest change in the RCP8.5 scenario, while the RCP2.6 and RCP6.0 scenarios keep similar rates of change. The mean changes in monthly precipitation are from 2.9% to 19.2%, −3.0% to 26.1%, 0.2% to 19.0%, and −3.4% to 35% under RCP2.6, RCP4.5, RCP6.0, and RCP8.5, respectively.

The temperature changes are positive and range from 0.2 to 0.6 °C/10 a at the BYBLK station (Figure 8). The average monthly changes (mean level) range from 0.37 to 0.51 °C/10 a, 0.46 to 0.63 °C/10 a, 0.47 to 0.68 °C/10 a, and 0.64 to 0.85 °C/10 a under the RCP2.6, RCP 4.5, RCP6.0, and RCP8.5 scenarios, respectively. Notably, temperature will rise more in winter and summer than in spring and autumn.

4.3. Climate Change Impact on Hydrological Processes

4.3.1. Water Balance

Annual streamflow in the Kaidu watershed ranges from 2.77 × 10⁹ m³ to 3.19 × 10⁹ m³ per year under future climate change scenarios. At both the DSK and BYBLK stations, under all RCP scenarios, annual mean runoff shows a clear decrease from 2041 to 2080. Generally, the annual runoff would change most significantly under the RCP8.5 scenario, stream runoff would decreases by 40% and 22% at the DSK and BYBLK stations, respectively. Minimal runoff changes would be found under the RCP2.6 scenario, annual mean runoff decreased by 21.1% and 14.7% at the DSK and BYBLK stations, respectively.

Such hydrological variation is mainly related to the relative variation between evapotranspiration and precipitation. Xu et al. [12] found that small variations in precipitation may lead to large changes in streamflow from upstream rivers. In Liu et al. [38]’s research on the Yarkant watershed, which is one of the other upstream catchments of the Tarim Basin, the amount of water dissipated by evaporation in the basin exceeds the amount of snow storage and streamflow to the lower reaches. Figure 9 presents the relative variations in annual average precipitation, ET and surface runoff at the subwatershed monitored by the BYBLK station and the entire Kaidu watershed, which is monitored by the DSK station. Generally, the changes in the water balance are more significant under the RCP4.5 and RCP8.5 scenarios and are relatively smaller under the RCP2.6 and RCP6.0 scenarios (Figure 9). The climate impact on annual mean precipitation is not substantial. The relative variations in annual precipitation range from 0.57% to 2.42% in the Kaidu watershed and from 0.60% to 2.49% in the area monitored by the BYBLK station. However, the annual mean ET shows a large increasing trend. The relative changes in ET range from 22.7% to 33.0% over the entire watershed and from 3.1% to 22.8% at the BYBLK subwatershed level. Therefore, water resources at both subwatershed and watershed scales will decrease.

Notably, there are large uncertainties in runoff assessment due to the uncertainties in GCM emission scenarios and the structure and parameters of hydrological models [60]. Zhang et al. [41] found that the uncertainties due to emission scenarios, parameters of the hydrological model, and the extreme flow distributions might be much smaller than the GCM uncertainty. In this study, only uncertainties from GCMs and emission scenarios were investigated, and the results proved that the uncertainty originating from the GCMs is larger than that related to the different emission scenarios. Moreover, when compared with the previous study by Liu et al. [35], large uncertainties from ET calculations are found due to different structures between the MIKE SHE and SWAT models, although annual precipitation and temperature showed similar relative changes under the IPCC AR4 and AR5 emission scenarios. With similar changes in precipitation and temperature, ET changed from 5% to 10% in Liu et al. [35]’s study. However, in this study, the change in ET reached 20% to 30%, and the mean ET ranges from 182.1 to 266.5 mm/year, which agrees with the observed ET in the Kaidu watershed [43].

4.3.2. Monthly River Flows

Since the seasonal variation in discharge is essential for flood control and water assignment, the monthly change in discharge was reviewed at both stations. For the DSK station, monthly river flows decrease in most simulated scenarios (Figure 10a). Significant reductions can be found in winter (i.e., December, January, and February), autumn (i.e., September, October, and November), and the summer (i.e., June, June, and August). The largest decrease in discharge reaches 50.3% under RCP8.5. Increases in simulated monthly runoff can be found during the snow melt period (i.e., April and May). The largest increase in monthly discharge is as high as 7.4% under the RCP8.5 scenario, and the increase of monthly discharge also found in May under RCP2.6 and RCP4.5 scenarios with the change rates of 6.5% and 5.3%.

Variations in the monthly discharge at the BYBLK station differ from the downstream DSK station (Figure 10b). Generally, monthly streamflow increases under most simulations. Moreover, the BYBLK station has relatively higher increases than the DSK station. Monthly runoff would increase from 3.0% to 60.5% in winter and autumn, and from 27.5% to 78.3% in summer. During the snow melting period, positive impacts (from 8.3% to 59.0%) occur primarily in March and April, while negative impacts (from 2.7% to 20.9%) are found in May and June.

Variations in seasonal distribution of river runoff are mainly due to the relative change in snow accumulation and melting. Figure 11a shown that the climate change impact on annual mean snowfall decreases in the future period due to significant temperature increases. However, the relative changes over the entire Kaidu watershed are larger than that in the BYBLK subwatershed. Mean annual snow storage decreases continuously in the Kaidu watershed primarily because of the decrease in snowfall resulting from the rising temperatures. The decrease occurs mainly in summer, from June to September. Compare with the snowfall from 1961 to 2000, the changes in annual snowfall in the Kaidu watershed would be −19.2%, −21.4%, −24.2%, and −26.4% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5 scenarios, respectively, and those in the BYBLK subwatershed would be −17.7, −19.4%, −22.3%, and −23.8% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5 scenarios, respectively. Thus, despite the higher temperatures, precipitation in the form of snow continues to occur, although less frequently, in winter.

Increased temperatures lead to earlier occurrence of snow melting, indicating an increasing trend in monthly streamflow in spring; the monthly mean streamflow in June, July, and August decreases because of the reduction in rainfall and relatively larger ET increments in summer. Drastically decreasing snowfall and earlier snowmelt results in less snow available for melting in summer and autumn months. As shown in Figure 11b, the mean annual snowmelt declines sharply under the future climate change scenarios. The changes in annual snow melt in the Kaidu watershed are from −18.69%, −20.5%, −23.4%, and −25.2% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5 scenarios, respectively, and are from −17.8%, −19.6%, −22.5%, and −24.0% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5 scenarios in the BYBLK subwatershed, respectively. Similar results were found by Liu et al. [34,35].

4.3.3. Extreme Flows

As shown in Figure 12, the SWAT model predicts that the peak flows would almost decrease under all RCP scenarios. The changes in extreme discharge are roughly from −26.6% to −2.5% at the DSK station and from −36.5% to 1.8% at the BYBLK station. At the BYBLK station, the model detects only a positive change in the peaks with a return period of 20 years under RCP6.0. Generally, changes in extreme flow are largest under the RCP4.5 scenario and show smaller changes under the RCP2.6 and RCP6.0 scenarios.

Figure 13 shows the changes in monthly maximum daily discharge at Figure 13a the DSK station and Figure 13b the BYBLK station under future climate change. The results indicated that, at the DSK station, monthly maximum daily discharge would increase in April, July, August, September, and October; the largest discharge increment can be found in April, and the change rate would be 56%, 97.0%, 88.8%, and 97.7% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5, respectively. At the BYBLK station, the increase of monthly maximum daily discharge is mainly found in March, April and October, and the change rate ranges from 15.3% to 31.0%, 30.2% to 53.1%, and 17.2% to 22.6%, respectively. The monthly extreme flow would decrease in winter, late spring, and early summer months at both stations, and the most significant negative change would reach 30% at the DSK station under the RCP8.5 scenario.

In flood control operations, the frequency and magnitude of peak flows are the most important details for designing a dam or reservoir [7]. Therefore, time series of the annual maximum discharge were generated to describe the extreme values of the GEV distribution. Subsequently, the frequency of extreme flows was analyzed.

Table 2. Parameter estimation, the type of GEV distribution and related statistical values of the peak flow series.

Gauge Station	Maximum Daily Runoff	Main Parameters			Fitted Distribution	Statistical Value
		$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\epsilon }$		K-S	R
DSK	Base line	339.51	97.33	0	Gumbel	0.112	0.92
	RCP2.6	290.17	69.63	0.240	Frèchet	0.066	0.97
	RCP4.5	261.43	88.41	0	Gumbel	0.195	0.94
	RCP6.0	279.22	73.15	0.261	Frèchet	0.108	0.87
	RCP8.5	264.20	90.10	0	Gumbel	0.189	0.95
BYBLK	Base line	105.90	46.44	0	Gumbel	0.272	0.93
	RCP2.6	104.53	28.12	0.269	Frèchet	0.053	0.96
	RCP4.5	84.78	40.55	0	Gumbel	0.201	0.87
	RCP6.0	90.39	31.74	0.360	Frèchet	0.159	0.88
	RCP8.5	103.27	27.99	0.275	Frèchet	0.106	0.87

presents the values of three parameters of the GEV distributions, the type of the GEV distributions and related statistical values of the peak flow series. At the DSK station, time series of peak flow are represented by a Gumbel distribution (type I of the GEV distribution) under the RCP4.5 and RCP8.5 scenarios, and can be described by a Frèche distribution (type II of the GEV distribution) under the RCP2.6 and RCP6.0 scenarios. For the BYBLK station, peak flow series are fitted by a Gumbel distribution under the RCP4.5 scenario, while the series match a Frèche distribution under the RCP2.6, RCP6.0 and RCP85 scenarios.

The designed floods at various return periods under the four RCP scenarios show distinct differences (

Table 3. Comparison of design floods (m³/s) under historical and future climate change scenarios.

Gauge Station	Scenario	Return Period (Years)
		5	10	20	50	100
DSK	Base line	482	560	643	738	800
	RCP2.6	412	503	596	713	739
	RCP4.5	387	455	521	593	667
	RCP6.0	406	513	599	721	810
	RCP8.5	376	466	527	589	679
BYBLK	Base line	175	206	250	286	316
	RCP2.6	152	191	229	269	310
	RCP4.5	139	172	208	233	279
	RCP6.0	157	203	251	325	386
	RCP8.5	155	200	233	276	328

). Based on historical records, the 5-, 10-, 20-, 50- and 100-year floods are 482 m³/s, 560 m³/s, 643 m³/s, 738 m³/s and 800 m³/s at the DSK station, respectively, and are 175 m³/s, 206 m³/s, 250 m³/s, 286 m³/s and 316 m³/s at the BYBLK station. However, variations in the design floods will be obviously different under future climate change scenarios. Generally, the future design flood would decreases from approximately 1% to 22% at both stations. However, at the BYBLK station, the design flood with 50- and 100-year return periods would increase under the RCP2.6 and RCP6.0 scenarios; the design flood with a 100-year return period at the DSK station would increase under the RCP6.0 scenarios. The changes in the design floods under the RCP4.5 and 8.5 scenarios are higher than those under the RCP2.6 and 6.0 scenarios.

5. Conclusions

This study investigated climate change impacts on future streamflow and peak flows using SWAT and GEV statistical analysis in the Kaidu watershed. SWAT was calibrated for the Kaidu watershed and the BYBLK subwatershed and reproduced the daily and monthly runoff reasonably well. Due to the uncertainty among GCMs, the average ensemble of 15 GCMs was employed to assess the hydrological impact under future climate change.

In the Kaidu watershed, the numbers of rainy days would be increased in most winter and spring months, and would be decreased in summer months under all RCP scenarios. The changes of monthly rainy days under the RCP2.6 scenario may have the largest range from −6% to 5% at the mean level. The rain intensities for different quantiles would not obviously change in summer, while the changes in winter would be significant. Mean annual precipitation in the future would exhibit an increasing trend under the four RCP scenarios. In addition, the mean changes in precipitation would be from 1.2% to 50%, 0.5% to 61%, 0.2% to 46%, and 1.0% to 78% under the RCP2.6, RCP4.5, RCP6.0, and RCP8.5 scenarios, respectively. The temperature changes would be positive, and the mean monthly changes would range from 1.86 to 2.57, 2.31 to 3.16, 2.35 to 3.49, and 3.20 to 4.23 degrees under the RCP2.6, RCP 4.5, RCP6.0, and RCP8.5 scenarios, respectively.

Under all RCP scenarios, annual mean runoff would decrease, and the rates of decrease are largest for RCP8.5 (i.e., 40%). A significant reduction in streamflow can be found in winter (i.e., December, January, and February), autumn (i.e., October and November) and the beginning of summer (i.e., June). The largest decrease in monthly discharge would reach 50.3% under the RCP8.5. Increases in simulated monthly runoff will be found in early April to May, and discharge increases in April will be obviously higher than those in other months. The changes in annual maximum daily discharge are roughly from −26.6% to −2.5%; monthly peak flow would increase significantly in early spring, and would decrease in winter, late spring, and early summer months. In addition, design floods in the future (i.e., 2041–2080), which were simulated by the SWAT model with downscaled historical GCM data, would decrease by 3% to 20% compared to those in the historical period (i.e., 1961–2000). Generally, under the impact of future climate change, the flood intensity in the study area would decrease and the extent of the drought would intensify.