Evaluating climate change impacts on streamflow variability based on a multisite multivariate GCM downscaling method in the Jing River of China

Projected hydrological variability is important for future resource and hazard management of water supplies because changes in hydrological variability can cause more disasters than changes in the mean state. However, climate change scenarios downscaled from Earth System Models (ESMs) at single sites cannot meet the requirements of distributed hydrologic models for simulating hydrological variability. This study developed multisite multivariate climate change scenarios via three steps: (i) spatial downscaling of ESMs using a transfer function method, (ii) temporal downscaling of ESMs using a single-site weather generator, and (iii) reconstruction of spatiotemporal correlations using a distribution-free shuffle procedure. Multisite precipitation and temperature change scenarios for 2011–2040 were generated from five ESMs under four representative concentration pathways to project changes in streamflow variability using the Soil and Water Assessment Tool (SWAT) for the Jing River, China. The correlation reconstruction method performed realistically for intersite and intervariable correlation reproduction and hydrological modeling. The SWAT model was found to be well calibrated with monthly streamflow with a model efficiency coefficient of 0.78. It was projected that the annual mean precipitation would not change, while the mean maximum and minimum temperatures would increase significantly by 1.6± 0.3 and 1.3± 0.2 C; the variance ratios of 2011–2040 to 1961–2005 were 1.15± 0.13 for precipitation, 1.15± 0.14 for mean maximum temperature, and 1.04± 0.10 for mean minimum temperature. A warmer climate was predicted for the flood season, while the dry season was projected to become wetter and warmer; the findings indicated that the intra-annual and interannual variations in the future climate would be greater than in the current climate. The total annual streamflow was found to change insignificantly but its variance ratios of 2011–2040 to 1961–2005 increased by 1.25± 0.55. Streamflow variability was predicted to become greater over most months on the seasonal scale because of the increased monthly maximum streamflow and decreased monthly minimum streamflow. The increase in streamflow variability was attributed mainly to larger positive contributions from increased precipitation variances rather than negative contributions from increased mean temperatures.


Introduction
In comparison with changes in the mean state, hydrological variability can cause more disasters such as flooding or drought and seriously threaten natural and social systems.Under the background of global warming, the expectation is that hydrological extremes would occur more frequently with greater severity because of changes in climate extremes.However, this has not been proven conclusively because of limited surface observations (IPCC, 2013), complex watershed properties (Andrés-Doménech et al., 2015), and lim-itations both in hydrological modeling and in the development of climate change scenarios (Wilby and Harris, 2006).Therefore, the extent to which hydrological variability is influenced by climate variability should be investigated thoroughly.
According to global-scale projections, hydrological variability will not change uniformly across the globe.For example, it is predicted that 30-year floods will occur more frequently over 50 % of the globe (Dankers et al., 2014) and that increased hydrological droughts will occur over 40 % of the analyzed land area (Prudhomme et al., 2014).In addition to the effects of catchment properties, the spatial variations of hydrological variability changes are due to those of climaterelated changes.For example, drought increase is generally located where precipitation decreases; however, drought can still increase in some areas with increased precipitation if stronger evaporation is driven by temperature increase (Prudhomme et al., 2014).Thus, the mechanisms by which climate variability influences hydrological variability should be analyzed.
The greatest uncertainties in impact assessments of climate change on hydrology originate from climate change scenario development, including general circulation models (GCMs), emission scenarios, and downscaling methods (Wilby and Harris, 2006;Kingston and Taylor, 2010;Chen et al., 2011).These uncertainties, to some extent, can be interpreted by introducing multiple climate models, emission scenarios, and downscaling methods.However, although the reconstruction of the spatial structure of climate has been considered widely in downscaling techniques, it has been incorporated rarely in hydrological modeling.When hydrological models are applied to large-scale catchments with multisite climate without consideration of the spatial structure, high flow in one subbasin could be offset by low flow in a neighboring subbasin (Wilks, 1998;Thyer and Kuczera, 2003;Clark et al., 2004;Wheater et al., 2005).Therefore, failure to feed hydrological models with appropriate spatiotemporally correlated climate could reduce hydrological variability and potentially misrepresent climate risks.
Parametric weather generators, including the Richardsontype model and the circulation-based model (Katz and Parlange, 1996), are used widely for hydrological modeling because of their simple implementation.The Richardson-type model directly perturbs the daily weather generator parameters based on changes in the corresponding monthly statistics (Wilks, 1999a), and the circulation-based model specifies the daily variation of parameters using regressions between local predictands and large-scale predictors (Wilby et al., 2002).After parametric adjustment, the two models are driven with correlated random numbers to capture the spatial structure exhibited in the daily weather data.With its direct parametric adjustments, the Richardson-type model in particular could be used more fruitfully for either impact assessments or sensitivity analyses.Specifically, it could be used either within a "top-down" framework of an impact study by downscaling GCM outputs to assess the potential impacts, or with "bottom-up" or stress-testing adaptation options to explore the sensitivity of hydrology to changed climate conditions by direct adjustments of daily weather generator parameters.
For Richardson-type weather generators, the generation of correlated random numbers is computationally intensive.Instead, an improvement of the Richardson-type model, which involves replacing the preprocessing steps of random number generation with a postprocessing procedure for recorrelating the generated data, has been elaborated in recent years because of its high efficiency and good performance.Currently, the algorithm improvement is only performed for multisite simulation of precipitation without consideration of multivariate correlation (Tarpanelli et al., 2012;Li, 2014).Further extension of this method to multisite multivariate downscaling would promote the application of weather generators for impact assessments of climate change.
Considering the importance of the evaluation of potential changes in hydrological variability, and the difficulty in applying the Richardson-type weather generator to multisite multivariate downscaling, the objective of this study was to extend an efficient multisite precipitation generator, i.e., a two-stage weather generator (TSWG; Li, 2014), to a multisite multivariate GCM downscaling method.The impacts of climate change on streamflow variability in a river basin on China's Loess Plateau were assessed by combining the generated climate change scenarios with a distributed hydrological model.Thus, this study developed a framework for maximizing the application of a parametric weather generator to multisite multivariate simulations, which was shown to be capable of producing information useful for water resource management.
2 Data and methodology

