**Yiheng Xiang 1, Lu Li 2,\*, Jie Chen 1,3,\*, Chong-Yu Xu 4, Jun Xia 1,3, Hua Chen <sup>1</sup> and Jie Liu <sup>1</sup>**


Received: 16 October 2019; Accepted: 11 November 2019; Published: 18 November 2019

**Abstract:** The impacts of climate change on water resources in snow- and glacier-dominated basins are of great importance for water resource management. The Snowmelt Runoff Model (SRM) was developed to simulate and predict daily streamflow for high mountain basins where snowmelt runoff is a major contributor. However, there are many sources of uncertainty when using an SRM for hydrological simulations, such as low-quality input data, imperfect model structure and model parameters, and uncertainty from climate scenarios. Among these, the identification of model parameters is considered to be one of the major sources of uncertainty. This study evaluates the parameter uncertainty for SRM simulation based on different calibration strategies, as well as its impact on future hydrological projections in a data-scarce deglaciating river basin. The generalized likelihood uncertainty estimation (GLUE) method implemented by Monte Carlo sampling was used to estimate the model uncertainty arising from parameters calibrated by means of different strategies. Future snowmelt runoff projections under climate change impacts in the middle of the century and their uncertainty were assessed using average annual hydrographs, annual discharge and flow duration curves as the evaluation criteria. The results show that: (1) the strategy with a division of one or two sub-period(s) in a hydrological year is more appropriate for SRM calibration, and is also more rational for hydrological climate change impact assessment; (2) the multi-year calibration strategy is also more stable; and (3) the future runoff projection contains a large amount of uncertainty, among which parameter uncertainty plays a significant role. The projections also indicate that the onset of snowmelt runoff is likely to shift earlier in the year, and the discharge over the snowmelt season is projected to increase. Overall, this study emphasizes the importance of considering the parameter uncertainty of time-varying hydrological processes in hydrological modelling and climate change impact assessment.

**Keywords:** Snowmelt Runoff Model; parameter uncertainty; data-scarce deglaciating river basin; climate change impacts; generalized likelihood uncertainty estimation

### **1. Introduction**

The intergovernmental on Climate Change (IPCC) stated that the temperature will continue to increase in the 21st century [1]. Potential impacts of global warming on water availability, especially in snow- and glacier-dominated regions are of great concern [2]. Due to the sensitivity and vulnerability of snow and glaciers due to climate change, previous studies have paid considerable attention to how climate change impacts snow accumulation [3–5], as well as to the timing and amount of snowmelt [6–8]. Water resources in snow- and glacier-dominated regions will suffer from the probable risk of declining amounts of snow and size of glaciers. In order to assess the climate change impacts on the water resources in snow- and glacier-dominated watersheds, it is common to use data from climate simulations to drive snow/glacier hydrological models [9–11]. In these models, snow- and glacier-related processes are usually simulated by energy budget methods [12–14] and/or temperature index methods [15,16], as well as/or multi-layer snow models with snow depth data assimilation methods [17–19]. The physically-based energy-budget methods and snow models incorporate energy fluxes occurring within the snowpack and interacting with the surface [20], while the empirical temperature-index method assumes a linear relationship between air temperature and the rate of snowmelt [21]. Since the energy budget method requires more meteorological data as input, it is often impractical to conduct at high-altitude snow and glacier basins with sparse observations. Temperature-index methods are more practical for snowmelt simulation in data-scarce mountain regions with insufficient data, as the air temperature is relatively widely available compared to other components [21].

The Snowmelt Runoff Model (SRM) is one of the most commonly used temperature-index models. It is a conceptual model developed to simulate and predict daily streamflow based on the degree-day method in high mountain basins where snowmelt is a major contribution [16]. It has been used to assess the climate change impacts on river discharge in mountainous environments [11,22–26]. For instance, Elias et al., [11] investigated the climate change impacts on streamflow timing and runoff volume utilizing the SRM in the Upper Rio Grande basin, Colorado and New Mexico, USA, and found that the large decrease in runoff volume in summer would pose challenges to water management in this region. Tahir et al., [25] applied the SRM to assess the climate change impacts on snowmelt runoff in the Upper Indus Basin and showed that the SRM is effective in estimating snowmelt runoff in high mountain regions. Xie et al., [26] proposed and applied a progressive segmented optimization algorithm (PSOA) for calibrating time-variant parameters of the SRM on the Manas River basin of Xinjiang, China. Their study showed that the PSOA can effectively calibrate the time-variant model parameters while avoiding the high increase in computational time caused by a significant increase of parameter dimensionality.

However, there are multiple sources of substantial uncertainties when applying the SRM for snowmelt runoff modeling and projection. Snowmelt modeling generally involves three uncertainty sources related to forcing data (i.e., input data and data for calibration), imperfections in the model structure and model parameters [27]. In climate change impact studies, uncertainties arise from different climate scenarios, as well as from parameter instability due to the possible changes in a basins' physical characteristics [28]. Among these uncertainties, the identification of model parameters is an essential step for any hydrological modeling, as there are parameters that cannot be derived directly from watershed characteristics but rather need to be inferred via a trial-and-error process [29,30].

Parameter uncertainty may be especially amplified in future climate change impact studies [31–33]. Therefore, quantifying uncertainties due to parameter instability, identifiability and non-uniqueness should be undertaken routinely, and robust methodologies should be implemented. Much attention has been given to uncertainty in model parameters in the hydrological literature [28,31,34–36]. Among the methods developed to deal with parameter uncertainty, the generalized likelihood uncertainty estimation (GLUE) method [37] has been widely used in recent decades [38–41]. The GLUE method is easy to implement and allows a flexible definition of the so-called likelihood function. The likelihood function, which is used to separate behavioural and non-behavioural solutions is particularly valuable for assessment of integrated, distributed models that operate with multi-variable, multi-site and multi-response criteria [39]. Moreover, a numerical Monte Carlo sampling technique [42] which is used to sample from the prior parameter distribution is usually coupled with GLUE for quantifying parameter uncertainties in hydrological modeling [40,43]. Conditioned on the available observational

data and associated uncertainty bounds, this approach produces a well-distributed set of parameter distributions, the acceptability of which are similarly good at producing model simulations.

According to previous studies [44–46], the determination of parameters is essential to utilizing the SRM. Among the SRM parameters, the runoff coefficients (*Cs* and *Cr*) expressing the losses as ratios (*Cs*/*Cr* corresponds to the ratio of runoff contributed by snowmelt/rainfall to precipitation) vary in time, and so are affected by different climatic conditions and land cover properties in a hydrological year. They are usually calibrated if the agreement between observed and simulated runoffs is not acceptable [47]. The calibration strategy most commonly used for time varying parameters is to calibrate them in several sub-periods during the whole calibration period, and then assume that these parameters are constant in each sub-period. For example, Senzeba et al., [48] used the SRM to simulate the snowmelt runoff in a seasonally snow-covered eastern Himalayan catchment in 2006, 2007 and 2009. Single values of *CS* and *Cr* were applied for each calibration year and the mean value from the three calibration years was used for a validation year (hereafter referred to as 1-period calibration). Zhang et al., [49] divided the hydrological year into two sub-periods (i.e., May–June and July–April) to calibrate the runoff coefficients of the SRM (hereafter referred to as 2-period calibration). In the study of Fuladipanah & Jorabloo [50], more than two sub-periods were divided in a hydrological year based on the variations of climatic conditions (hereafter referred to as multi-period calibration).

In the last few decades, the SRM has been applied successfully to mountainous watersheds all over the world, and it was recommended by Abudu et al., [46] that the SRM has a potential to forecast streamflow and to evaluate the effects of climate change on runoff in the mountainous watersheds of northwestern China, especially in data-scarce watersheds in high-elevation regions. However, the parameter uncertainty of the SRM restricts its application in this respect, as its parameter stability is relatively limited due to the highly dynamic behavior of snowpack and snowmelt.

Accordingly, this study aims to evaluate the parameter uncertainty for SRM simulation by using different calibration strategies, and to assess its impact on future hydrological projections. More specific objectives are to: (1) estimate the uncertainty of SRM simulation based on parameter uncertainty (i.e., the overfitting problem) and different calibration strategies (i.e., yearly calibration and multi-year calibration); and (2) assess the uncertainty of future SRM discharge projections determined by utilizing different calibrated parameter sets in a data-scarce deglaciating river basin.

In order to have a wide review of calibration strategies, we apply yearly calibration (which obtains optimal parameter sets from year to year) and multi-year calibration (which gets one optimal parameter value set for many years) with different sub-periods. The GLUE method is used to assess the flexibility and sensitivity of the parameters and to model the uncertainty arising from calibrated parameters.

#### **2. Study Area and Data**

#### *2.1. Study Area*

The data-scarce deglaciating Yurungkash watershed was chosen to implement the proposed study (Figure 1). It is one of the tributaries of the snow- and glacier-fed Tarim River Basin, the biggest inland river basin in China. The Yurungkash watershed is located approximately between 80◦–82◦ N latitude and 74◦–75◦ E longitude in the Xinjiang Uygur Autonomous Region and covers a surface area of about 14,600 km2. It originates from the Akshay Glacier located on the northern slope of Kunlun Mountain, and its elevation ranges from about 1280 m above sea level (asl) up to 6780 m asl. The mean altitude of the watershed is about 4680 m asl and nearly 60% of the area lies above 5000 m asl.

**Figure 1.** The location and topography of the deglaciating Yurungkash watershed with the meteorological station (Hotan) and the discharge station (Tunguz Locke).

The studied watershed is typically snow- and glacier-characterized, with about 40% of the watershed covered by permanent snow and ice. There is no meteorological station installed within the watershed; Hotan station is the closest one (Figure 1), about 600 km far away from the outlet of the river basin, located at 37.13◦ N, 79.93◦ E and 1375 m asl. The mean annual temperature of the study region is about –8 ◦C, and it receives a mean total annual precipitation of 260 mm with 65% falling between June and September. The Tunguz Locke station is located at the basin outlet (36.81◦ N, 79.92◦ E) and records a mean annual flow of 80 m3/s. The annual runoff presents a large intra-annual variability with about 90% occurring in the snowmelt season (April–September).

#### *2.2. Data*

The Shuttle Radar Topography Mission (SRTM), comprised of a modified radar system and offering high-resolution DEM data [51], was used to derive the topography of the study area and delineate the elevation zones. The Yurungkash watershed was divided into six elevation zones and the properties of these zone areas are listed in Table 1. The mean elevation of each elevation zone was calculated by using hypsometric curves, which were also used to derive the temperature of different elevation zones by extrapolation of base station temperature. Daily air temperature observed at the Hotan station and daily precipitation from the China Ground Rainfall Daily Value 0.5◦ × 0.5◦ Lattice Dataset (CGRD) [52] for the 2003–2012 period were used for model simulation. Two different datasets were used for temperature and precipitation because precipitation at the Hotan station cannot well represent the whole study watershed, due to the large spatial variability of precipitation in mountain basins. The CGRD was generated by interpolating observed precipitation from more than 2000 meteorological stations in China, utilizing the thin plate spline approach which was developed to interpolate and smooth scattered data in geosciences [53]. The reliability of this dataset has been proven by Zhao et al., [54]. Daily runoff data at the Tunguz Locke station for the 2003–2012 period was used for SRM calibration and validation.


**Table 1.** Characteristics of elevation zones.

The daily snow-covered area (SCA) is another model input, derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) [55] snow cover product. The MODIS 8-day composite snow cover data product (MOD10A2) is generated by compositing observations from the MODIS Snow Cover Daily L3 Global 500 m Grid (MOD10A1) data set. For every eight-day period, the grid is mapped as snow if snow is observed on any single day, and cloud cover is reported only if the grid was cloud-obscured for all eight days. The MOD10A2 can better eliminate the influence of clouds and improves the precision of snow identification over that of the MOD10A1 [56]. The daily SCA derived from the MOD10A2 by using linear interpolating was used for model simulation in this study. The mean monthly watershed-averaged temperature (T), precipitation (P) and snow-covered area (SCA) for the 2003–2012 period is presented in Figure 2. This figure indicates that the T is above 0 ◦C during May–September. The watershed receives much more precipitation in July and August than in other months. The SCA decreases from April and increases after August, and then decreases again after October. The minimum SCA of the watershed is greater than 40%, which is the area of its permanent snow/ice.

**Figure 2.** The watershed-averaged temperature (T), precipitation (P) and snow-covered area (SCA) (Period 2003–2012). The T, P and SCA were derived from the Hotan station, China Ground Rainfall Daily Value 0.5◦ × 0.5◦ Lattice Dataset and MODIS product, respectively.

For runoff projection, General Circulation Model (GCM) simulations from phase 5 of the Coupled Model Intercomparison Project [57] were used to drive the SRM. Although it is desirable to use as many GCMs as possible, the selection of a representative subset is usually necessary due to the high computational costs of model simulations.

The k-means clustering method [58], which has been used in several previous studies [59–61], was applied to select subsets from 26 GCMs under representative concentration pathway 4.5 (RCP 4.5) and RCP 8.5, respectively. The k-means clustering method divided the ensemble of 26 GCMs into a number of clusters to minimize the within-cluster sum of squared error (SSE). The SSE was characterized by the Euclidean distance based on the standard changes in temperature and precipitation

simulated by GCMs between the 2003–2012 and the 2041–2050 periods. A value of four clusters was determined to offer good group separation as well as a manageable number of simulations. Figure 3 presents the selected subsets of GCMs under RCP 4.5 and RCP 8.5. In each cluster, the colour marker which is closest to the centroid is the selected one. Table 2 presents the detailed information of the seven GCMs used in this study.

**Figure 3.** Sub-sets of the four selected General Circulation Model (GCMs) (color-filled markers) under representative concentration pathway (RCP) 4.5 and RCP 8.5 found by the k-means clustering method from a total of 26 GCMs (non-color-filled markers). These are based on the standard changes in temperature and precipitation simulated by GCMs between the 2003–2012 and the 2041–2050 periods.


**Table 2.** Seven General Circulation Model simulations' details of CMIP5.
