*3.1. SRM Modeling*

The SRM is a conceptual, simple degree-day-based model that was designed to simulate and forecast daily streamflow resulting from snowmelt and rainfall in mountain basins. Since it was developed by Martinec et al., [16], the SRM has been successfully applied to many mountain basins of various sizes ranging from 0.76 to 917,444 km2, and elevations ranging from 0 to 8848 m asl [62]. In recent years, the SRM has also been applied to investigate the impacts of climate change on runoff [7,11,25].

In this study, the SRM is applied to simulate and predict daily runoff for the Yurungkash watershed. It requires three daily input variables (temperature, precipitation and snow-covered area) and has seven parameters (runoff coefficients Cs and Cr, the degree-day factor, critical temperature, rainfall contribution area, recession coefficient and the time lag). A detailed description of the SRM model can be found in the Appendix A.

#### *3.2. SRM Calibration Methods*

Among the SRM parameters, the runoff coefficients (*Cs* and *Cr*) vary with time, and so they are the parameters that are most sensitive to different climatic conditions and land cover properties over a year-long time frame [47]. Therefore, they are usually the parameters chosen for calibration in SRM modelling. In order to have a broad perspective for calibrating runoff coefficients, both the yearly and multi-year calibration strategies with 1, 2, 4, 5, 6, and 7 sub-periods were used for comparison. X-period calibration means that a hydrological year is divided into × sub-periods based on similar characteristics in climatic data for runoff calibration. In total, our approach used 12 calibration strategies. Additional details, including the beginning and ending dates of each sub-period of the different calibration strategies, are presented in Figure 4. For the 2-period calibration, the ending date of the last period is 31st may in the following year. For the other calibration strategies, the ending date of the last period is 31st March in the following year. Apart from *Cs* and *Cr*, other model parameters are presented in Table 3, which were determined as described in the appendix. Calibration was carried out in the odd-numbered years and validation was performed in the even-numbered years for the 2003–2012 period.

**Figure 4.** The beginning and ending dates of each period in a hydrological year for different calibration strategies. Each color in each calibration strategy represents a sub-period. For example, there are 2 sub-periods in the 2-period calibration, the 1st period is June 1 August 31 and the 2 period is September 1 May 31. *a* represents June 20, *b* represents July 15 and *c* represents August 5.

**Table 3.** Parameters determined by hydrological judgment as described in the appendix.