Data description
To project the impact of climate change on hydrology, two datasets are essential: one for climate change scenario development and the other for hydrological simulation.For climate change scenario development, the climate forcing data used here included historical daily precipitation (P ) and maximum and minimum temperatures (T max and T min ) from 18 meteorological stations for the period 1961-2005, and simulated and projected monthly P , T max , and T min for the periods 1961-2005 and 2011-2040, respectively, from five Earth System Models (ESMs).The period 2011-2040 was chosen for the impact study for two reasons: (i) the results from the near-term horizon can be used directly for adaption, and (ii) the uncertainties are minimized, because they increase over time as a reflection of the inherent uncertainties of the ESMs and emission scenarios (IPCC, 2013).
Five ESMs (CanESM2, CSIRO-Mk3.6.0,GFDL-CM3, HadGEM2-ES, and MPI-ESM-LR) were run under a historical emission scenario and four representative concentration pathways (RCPs 2.6, 4.5, 6.0, and 8.5) based on the Fifth Assessment Report of the Intergovernmental Panel on Climate Change.The selection of the five ESMs was based on recommendations following an evaluation of climate model skill (Guo et al., 2013;Chen and Frauenfeld, 2014), and these ESMs are considered to represent the latest accomplishments in climate modeling science and technology.The selected ESMs can provide data for almost all emission scenarios, except CanESM2 and MPI-ESM-LR that have no data for RCP6.0 (Table 1).The four RCPs are named for the radiative forcing values for the year 2100 (i.e., 2.6, 4.5, 6.0, and 8.5 W m −2 , respectively).These RCPs cover most possible future greenhouse emission scenarios, and the ESMs associated with these RCPs have projected significant temperature rises (IPCC, 2013).
The dataset for the hydrological simulation was from the Jing River catchment on China's Loess Plateau (Fig. 1).As the SWAT (Soil and Water Assessment Tool; Arnold et al., 1998) was used to simulate the hydrological cycle, data and/or maps related to climate, soil, vegetation, and hydrology were essential.These datasets were collected from the Data Sharing Infrastructure of the Loess Plateau, and they comprised daily weather data from 18 stations, maps of soil types and properties, a land use map for 1986, and monthly streamflow at the catchment outlet (the Zhangjiashan station).
The Jing River catchment was selected as the study area because it is a typical catchment with high intra-and interannual variability of climate and runoff.This variability has been threatening the management of water resources and soil erosion.The catchment has an area of 45 421 km 2 and it is located within a transition zone between subhumid and semiarid climates.Although the area-averaged annual pre- cipitation was only 542.1 mm, 55 % of the precipitation fell during the flood season between July and September .Several extreme rainfall events during 1961-2000 generated severe soil erosion with an estimated soil loss of 5015 t km −2 year −1 .In addition, because of its dry climate and a runoff ratio of 7 %, the catchment is often subject to severe water shortages.Obviously, the water-related problems within the Jing River catchment are correlated strongly with climatic and hydrological variability.Therefore, the potential impacts of climate change on the hydrological variability of this catchment should be evaluated further.

Multisite and multivariate downscaling
The multisite multivariate ESM downscaling was performed via three steps.In the first step, the ESM outputs were downscaled spatially from a grid scale to a station scale.The second step disaggregated the monthly data further to daily weather series, and the third step reconstructed the multisite and the multivariate correlations.A schematic of the methodology is presented in Fig. 2, and a brief introduction of the methods is given below.
The first and second steps performed spatial and temporal ESM downscaling for single sites.A popular technique for accomplishing this is to use the transfer function method for the spatial downscaling, and then to employ a weather generator for the temporal downscaling.Here, we used the parametric quantile mapping method and a Richardson-type weather generator, respectively, for these two steps (Zhang and Liu, 2005;Li et al., 2011).For the first step, linear and nonlinear transfer functions were fitted to the rank-ordered monthly observations and the ESM data for each calendar month for the period 1961-2005, and then these were applied to the period 2011-2040 to calculate the monthly mean and variance.The nonlinear function was used to transform the ESM monthly precipitation values that were within the range in which the nonlinear function was fitted, while the linear function was used for those values that were outside of the range.Temporal downscaling was then implemented by adjusting the precipitation-and temperature-related parameters of a single-site weather generator (SSWG) calculated from the baseline period.In the SSWG, both the occurrence and the amount of precipitation were simulated using a firstorder two-state Markov chain and a skewed normal distribution based on our previous evaluation (Li et al., 2014), while the temperature was generated using a normal distribution.Thus, the precipitation-and temperature-related parameters for adjustment included the transitional probabilities of a wet day following a wet day (P w/w ) and of a wet day following a dry day (P w/d ), together with the mean and variance of daily precipitation of wet days and the mean and variance of the maximum and minimum temperatures.These were adjusted according to relationships developed based on observations, the detailed procedure of which can be found in Zhang and Liu (2005) and in Li et al. (2011).The adjusted parameters were used to drive the SSWG to obtain climate change scenarios for a 100-year period.To test the performance of our method for the reconstructions of multisite and of multivariate correlations, the temperature generation was not taken as dependent on the dry/wet status in the weather generators, and no intersite or intervariable correlations were taken into account during the single-site downscaling.
For the third step, a method for pairing independent variables to induce the desired rank correlations (Iman and Conover, 1982), which was used successfully in our previous study for multisite precipitation simulation using a TSWG (Li, 2014), was introduced to obtain the multisite and the multivariate correlations.The theoretical basis is described as follows.To assign a desired correlation matrix [C] to a random row vector [X], two steps should be performed.First, with the desired correlation matrix [C].Specifically, for this study, the daily time series of P , T max , and T min , which were generated by the SSWG for all stations and for each month, were placed in one matrix where the rows represented days and the columns represented stations and variables.Then, the ranks of each column were converted to a standard normal distribution by calculating the van der Waerden scores, which can be calculated by −1 {i/(n + 1)}, where −1 is the inverse function of the standard normal distribution, and i represents the rank of each column.The score matrix was multiplied by the decomposed correlation matrix [R ] to obtain a new matrix with the target correlation coefficients.The ranks in the new matrix were used to shuffle the raw matrix, and during this shuffle procedure, the non-positive correlation matrices and the tied ranks due to dry days were adjusted.The non-positive correlation matrices were amended using a spectral decomposition procedure (Rebonato and Jäckel, 2000).The tied ranks were resolved by assigning small values of less than the threshold definition of a wet event (0.1 mm in this study) to dry days.As the data rearrangements perturb the occurrence structure, occurrence adjustment should be undertaken according to those from the SSWG.Detailed descriptions of the procedures can be found in Li (2014).

Hydrological simulation
The SWAT, which is a physically based distributed hydrological model for studying the impact of environmental changes on hydrology (Arnold et al., 1998), was employed to evaluate the response of streamflow to climate change.First, the SWAT was calibrated using observed data, and then it was used to simulate hydrological processes for the period 2011-2040.The two main components in the hydrological cycle, i.e., runoff and potential evapotranspiration (ET 0 ), were simulated using the curve number method (USDA- SCS, 1972SCS, ) 1960SCS, 1965SCS, 1970SCS, 1975SCS, 1980SCS, 1985SCS, 1990SCS, 1995SCS, 2000  and the Hargreaves method (Hargreaves et al., 1985), respectively.The Hargreaves method was preferred over the Penman-Monteith and Priestley-Taylor methods because the climate projection in this study considered only changes in temperature.As the Hargreaves-based ET 0 depends mostly on temperature, the changes in the projected temperature can thus be better translated into evaporation losses.
Monthly streamflow for the period 1960-1970, recorded at the Zhangjiashan gauge station, was used to calibrate the SWAT.The period 1960-1970 was chosen because of the minimal human activity and the small changes in climate.Soil conservation measures as well as the other human activities were minor before 1970 and thus, they had minimal impact on rainfall-runoff relationships.The period before 1970 is used by the Yellow River Conservancy Commission as the baseline against which to assess the effects of soil conservation measures.Therefore, 1961-1970 was considered a reasonable period for SWAT calibration and validation.
The period 1960-1964 was used for model calibration and the period 1965-1970 was used for model validation.Automated calibration/validation and uncertainty analyses were undertaken using the Sequential Uncertainty Fitting version 2 in the SWAT Calibration and Uncertainty Programs (Abbaspour et al., 2007).After sensitivity analysis, the parameters most responsible for runoff simulation were identified and they were calibrated using the objective function of the Nash-Sutcliffe efficiency coefficient.In this study, the sensitive parameters comprised the curve number (CN), base flow recession coefficient (ALPHA_BF), soil evaporation coefficient (ESCO), available water capacity (SOL_AWC), and groundwater delay time (GW_DELAY).The Nash-Sutcliffe efficiency coefficients for the two periods were both 0.78, indicating that the SWAT could satisfactorily simulate the streamflow of the studied catchment (Fig. 3).
Although soil and water conservation projects, implemented in the study basin since the 1990s, have affected runoff processes and streamflow amounts considerably, they were not taken into account in our SWAT simulations, which might have introduced simulation errors for the period 1991- 2005.To exclude the impacts of human activities, the natural runoff represented by the SWAT-simulated runoff for the period 1960-2005 was used hereafter as the baseline against which to emphasize climate-induced changes in runoff.
During the hydrological simulation for both current and future periods, the land surface conditions were assumed invariant.According to our analysis, the land use pattern in 2010 was similar to 1986, although there was some variation during 1986-2010 (Li et al., 2016).Furthermore, the current vegetation coverage is approaching the limit sustainable by the available water resources (Feng et al., 2016).Thus, vegetation coverage cannot increase in the future and according to the land use planning of the local government, it will probably remain invariant.Accordingly, the assumption of invariant land use pattern in this study was considered reasonable.

Performance of correlation reconstruction method
To validate the performance of the proposed method, the downscaled parameters related to changes in precipitation and temperature from RCP2.6 in CanESM2 were used to generate multisite multivariate climate change scenarios.As no correlation was considered during the single-site ESM downscaling, the intersite and intervariable correlations fluctuated around zero (Fig. 4a).However, after rearranging the structure of the data matrix, the multisite multivariate correlations were reproduced well (Fig. 4b).In addition, the averages and standard deviations of monthly precipitation and monthly mean temperatures were reproduced well because the shuffle procedure did not change them from the SSWG (Fig. 4c-h).The slightly underestimated standard deviations of monthly precipitation were caused by inherent weakness of the weather generator, i.e., an underestimation of low-frequency variability.The statistics of the monthly mean temperatures were almost the same as the ESM-downscaled parameters.The above results imply that the proposed method could be considered effective in correlation reconstruction, and that it could satisfactorily reproduce the statistics including low-frequency variability.
To undertake an impact assessment of climate change, the projected climate changes should be transferred to a hydrological simulation.The observed precipitation and temperature for the period 1961-2005 were used to generate 100year climates using the SSWG and TSWG, which were used to drive the SWAT.Then, the simulated hydrological statistics were compared with observations to ensure the correlation reconstruction method did not introduce errors (Fig. 5).Obviously, the three climate series produced similar monthly mean streamflows.The SSWG underestimated the variances and maxima of the monthly streamflow but overestimated its minima, whereas the TSWG produced similar variances and extremes of the monthly streamflow except for a few months.The above results suggested that the developed multisite multivariate climate change scenarios based on the TSWG could be used effectively for simulating hydrological variability and extremes.

Projected climate changes
Compared with that during 1961-2005, the climate projected for 2011-2040 appeared drier and warmer under most scenarios (Table 2), and the trend was more significant under higher RCPs.Averaged over all scenarios, the annual mean precipitation decreased by −1.3 ± 4.4 %, while T max and T min increased by 1.6 ± 0. change (Fig. 6).Precipitation decreased significantly from August to October but it increased from November through to March and into May (Fig. 6a), and temperature increased significantly across all seasons (Fig. 6b and c).Therefore, during 2011-2040, a drier climate would be expected during the flood season, whereas a wetter climate might exist through winter and spring.
The variances in precipitation and temperature during 2011-2040 relative to 1961-2005 tended to increase under most scenarios (Table 3).Averaged over all scenarios, the variance ratios of P , T max , and T min were 1.15 ± 0.13, 1.15 ± 0.14, and 1.04 ± 0.10, respectively.The significance test (p = 0.05) further confirmed the variance increase for P and T max , which suggested that the future climate would be more variable than the present climate.
The variance of monthly precipitation tended to increase under most scenarios and for most months (Fig. 6d); however, the upward trends in variance were only significant for 6 months (i.e., the 1st, 3rd, 5th, 7th, 11th, and 12th months).For temperature, monthly variances increased over the first half of the year and decreased during the second half of the year (Fig. 6e and f), and the significance test further showed that T max variances increased significantly in the first half of the year, except for March and May, while T min variances increased significantly from March to May.Overall, the increase in climate variability was significant during the first half of the year.

Projected changes in hydrological variability
The streamflow in the Jing River was simulated with the SWAT using the 18 climate change scenarios as external forcing (Table 4).Averaged over all scenarios, the annual mean streamflow decreased insignificantly by 1.0 ± 15.0 %.The seasonal patterns of streamflow changes were similar under all RCPs (Fig. 7a).Monthly mean streamflow decreased from September through to November and it increased during the other months; the greatest increase occurred during winter and spring.The t test further showed that both the upward trend from November through to June and the downward trend in September and October were significant, while no significant trend was detected for July and August.
The variances of annual streamflow during 2011-2040 relative to 1961-2005 increased under most scenarios (Table 4).Averaged over all scenarios, the variance ratios of annual streamflow were 1.25 ± 0.55 (p = 0.03), which implies that the interannual variability of future streamflow would become more significant.The variances of monthly streamflow had similar seasonal patterns under all RCPs (Fig. 7b); they increased significantly from November through to August except for January and June, but they decreased in October (p = 0.05), which implies that intra-annual variability would also become greater during 2011-2040.
The maximum/minimum monthly streamflow increased/decreased significantly by 33 ± 22 % and −33 ± 11 %, respectively (Table 5).The monthly maxima increased for most months except October (Fig. 7c, and October was excluded by the t test).The monthly minima decreased for most months except for an increase in January Hydrol.Earth Syst.Sci., 21, 5531-5546, 2017 www.hydrol-earth-syst-sci.net/21/5531/2017/  and February (Fig. 7d), and the downward trend in April and from June through to November and the upward trend in January, February, and May were confirmed by the t test.
The combined effects of the upward trend in the maxima and the downward trend in the minima led to the increase in variability.

Discussion
4.1 Why was the proposed GCM downscaling method used?
This study extended the TSWG model for multisite precipitation simulation (Li, 2014) to a multisite multivariate GCM downscaling method to generate climate change scenarios, which help better project the potential changes in hydrolog- Table 5. Sensitivity of streamflow to changes in means and variances of precipitation and temperature.Pm and Pv -changes in mean and variance of precipitation; Rm, Rv, Rx, Rn -changes in mean (%), variance (%), maximum (%), and minimum (%) monthly streamflow; TXm, TXv, TNm, TNvchanges in mean and variance of Tmax and T min ; k, p -slope and significance level of linear regression.The data in bold font indicate that the relationships are significant (p < 0.05).
ical variability with distributed hydrological models.Compared with the TSWG, the extension first used a larger data matrix to include daily maximum and minimum temperatures, and it then adjusted the input correlation matrix with different pairs of variables.With these adjustments, the proposed method realistically reproduced the climate statistics, and intersite and intervariable correlations (Fig. 4).Furthermore, the computing efficiency of the model was not affected either with or without additional variables.
As stated previously, many downscaling methods require considerable effort regarding parameter estimation and verification.Even though the application of such methods might not be limited by computing power, approaches that have a simple structure and that are easy to implement remain highly desired (Clark et al., 2004;Bárdossy and Pegram, 2012;Li, 2014;Srivastav and Simonovic, 2015;Scheuerer et al., 2017).The proposed method first generates daily series weather for the future period at single sites.It then reconstructs the intersite and intervariable correlations using a distribution-free shuffle procedure.Accordingly, the efficiency of the method lies in two aspects: the easily controlled accuracy because of the availability of multiple algorithms for the generation of climate change scenarios, and the efficient method for correlation reconstruction.
The first step in this study generated climate change scenarios using a Richardson-type weather generator, and the accuracies of the climate statistics were controlled satisfactorily (Fig. 4).However, unlike other methods restricted to specific algorithms, the proposed method is open to any algorithm used for single-site climate scenario generation, e.g., the empirical method of LARS-WG (Semenov and Barrow, 1997) and the circulation-based weather generator of SDSM (Wilby et al., 2002).Furthermore, as the second step does not perturb the statistics of climate from the first step, the advantages of the methods for single-site climate change scenario generation can be adopted.The second step simultaneously reconstructs the multisite and multivariate correlations with one simple shuffle, which is a major advantage of the postprocessing method.This method constitutes a good option for methods that consider only the intersite correlation of one variable, e.g., precipitation or temperature (Kottegoda et al., 2003;Harpham and Wilby, 2005;Kleiber et al., 2012), or that reconstruct the correlations for each variable and for each process, e.g., the preprocessing multisite weather generator (Wilks, 1998(Wilks, , 1999b)).
The method for multisite multivariate correlation reconstruction slightly perturbs the precipitation occurrence (Li, 2014) ods (Clark et al., 2004;Bárdossy and Pegram, 2012).Nevertheless, the suitability of the method for sensitivity analyses or impact assessments in relation to streamflow modeling on the monthly scale in a large catchment (area: 45 421 km 2 ; 18 weather stations) has been demonstrated.To extend the validation of its applicability on different spatial and temporal scales, the method should be applied to daily streamflow simulations and to other watersheds with different catchment areas and climates.

How does climate change influence streamflow variability?
The potential changes in hydrological extremes and regimes have been discussed widely (Bürger et al., 2011;Hirabayashi et al., 2013;Arnell and Gosling, 2014;Dankers et al., 2014;Prudhomme et al., 2014;Schewe et al., 2014); however, the links between climate change and hydrological variability have not been explored fully.Most previous studies have compared the trends or mean changes of climatic variables with those of hydrological variables only qualitatively, providing only a basic understanding of how climate change might influence hydrological variability (Bawden et al., 2014).Nevertheless, the effects of climate change on hydrological extremes are not necessarily a simple function of change in precipitation (Hirabayashi et al., 2013;Arnell and Gosling, 2014;Schewe et al., 2014); instead, multiple climatic and hydrological variables should be examined.
The multiple combinations of ESMs and RCPs for each month used in this study presented the potential for discussing the links between climate change and hydrological variability.Specifically, this study presented the means and variances of monthly climate and streamflow, which could be used to interpret how the changes in mean and extreme climate might influence streamflow variability.Univariate linear regression analysis was used to obtain the sensitivity coefficients between climate changes and streamflow responses.Specifically, based on the data of 18 scenarios accumulated over 12 months, the changes in the means and variances of monthly P , T max , and T min were used to develop relationships with the changes in the means, variances, or extremes of monthly streamflow.The slope of the linear equation was used as the sensitivity coefficient.
According to the sensitivity analysis of precipitation changes during 2011-2040, changes in either the means or in the variances had significant positive correlation with streamflow changes (p < 0.01; Table 5).Precipitation changes had smaller impacts on the mean streamflow than on the streamflow variances; a 1 % increase in the means or variances of precipitation increased mean streamflow by 0.3-0.6 %, while it increased streamflow variances by 1.5-1.9%.Changes in either the means or variances of precipitation had similar impacts on hydrological extremes; a 1 % increase in precipitation increased extreme monthly streamflow by about 1 %.Changes in mean temperatures were correlated negatively with the means and variances of streamflow, while changes in temperature variances were correlated positively with the variances and extremes of streamflow (Table 5).The impacts of changes in mean temperatures were much greater than temperature variances.For example, the streamflow means decreased by 7.9 % (p < 0.01) and increased by 0.1 % (p = 0.09) because of a 1 • C increase in T max means and variances, respectively.However, the changes in temperature means had greater impacts on streamflow variances than on streamflow means.For example, a 1 • C increase in T max means decreased streamflow means significantly by 7.9 % (p < 0.01), while it decreased streamflow variance by 24.8 % (p = 0.06).Considering the significance level, the main factors controlling streamflow variations were found to be the changes in P variances and T max means.
The factors identified in this study as driving streamflow changes, especially for streamflow variances and extremes, differed from studies on historical periods.For example, the trends in the historical streamflow of the Yangtze River in China (Chen et al., 2014) and of the Athabasca River basin in Canada (Bawden et al., 2014) were found to be correlated strongly with precipitation but not related to trends in temperature.In contrast, the results of this study were similar to studies based on hydrological projections for future periods, in which changes in precipitation and temperature both contributed to hydrological variability (Leng et al., 2014;Qiao et al., 2014).In particular, the potential changes in global flood exposure or water scarcity have been found to be correlated positively with increases in temperature (Hirabayashi et al., 2013;Arnell and Gosling, 2014;Schewe et al., 2014).The temperature increase in the future is projected to be greater than in the current period (IPCC, 2013), and its magnitude might exceed the threshold of the temperature that influences streamflow.The different driving factors between the current and future periods highlight the reliability of the results obtained from the sensitivity analysis of this study.

What are the uncertainties and limitations of the projected hydrological changes?
In assessing the impacts of climate change on hydrology, uncertainties are usually attributable to the climate models, emission scenarios, downscaling, and hydrological modeling.In most cases, the uncertainties linked with the climate model structure are greater than those associated with the hydrological model or the downscaling method (Wilby and Harris, 2006;Arnell, 2011;Chen et al., 2011;Gosling et al., 2011).For this study, as the climatic statistics and multisite multivariate correlations were reproduced well and the SWAT was calibrated reasonably well, their uncertainties were considered small.However, the uncertainties from the ESMs and RCPs were comparatively large.Averaged over all the ESMs and RCPs, the annual mean precipitation would change by −1.3 ± 4.4 % (Table 2), and the large uncertainty in the climate projection caused an even larger uncertainty in the runoff change.Specifically, the annual mean runoff changed by 1 % ± 15 %, while the monthly maximum and minimum streamflows changed by 33 ± 22 % and −33 ± 11 %, respectively (Table 4).Averaged over all projections for a certain ESM or RCP, the uncertainty related to the ESMs and RCPs was large (Fig. 8).For example, three ESMs (in the horizontal dash-dotted rectangle in Fig. 8) projected similar changes in precipitation but different changes in temperature.Similarly, three ESMs (in the vertical dashed rectangle in Fig. 8) projected similar temperature changes but quite different precipitation changes.The sensitivity of the uncertainties linked to the ESMs was reduced in the hydrological response (Fig. 9), where four out of five ESMs projected similar changes in the means or variances of runoff.Figure 9 shows that the CanESM2 model contributed greatly to the uncertainty envelope of runoff changes.Thus, over-reliance on a single ESM could lead to inappropriate projection of climate changes.By adopting arithmetic ensemble averages, the uncertainty from ESMs in an impact analysis could be reduced to some extent (Leng et al., 2014).Different from ESMs, the sensitivity of the uncertainties related to the RCPs increased in the impact assessment, where the changes in annual mean streamflow covered greater absolute ranges than those of precipitation and temperature.The projected directions of change in streamflow differed even for the different RCPs.
All the simulations in this study were driven by the ESM data.The use of such an offline approach means terrestrial feedback cannot be taken into account.However, over this inland watershed, land-atmosphere interactions are usually very significant (Li et al., 2010), which could affect the streamflow's magnitude and its variability.Therefore, there is need for further evaluation of the effects of such interactions on runoff processes using a coupled high-resolution climate-Hydrol.Earth Syst.Sci., 21, 5531-5546, 2017 www.hydrol-earth-syst-sci.net/21/5531/2017/ hydrology model, which is a research direction planned for the near future.

Conclusions
Extreme hydrological events such as flooding and drought have great societal importance, especially for a densely populated watershed with fragile environments such as the Jing River (China).In this study, near-term projections of streamflow variability at a watershed scale were produced via a downscaling approach, which has been evaluated fully in both previous research and this study, testifying to the reliability of our results.In addition, the uncertainty ranges of the projections of this study were quantified.Thus, the research results obtained provided a solid basis for deriving clear insights into the resource and hazard management of water supplies.
The projections for the Jing River catchment showed that both annual mean precipitation and streamflow would not change in the future, which implies that the gross amount of the water resource might remain similar to the present level.Further analysis showed that precipitation variances had greater impacts than precipitation means on the streamflow, and that high flows tended to be sensitive to precipitation changes, while low flows were determined by temperature changes.The rising variance in streamflow suggests greater occurrence of floods during the summer and severe water shortages during the winter and spring.This implies that the available water resource could possibly decrease and that soil erosion could become more severe.
Changes in the statistical parameters of precipitation and/or temperature (i.e., the means and variances) played different roles in the mean state and variability of streamflow.The changes in precipitation/temperature were correlated positively/negatively with changes in streamflow.Specifically, the increase in streamflow variances was attributed mainly to the increase in precipitation variances and temperature means, but the positive contribution from increased precipitation variances was larger than the negative contributions from increased temperature means.Competing interests.The authors declare that they have no conflict of interest.Special issue statement.This article is part of the special issue "Coupled terrestrial-aquatic approaches to watershed-scale water resource sustainability".It is not associated with a conference.

Figure 1 .
Figure 1.Location of the Jing River catchment with superimposed digital elevation model (DEM).

Figure 3 .
Figure 3. Observed and simulated runoff at the Zhangjiashan station on the Jing River during 1961-2005 (CAL and VAL indicate the periods used for model calibration and validation, respectively).
Figure 4. (a) and (b) Observed (OBS) versus generated multisite and multivariate correlations of daily precipitation (P ), maximum temperature (T max ), and minimum temperature (T min ).(c)-(h) Average and standard deviations (SD) of monthly P , T max , and T min .

Figure 5 .
Figure 5. Statistics of observed (OBS) and simulated (by SSWG-and TSWG-generated climate) monthly streamflow.Observation period was 1961-1990 and the simulated runoff was from the 100-year climate generated from the statistical parameters from the observed climate.

Figure 6 .
Figure 6.Changes in the averages and variances of monthly climate during 2011-2040 relative to 1961-2005.

Figure 7 .
Figure 7. Changes in average and variability of monthly streamflow during 2011-2040 relative to 1961-2005.

Figure 9 .
Figure 9. Projected changes ( ) in precipitation (P ), maximum temperature (T max ), and streamflow (Q) averaged over the ESMs or the RCPs.
Data availability.The observed climate and streamflow data are available in the Data Sharing Infrastructure of Loess Plateau (2004) (http://loess.geodata.cn/),and the projected climate data were collected from the World Climate Research Programme's Working Group on Coupled Modelling (1989) (https://pcmdi.llnl.gov/?cmip5/).

Table 1 .
Data used for future climate change scenario construction (his; historical).Figure 2. Schematic of the multisite multivariate GCM downscaling and the structure of this study.

Table 3 .
Variance ratios of precipitation and temperature during 2011-2040 relative to 1961-2005.Mean-each/all RCP, average changes for all ESMs under one/all RCP; p-each/all RCP, significance of t test for all ESMs under one/all RCP.The data in bold font indicate that changes are significant (p < 0.05).

Table 4 .
Changes in average and variance of annual streamflow and mean monthly extreme streamflow during the 2020s relative to 1961-2005.Mean-each/all RCP, average changes for all ESMs under one/all RCP; p-each/all RCP, significance of t test for all ESMs under one/all RCP.The data in bold font indicate that changes are significant (p < 0.05).
P m &R x P v &R x P m &R n P v &R n TX m &R x TX v &R x TX m &R n TX v &R n TN m &R x TN v &R x TN m &R n TN v &R n , which is a weakness inherent in postprocessing meth-Figure 8. Projected changes ( ) in annual mean precipitation (P ) and maximum temperature (T max ) averaged over ESMs or RCPs.