The Nash-Sutcliffe coefficient (NSE) [63] and the relative volume error (VE), two classical efficiency criteria in hydrological modelling, were used to evaluate the performance of the SRM. The NSE and the VE are computed as follows:

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} \left(Q\_i - Q\_i'\right)^2}{\sum\_{i=1}^{n} \left(Q\_i - \overline{Q}\right)^2} \tag{1}$$

$$\text{VE}[\%] = \frac{V\_R - V\_R'}{V\_R} \cdot 100 \tag{2}$$

where *Qi* is the measured daily discharge, *Q*- *<sup>i</sup>* is the simulated daily discharge, *Q* is the average measured discharge of the given time period or snowmelt season, *VR* is the measured yearly or seasonal runoff volume, *V*- *<sup>R</sup>* is the simulated yearly or seasonal runoff volume, and *n* is the number of daily discharge values.

#### *3.3. Bias Correction of GCM Outputs*

\_ENREF\_12GCM outputs are too coarse and biased to directly drive an SRM in climate change impact studies; therefore, a downscaling or bias correction process is needed [64]. In this study, the watershed-averaged value of the GCM–simulated precipitation was first calculated by the Thiessen Polygon method [65]. The GCM–simulated temperature was interpolated to the Hotan station using the Inverse Distance Weighted method [66] in order to be extrapolated later to the higher elevation zones. The precipitation and temperature were then corrected, respectively based on the watershed-averaged precipitation from CGRD and temperature at the Hotan station, by a distribution-based bias correction method, called the Daily Bias Correction (DBC) method in Chen et al., [67]. The DBC is a hybrid method combining the Local intensity scaling (LOCI) method [68] to correct the precipitation occurrence and the Daily translation (DT) method [69] to correct the frequency distributions of precipitation amounts and temperatures. Here are the two steps of this DBC method:


#### *3.4. Future Snow Covered Area Projection*

A changed climate would cause a change in the SCA in mountain snow basins. The method for future SRM simulation developed by Rango & Martinec [70] was used to evaluate the impact of future temperature and precipitation on SCA in this study. This method estimates the time shift of the SCA in the present climate over the snowmelt period to produce the SCA in a future climate. The characteristics of the watershed (i.e., the number of elevation zones, the area and the mean elevation of each zone) and the SCA parameters (i.e., the degree-day factor, the temperature lapse rate and the critical temperature) determined in the present climate, as well as the bias corrected temperature and precipitation in the future climate are required for future SCA projection. More detailed procedures for calculating SCA can be found in Rango and Martinec [70].

#### *3.5. GLUE Method for Uncertainty Estimates*

The GLUE method [37,71] coupled with Monte Carlo sampling [42] was used in this study to investigate the parameter uncertainty of the SRM. In this method, a large number of model runs were made by using many different parameter sets randomly selected from a priori probability distribution. The acceptability of each run was evaluated against observed values and, the posterior parameter distribution was extracted based on a certain subjective threshold to maintain a good population of behavioral solutions. The parameter set with the corresponding model run considered to be non-behavioral was removed from further analysis. In this study, the GLUE scheme associated with the SRM includes the following steps:


The sample sizes of the posterior parameter distributions were all over 1000. These posterior parameter sets were then used in validation and projection. It should be mentioned that for the yearly calibration, the posterior parameter sets of each year were used individually in the validation and projection.

The following three indices were used to evaluate the resolution, reliability and efficiency of the uncertainty interval estimates:

1. The Average Relative Interval Length (ARIL) [34] measures the resolution of the predictive distributions, which is defined over the entire simulation time period as:

$$ARIL = \frac{1}{n} \sum \frac{Q\_{Upper,t} - Q\_{Lower,t}}{Q\_{obs,t}} \tag{3}$$

where *n* is the number of days in the observed record, *QUpper*,*<sup>t</sup>* and *QLower*,*<sup>t</sup>* are the upper and lower simulated discharges of the 95% confidence interval, respectively, and *Qobs*,*<sup>t</sup>* is the observed discharge.

The percentage of observations bracketed by the Confidence Interval (PCI) [72] measures the amount of the observed data within a range of simulated intervals, defined as:

$$PCI = \frac{Q\_{\text{obs},in}}{n} \tag{4}$$

where *Qobs*,*in* is the number of observed discharges that are contained within the 95% confidence interval. A smaller ARIL means a narrower uncertainty interval and a larger PCI means that the interval has a greater reliability; and

2. The percentage of observations bracketed by the Unit Confidence Interval (PUCI) [35] based on ARIL and PCI, which is defined as:

$$\text{PUCTI} = \frac{\left(1 - Abs(PCI - 0.95)\right)}{ARIL} \tag{5}$$

Larger PUCI values represent smaller uncertainty of 95% confidence interval of discharge in general.

For the climate change impacts on runoff projection for the 2041–2050 period, the sub-periods' division in a hydrological year may be different due to the changes in climatic conditions and land

cover properties, which may bring further uncertainty to hydrological projections. This study assumed there is no change of sub-period division in the future climate. To better represent the changed climate, the value of *Cs*, which reflects the decline of the snow coverage and the stage of vegetation growth, was shifted earlier by 30 days in the climate run, as recommended by previous studies [11,16,73]. The climate change impacts on hydrology were assessed at the 95% confidence intervals of average annual hydrographs, annual discharge, and flow duration curves.

### **4. Results**

#### *4.1. Parameter Uncertainty*

The SRM was calibrated based on the yearly and multi-year calibration strategies, with different sub-periods (i.e., 1-, 2-, 4-, 5-, 6- and 7-period strategies) in the odd-numbered years, and validation was performed in the even-numbered years during the 2003–2012 period. Calibration performances with NSE > = 0.55 and −10% < VE < 10%, as well as the corresponding validation performances are shown in Figure 5 for different calibration strategies. This figure shows that: (1) for both yearly and multi-year calibration strategies, the maximum NSE values at the calibration period increase from 1– to 5-period approaches, while there is no obvious improvement for the 6- and 7-period strategy; (2) the maximum NSE values are larger at the calibration period for yearly calibration, but smaller at the validation period than those for multi-year calibration; and (3) there are much more samples with NSE values below 0.55 and absolute values of VE above 10% when using more than two sub-periods at the validation period, especially for yearly calibration.

**Figure 5.** Scatter plots of the Nash-Sutcliffe coefficient (NSE) (upper) and the relative volume error VE (lower) of calibration versus validation for different calibration strategies. The left subplots present the yearly calibration and the right subplots present the multi-Year calibration.

#### *4.2. Confidence Interval of Discharge*

Three uncertainty indices for evaluating 95% confidence intervals of discharge derived from different calibration strategies are given in Table 4. The following features can be observed: (1) both ARIL and PCI increase with the increase of sub-periods for both yearly and multi-year calibration strategies, making it difficult to identify which calibration strategy is better; (2) considering the PUCI, the 1– and 2-period calibrations have a larger value of PUCI than the other sub-periods, which indicates their results are more reliable and less uncertain; and (3) the PUCI values of the multi-year calibrations are generally larger than those of the yearly calibrations, which means the former approach is better.

**Table 4.** The Average Relative Interval Length (ARIL), the percentage of observations bracketed by the Confidence Interval (PCI) and the percentage of observations bracketed by the Unit Confidence Interval (PUCI) values of the different calibration strategies for the 95% confidence interval of discharge by different calibrations.


The 95% confidence intervals of the simulated daily discharge are presented in Figure 6. The results from the 5-, 6- and 7-period calibrations are not shown in this figure because they are very similar to those from the 4-period calibration. Figure 6 shows that the confidence intervals of discharge derived by yearly calibration are slightly larger than those derived by multi-year calibration. The confidence intervals also become larger with the increase in sub-periods, which is consistent with the higher PCI values in Table 4. Furthermore, discharges can be observed lying outside the confidence intervals, which indicates that there are some other uncertainty sources related to forcing data and imperfections in the model structure that impact SRM discharge simulation in addition to the parameter uncertainty.

**Figure 6.** The 95% confidence intervals of daily discharge simulated by posterior parameter sets (color envelope), compared with observed runoff (black lines).

#### *4.3. Projected Future Temperature, Precipitation, and SCA*

By the degree-day method, daily temperature in a future climate modifies the SCA in the present climate to the approximate future SCA in the snowmelt season (April–September). Figure 7 depicts the daily watershed-averaged temperature, precipitation and SCA in snowmelt season from four selected GCMs under RCP 4.5 and RCP 8.5 for the 2041–2050 period (hereafter 'seasonal' means the average value from the snowmelt season). Four GCMs suggest an increase in temperature over the catchment between 1.2 to 2.6 ◦C under RCP4.5, and an increase of between 1.2 to 3.3 ◦C under RCP8.5. The precipitation over the catchment is also projected to increase by 12.0% to 39.8% under RCP4.5, and by 12.2% to 61.8% under RCP8.5. The SCA is projected to decrease in response to the high temperatures in late June. The figure shows that the predicted SCA decreases rapidly from *a* (June 21st, 58.0%) to *b* (July 31st, 40.1%) under RCP4.5 and from *c* (June 11st, 60.3%) to *d* (July 25th, 42.8%) under RCP8.5. The relative decrease in seasonal SCA ranges from 5.7% to 17.7% under RCP4.5, and from 1.2% to 14.8% under RCP8.5.

**Figure 7.** Projected mean monthly temperature (T), precipitation (P) and snow coverage area (SCA) (Period 2041–2050) compared with observations (Period 2003–2012). The pink shading represents the four chosen GCMs, whose mean value is represented by the red line. The same representation is used for the precipitation, indicated in blue. *a* represents (June 21st, 58.0%), *b* represents (July 31st, 40.1%), *c* represents (June 11st, 60.3%), *d* represents (July 25th, 42.8%).

#### *4.4. Uncertainty in the Discharge Projections*

The daily discharges are simulated by the SRM using bias corrected temperature and precipitation from the chosen GCMs under RCP4.5 and RCP8.5, as well as the projected SCA, utilizing the posterior parameter sets derived by the yearly and multi-year calibration strategies with 1-, 2- and 4-period segments. Figure 8 presents the 95% confidence interval of the mean daily discharge for the snowmelt season of the 2041–2050 period. The observed mean daily discharge of the 2003–2012 period is also presented for comparison. The confidence intervals represent uncertainties from four different GCMs as well as from the calibrated parameters. The following results can be observed: (1) the median values of the future discharge simulated by the yearly and multi-year calibration strategies are higher than those of the historical discharge, with the onset of snowmelt runoff shifted earlier. (2) The parameter uncertainty of the yearly calibration is larger than that of the multi-year calibration, as indicated by the wider pink intervals. (3) The discharge confidence interval is wider when there are more sub-periods, in particular for the 4-period strategy, which is much wider than the others.

**Figure 8.** The 95% uncertainty intervals for the GCM–srm predicted daily discharge (Period 2041–2050). The bottom and top of the colored intervals represent the 2.5th and 97.5th percentile values, respectively. In each subplot, the pink and blue shadings represent the yearly and the multi-year calibration strategies, respectively. Median values of the pink and blue envelopes are shown as the red and blue curves, respectively.

Figure 9 shows the 95% confidence intervals of the mean discharge in the snowmelt season from four GCMs under RCP4.5 and RCP8.5 for the 2041–2050 period. The observed mean discharge during the snowmelt season for the 2003–2012 period is also shown. Generally, the projections suggest increases in the mean discharge for the future period, especially under RCP8.5. The confidence intervals of mean discharge from the 1- and 2-period calibrations are smaller than those from the 4-period calibration. In order to compare different sources for uncertainty, these discharge confidence intervals were grouped into five sources (RCPs, GCMs, multi-year/yearly calibration strategies, division of sub-periods, and calibrated parameters). For example, to investigate the uncertainty related to the RCP, discharges were grouped by two RCPs. Each group includes discharges from four GCM, and yearly and multi-year calibration strategies with different sub-periods. The results show that uncertainty ranges related to the choice of RCP, GCM, multi-year/yearly calibration, sub-period and calibrated parameter are [283.11, 308.25], [260.92, 347.93], [267.30, 318.49], [201.86, 229.22] and [165.18, 307.91] with the unit of m3/s, respectively. In terms of the uncertainty range, the largest uncertainty is found to be related to the SRM parameters.

**Figure 9.** The 95% confidence intervals for the GCM–SRM predicted average discharge of the snowmelt season (Period 2041–2050). The bottom and top of the colored bars represent the 2.5th and 97.5th percentile values, respectively. The black bar represents the average snowmelt season discharge during 2003–2012.

The 95% confidence intervals of the flow duration curve (FDC) for four GCMs under RCP8.5 for the 2041–2050 period are shown in Figure 10. Only the high-representative concentration pathway of RCP8.5 is presented to illustrate the uncertainties of the FDC and discuss the differences between GCMs and calibration strategies. Figure 10 reveals that the low, mid, and high flows are all predicted to increase in the Yurungkash watershed for the 2041–2050 period. The highest flow, which occurs less than 10% of the time, is predicted to increase or decrease for different GCMs. This suggests that the number of days with flows exceedingly above 10% will increase in the future, but there will be no significant change in the number of days with the highest flows.

**Figure 10.** The 95% confidence intervals for the GCM–SRM predicted flow duration curves under RCP8.5 (Period 2041–2050) compared with the observations (Period 2003–2012) shown as black line. The dotted and solid lines represent the 2.5th and 97.5th percentile values, respectively.

#### **5. Discussion**

It is widely acknowledged that calibration plays an important role in any hydrological modelling due to parameters that are spatially and temporally heterogeneous or that cannot be measured. Hence, the problem of parameter optimization and the associated uncertainty must be addressed [74]. Our results show that the simulation performance of the SRM is sensitive to calibration strategies with different sub-periods for calibrating the time-varying parameters. In terms of the NSE and the VE, the model's calibration performance becomes better with more sub-periods. However, there are more NSE values below 0.55 when the sub-period is more than 4, indicating higher uncertainty. According to previous studies [32,75], increasing the number of calibrated parameters could yield better calibration performance, while also bringing larger uncertainties at the validation period. An overfitting problem may arise from over-parameterization when the amount of information contained in a hydrograph used for calibration is not enough to estimate a larger number of parameters [76]. An over-fitted model results in a smaller error for calibration but a larger error for validation. Hence, the 1– and 2-period strategies appear to be the most rational choice for calibration when considering the balance between simulation performances and the overfitting problem. In addition, the multi-year calibration would be more stable than the yearly calibration for SRM modeling, as the former performs better at the validation period.

In addition, much attention has been paid to parameter uncertainty assessment in hydroclimatic modeling recently [33,77]. In this study, both the yearly and multi-year calibration strategies with 1-, 2- and 4-periods were used for hydrological climate change impact assessment. Although the runoff projection contains large uncertainty, it shows that the onset of snowmelt runoff shifts earlier and the mean discharge of the snowmelt season increases. This result is similar to that of Wang et al., [78], in which the response of snowmelt runoff to climate change was investigated in a basin located in northwest China.

Particularly, it is noticeable that the predicted hydrographs derived from the 4-period calibration do not shift earlier for peak river runoff. This is not consistent with the common view that the effects of a warming climate will lead to a peak flow shifting from summer and/or autumn to winter and/or early spring [2]. Elias et al., [11] applied a forward shift of 30 days for *Cs* based on a 2-period calibration in the climate change impact assessment and predicted higher streamflow in March and April, which is similar to the results derived from the 1- and 2-period calibrations in our study. The result of the 4-period calibration is different from our expectations because it shifts a large value of *Cs* is shifted to August, which leads to the peak flow occurring in August under climate change. This result indicates that the shift in runoff timing predicted by more sub-periods contains greater uncertainty. In addition, the predicted hydrograph and mean discharge show much wider intervals with the 4-period calibration than with the 1- and 2-period calibrations, which also indicates larger uncertainties. Hence, the 1- and 2-periods are recommended for SRM hydrological climate change impact assessment. This work also suggests that a multi-year calibration strategy would be more stable than a yearly calibration strategy in climate change impact studies.

It is necessary to point out some of the limitations of this study. This study adopted the GLUE method to estimate parameter uncertainty, one drawback of which is that the derived parameter distributions and uncertainty intervals are sensitive to some subjective factors, e.g., the threshold values of the estimation criteria [40,43]. Lower threshold values could result in a wider uncertainty interval of the posterior distribution [34] which may explain why the parameter uncertainty of SRM has the largest impact on the future hydrological projection compared to other uncertainties in this study. This result is similar to that of Tian et al. [79], in which parameter was found to be the major source of uncertainty for simulating future high flows over a southern river basin in China. However, this result is different from that of previous studies [64,80], which showed that hydrological model parameter uncertainty was relatively small compared to other uncertainty sources (RCP and GCM). In previous studies, hydrological model parameters were calibrated by a certain optimization algorithm, which is not recommended for the SRM [47]. However, the use of thresholds in this study cannot guarantee that all parameters are optimized. Nevertheless, the results indicate that the uncertainty associated with a parameter should be highlighted in SRM modeling. It is acknowledged that the parameter distribution and uncertainty intervals derived by GLUE have no clear statistical meaning [81]. Comparison of the results by different threshold values is needed to test the sensitivity of the GLUE, as recommended by Jin et al., [34].

Another limitation of this study is that only different calibration strategies were taken into account in the parameter estimation. It would be more comprehensive to calibrate parameters over different periods, such as wet and dry climate periods [82], as well as by a multi-objective calibration approach [83]. In addition, Seiller et al., [77] showed that the diagnosis of the impacts of climate change on water resources is very much affected by the objective functions used for hydrological model calibration. Investigation of the parameters' nonstationarity in a changed climate could be another avenue for future studies [84].

Also, it is necessary to mention that although temperature-index models are more practical, they are often less accurate than full energy balance models [85]. In addition, the skill of temperature-index models is likely to decrease for future climate change impact studies due to the changes in the energy balance which is not a direct result of temperature (e.g., changes in albedo). The degree day factor has been shown to change over time scales of a decade or even over a single melt season [85–87].

#### **6. Conclusions**

This study evaluated the parameter uncertainty for SRM simulation by different calibration strategies and its impact on future hydrological projections in a data-scarce deglaciating river basin. Different sub-periods were compared and the uncertainty of future discharge derived by the calibrated parameter sets was assessed. The following conclusions can be drawn:


Overall, this research highlights the parameter uncertainty of SRM and its impact on hydrological climate change assessment in a data-scarce deglaciating river basin. The results imply that for hydrological models, like the SRM, time-variant parameters could result in over-parameterization issues and bring uncertainties in the hydrological simulation. Both of these have a significant effect on climate change impact assessment, and therefore should be routinely considered when using these models.

**Author Contributions:** Data curation, Y.X.; Funding acquisition, L.L. and J.C.; Investigation, Y.X.; Methodology, Y.X. and L.L.; Writing – original draft, Y.X.; Writing – review & editing, L.L., J.C., C.-Y.X., J.X., H.C. and J.L.

**Funding:** This research was funded by [the State Key Laboratory of Water Resources and Hydropower Engineering Science funding] grant number [2017SWG02], [the National Natural Science Foundation of China] grant number [51779176, 51525902, 51539009], and [the Thousand Youth Talents Plan from the Organization Department of CCP Central Committee (Wuhan University, China)]. The APC was funded by [the State Key Laboratory of Water Resources and Hydropower Engineering Science funding] grant number [2017SWG02].

**Acknowledgments:** This work was partially supported by the State Key Laboratory of Water Resources and Hydropower Engineering Science funding (No. 2017SWG02), the National Natural Science Foundation of China (Grant No. 51779176, 51525902, 51539009) and the Thousand Youth Talents Plan from the Organization Department of CCP Central Committee (Wuhan University, China). The authors would like to acknowledge the contribution of the World Climate Research Program Working Group on Coupled Modeling, and to thank climate modeling groups for making available their respective climate model outputs. The authors wish to thank the China Meteorological Data Sharing Service System and the Xinjiang Tarim River Management Bureau for providing the dataset for the Yurungkash river basin.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

### *Appendix A.1 SRM Modelling*

The following equation shows the daily discharge estimated by superimposing snowmelt and rainfall on the calculated recession flow:

$$Q\_{n+1} = \left[\mathbf{C}\_{\rm tr} \times a\_{\rm tr}(T\_n + \Delta T\_n)\mathbf{S}\_n + \mathbf{C}\_{\rm tr} \times P\_n\right] \times A \times a \times (1 - K\_{n+1}) + Q\_{\rm tr}K\_{n+1} \tag{A1}$$

where:

*Q* = the average daily discharge (m3/s)

*Cs*/*Cr*= the runoff coefficient expressing the losses as a ration of runoff to precipitation, with *Cs* referring to snowmelt and *Cr* referring to rainfall

*a* = the degree-day factor (cm/ ◦C/day)

*T* = the number of degree-days (◦C day)

Δ*T* = the adjustment by temperature lapse rate when extrapolating the temperature from the station to the mean elevation of the zone (◦C day)

*S* = the ratio of the snow covered area to the total area (%)

*P* = the precipitation contributing to runoff (cm), which is determined by a preselected threshold temperature, *Tc*, to be rainfall or snowfall. The contribution of rainfall is immediate while snowfall will be kept on storage until melting conditions occur.

*A* = the area elevation zone (km2)

<sup>α</sup> <sup>=</sup> 10,000/86,400, the coefficient converting data from a runoff depth (cm×km2/day) to discharge (m3/s)

*Kn*+<sup>1</sup> <sup>=</sup> the recession coefficient, *Kn*+<sup>1</sup> <sup>=</sup> *Qn*+<sup>1</sup> *Qn n* = the day number

It is recommended to subdivide the basin into several elevation zones if the elevation range exceeds 500 m. The simulated discharge of each zone is summed to estimate the runoff of the basin. The study area was divided into 6 elevation zones, characters of which are present in Table 1. Model structure are

shown above, the determination of input variables and parameters are declared below:
