Future Climate Change and Its Impact on Runoff Generation from the Debris-Covered Inylchek Glaciers, Central Tian Shan, Kyrgyzstan

The heavily debris-covered Inylchek glaciers in the central Tian Shan are the largest glacier system in the Tarim catchment. It is assumed that almost 50% of the discharge of Tarim River are provided by glaciers. For this reason, climatic changes, and thus changes in glacier mass balance and glacier discharge are of high impact for the whole region. In this study, a conceptual hydrological model able to incorporate discharge from debris-covered glacier areas is presented. To simulate glacier melt and subsequent runoff in the past (1970/1971–1999/2000) and future (2070/2071–2099/2100), meteorological input data were generated based on ECHAM5/MPI-OM1 global climate model projections. The hydrological model HBV-LMU was calibrated by an automatic calibration algorithm using runoff and snow cover information as objective functions. Manual fine-tuning was performed to avoid unrealistic results for glacier mass balance. The simulations show that annual runoff sums will increase significantly under future climate conditions. A sensitivity analysis revealed that total runoff does not decrease until the glacier area is reduced by 43%. Ice melt is the major runoff source in the recent past, and its contribution will even increase in the coming decades. Seasonal changes reveal a trend towards enhanced melt in spring, but a change from a glacial-nival to a nival-pluvial runoff regime will not be reached until the end of this century.


Introduction
Fresh water and its availability for human use are very important issues in Central Asia. In the Tarim catchment in northwestern China, the mean annual precipitation is just 87 mm [1]. In this region, cultivation of land is only possible using irrigation [2]. However, agricultural land area more At Langtang Himal, Nepal, significantly lower thinning rates at debris-covered glaciers compared to clean ice glaciers were measured [34].
Supraglacial lakes and ice cliffs [39,40] are very common on debris-covered glaciers and need to be considered. Several field studies and modelling attempts have shown that these are areas of increased ice melt [33,[41][42][43][44]. At ice cliffs, high melt rates result from the low albedo of the steep, dust covered slopes and thus are highly dependent on the shortwave radiation and therefore on the ice cliff orientation [45]. At supraglacial lakes, mass losses are mainly caused by lateral water line-melting and calving [46]. At the lake bottom, heat convection is inhibited by abundant debris of particle size between clay and silt with an almost negligible permeability [41]. For this reason, melt at the bottom is strongly dependent on the debris' thermal conductivity and relatively small.
On Northern Inylchek Glacier, melt rates below debris cover were observed by the authors during field investigations in 2008, 2012 and 2013. Using these findings, a model based on HBV-ETH was further developed to also include sub-debris melt [47].
In this study, we present a spatially distributed conceptual model, which is able to cope with supra-glacial debris cover, as well as ice cliffs and supraglacial lakes. This model is now referred to as HBV-LMU.
The following scientific questions will be addressed: (1) How will global climate change affect the regional climate in Central Tian Shan and how large is the uncertainty between the used scenarios? (2) How do debris cover, ice cliffs and supraglacial lakes influence the melt rates? (3) How will melt and runoff at the Inylchek glaciers change in the future considering different climate scenarios?

Study Area
The Inylchek glaciers are located in the Jengish Choqusu (Pobeda)-Khan Tengri massif in the Central Tian Shan (Figure 1), which is mainly built up by marine Paleozoic formations. These are folded schists, calcareous schists, siltstone and sandstone [48]. As a consequence, the clastic material of the moraines around and on the glaciers are weathered products of these metasediments. On such siliceous parent material develop acidic mountain soils that are generally poorly developed and skeletal. So-called leptosols are shallow and have a small water retention capacity, therefore they have a very small hydrological impact. In the Inylchek catchment, soils develop only on flat parts that were not glacier covered during the Little Ice Age.
The climate in the study area is generally continental with an amplitude of monthly air temperatures of around 30 K, but with increasing elevation, precipitation increases and temperature differences between summer and winter decrease [49]. At Koylyu meteorological station (2800 m a.s.l.) outside the catchment (Figure 1), annual mean temperature is −1.9 • C, with a minimum of −19.2 • C in January and a maximum of 10.3 • C in July. Westerlies from the Atlantic are the main source of humidity and the precipitation climate is mainly controlled by the interaction between the Siberian anticyclone during winter (dry conditions) and cyclonic activity from the west during summer (wet conditions). Thus, the Central Tian Shan experiences a distinct summer maximum of precipitation [49]. At Koylyu station, annual precipitation is 313 mm, 54% of which occurs in summer (Jun-Aug).
Located above the tree line, vegetation is limited to alpine tundra on a few lower and flatter locations of the catchment. Pioneer species like cushion plants exist in higher and steeper terrain, while in the snow zone no higher plants exist. Overall, vegetation is very scarce and transpiration can be neglected in the water cycle.
Until the second half of the 19th century, one joined, single Inylchek Glacier covered the Inylchek valley and its upper tributaries. Due to the warming trend following the Little Ice Age, the northern glacier part lost connection to the main part in the south and the glacier was separated into the Northern and the Southern Inylchek Glacier [50]; in the following called Northern and Southern Inylchek ( Figure 1). Southern Inylchek acts as an ice dam towards the valley of Northern Inylchek and the runoff of Northern Inylchek is blocked. This water, together with calving ice from the dam, forms an ice-dammed lake called Lake Merzbacher. The lake has an approximate size of 4 km 2 and a maximum depth of 100 m near the ice dam [51]. Lake Merzbacher empties nearly every year, causing a glacial lake outburst flood (GLOF), generally occurring in July or August.
Water 2018, 10, x FOR PEER REVIEW 4 of 26 Until the second half of the 19th century, one joined, single Inylchek Glacier covered the Inylchek valley and its upper tributaries. Due to the warming trend following the Little Ice Age, the northern glacier part lost connection to the main part in the south and the glacier was separated into the Northern and the Southern Inylchek Glacier [50]; in the following called Northern and Southern Inylchek ( Figure 1). Southern Inylchek acts as an ice dam towards the valley of Northern Inylchek and the runoff of Northern Inylchek is blocked. This water, together with calving ice from the dam, forms an ice-dammed lake called Lake Merzbacher. The lake has an approximate size of 4 km 2 and a maximum depth of 100 m near the ice dam [51]. Lake Merzbacher empties nearly every year, causing a glacial lake outburst flood (GLOF), generally occurring in July or August. Northern Inylchek experienced a surge in 1997 and slightly retreated again afterwards [53]. Today, Northern Inylchek covers an area of 152 km 2 , while Southern Inylchek has a glacier area of 418 km 2 . Their common catchment has a total area of about 800 km 2 [54]. Both glaciers are heavily debris covered. The debris thickness ranges from a few millimeters up to several decimeters at the end of the glacier tongues. As at many debris-covered glaciers, the debris areas are interspersed with glacial-karst structures. The melt water of both glaciers forms Inylchek river, which drains into Sary-Dshaz river, the main tributary of Aksu river.

Hydrological Model
The HBV-LMU model used in this study is a raster-based temperature index model with a grid size of 180 × 180 m. It consists of four storages (Figure 2), is controlled by 23 free parameters (Table  1) and calculates accumulation, melt and runoff on a daily time step. Northern Inylchek experienced a surge in 1997 and slightly retreated again afterwards [53]. Today, Northern Inylchek covers an area of 152 km 2 , while Southern Inylchek has a glacier area of 418 km 2 . Their common catchment has a total area of about 800 km 2 [54]. Both glaciers are heavily debris covered. The debris thickness ranges from a few millimeters up to several decimeters at the end of the glacier tongues. As at many debris-covered glaciers, the debris areas are interspersed with glacial-karst structures. The melt water of both glaciers forms Inylchek river, which drains into Sary-Dshaz river, the main tributary of Aksu river.

Hydrological Model
The HBV-LMU model used in this study is a raster-based temperature index model with a grid size of 180 × 180 m. It consists of four storages (Figure 2), is controlled by 23 free parameters (Table 1) and calculates accumulation, melt and runoff on a daily time step.
The model is based on the HBV3-ETH9 model version [55], which is a classic conceptual runoff model with a glacier routine that has been widely applied, also in Central Asia (e.g., Reference [23]). This general model approach was chosen because of its low data demand. Only precipitation and air temperature on a daily time step serve as input. Furthermore, the model showed a robust performance in earlier studies [20,21]. In order to reproduce sub-debris melt, the semi-distributed HBV-ETH has been set up partially distributed on a raster basis and was extended by several features. These are an automatic calibration procedure, allowing to introduce additional calibration criteria, and a routine to account for snow redistribution [56]. Finally, a routine to quantify melt below supraglacial debris, developed at Inylchek glaciers, has been implemented [47]. In this contribution, the calculation of melt at ice cliffs and supra-glacial lakes has also been included (see Section 2.2.2.) in the new model version (HBV-LMU). Owing to the special situation in the Inylchek catchment, the model also considers the buffering effect of Lake Merzbacher, but it cannot describe the outburst flood (see Section 2.2.4).  Table 1 for abbreviations. The model is based on the HBV3-ETH9 model version [55], which is a classic conceptual runoff model with a glacier routine that has been widely applied, also in Central Asia (e.g., Reference [23]).  Table 1 for abbreviations.  Based on large-scale input fields of accumulation, the spatial variation is calculated according to Huss et al. [19]. Steep slope angles and convex curvatures diminish accumulation while at concave curvatures, accumulation is increased. This is controlled by the parameters CURV, SPREAD, SMIN and SMAX in the model.
Melt of snow and bare ice (M S , M BI ) is calculated if the daily mean temperature (T) exceeds the threshold temperature (T0) based on an extended degree-day approach [57]. There, the degree-day factor is modified to consider shadowing and sun height at each grid cell. To do so, potential clear sky solar radiation (I) at every day and grid cell was calculated and added to the degree-day factor, which is now called melt factor MF, due to the modification. The relative influence is controlled individually for ice and snow areas using the parameters R Snow and R Ice .
Snow melt infiltrates the remaining snow cover. The amount of stored water depends on the actual snow water equivalent and a parameter for water-holding capacity of snow (CWH). This water can refreeze and will be added to the snow cover if the daily mean temperature falls below T0, controlled by a refreezing coefficient (CRFR).

Ablation at Debris-Covered Areas, Ice Cliffs and Supraglacial Lakes
Ablation was measured at 26 ablation stakes from 14 July to 4 August 2012 in order to calculate empirical degree-day factors. Air temperature was recorded at a weather station located at the ice dam of Northern Inylchek close to Lake Merzbacher. Measured temperatures were extrapolated from the weather station located at the ice dam of Northern Inylchek to higher elevations using a lapse rate of 0.008 K m −1 . This lapse rate has been determined by Han et al. [58] and is based on measurements at Koxkar glacier, which is also located in the Central Tian Shan. Ablation measurements were conducted similar to previous measurements at Koxkar Glacier [31].
The empirical relationship between debris cover thickness (DCT) and measured ablation was determined. The results are comparable to other studies [31][32][33]35]. Based on the relationship between DCT and ablation in dependence of the air temperature, the degree-day factors for debris-covered areas (DDF DI ) are derived in two different ways, depending on the DCT. For very thin debris layers, increased melt is approximated by a linear function up to the thickness of maximum melt (in our case 0.0058 m, Equation (2)). This formula also provides the empirical degree-day factor for bare ice (DDF BI ) for the case DCT = 0. For debris layers thicker than 5.8 mm, the degree-day factor is calculated by a power law regression (Equation (3)): DDF DI/BI = 0.0953 × DCT + 0.0063 : DCT ≤ 5.8 mm, DDF DI = DCT −0.3625 × 0.0011 : DCT > 5.8 mm.

of 26
The spatially distributed melt amount M DI (Equation (1)) is then determined by multiplying the potential bare ice melt M BI with the ratio of DDF DI to DDF BI at every grid cell: Melt at ice cliffs was also determined by ablation stakes in 2012 and empirical degree-day factors to model melt at ice cliffs in three different aspects were calculated (DDF [m d −1 • C −1 ]: South: 0.0085, North: 0.0053 and East/West: 0.0071), assuming that the majority of ice cliffs have a characteristic slope. A similar ratio to the bare ice conditions as above was calculated for the ice cliff situation and scaled with the respective bare ice melt.
At supraglacial lakes, lateral melting and calving processes could not be considered in our model, because no parameterization of these processes are available. As melt rates at the bottoms of supraglacial lakes only depend on the debris' thermal conductivity [41], they can be calculated based on the same relation of degree-day factor and DCT, assuming constant conductivity for the debris layer and a similar debris thickness in the lakes as in the debris-covered surroundings. In addition, a constant water temperature of 1 • C is considered [59].
Both ice cliffs and supraglacial lakes are mostly too small to be captured by the grid size of our model. For this reason, the percental areal share of ice cliffs and lakes to each grid cell was used as model input and considered for melt calculation.

Soil Moisture Storage and Runoff Generation
Melt water from snow and ice, as well as liquid precipitation, are summed up and form the hydrological input into the soil moisture routine ( Figure 2). There, they fill the soil moisture storage (SSM). Loss due to evaporation from this storage was calculated based on the filling level of the storage and the potential evapotranspiration (ETpot) and a threshold value LP. ETpot is assumed if the SSM fill level exceeds the threshold LP. Infiltration of the remaining water into the upper storage was governed by SSM, the maximum soil storage capacity (FC) and a coefficient (BETA). From this upper storage, surface runoff (Q0) was calculated as soon as the storage filling level (SUZ) exceeded a threshold value (LUZ). The amount of Q0 is controlled by the parameter k0. Interflow from this storage was controlled by the filling level SUZ and the parameter k1. Constant infiltration into the lower storage was controlled by the parameter CPERC (see Table 1). Linear drainage from this storage was calculated based on its filling level SLZ and the parameter k2. The daily runoff of the catchment finally consists of the sum of the runoff volumes of Q0, Q1 and Q2.

Consideration of Lake Merzbacher
To simulate runoff at the Inylchek glaciers, it is necessary to consider the special catchment situation caused by Lake Merzbacher. First modelling attempts [54] revealed that, even in the phase of lake filling, the blocking ice dam is not a completely closed barrier. A permanent outflow from the lake needs to be considered by the model. In this study, this outflow is calculated based on the individual lake filling volume multiplied by an outflow variable. To determine this variable, individual calibration runs were conducted and the simulated lake filling volumes were compared to outburst volumes determined by Ng and Liu [60], for more details see Mayr et al. [54]. This revealed a best fitting daily drainage of 0.35% of the actual lake volume.

Meteorological Input Data
Meteorological time series for air temperature and precipitation are available from Koylyu meteorological station (42 • 20 N, 79 • 00 E, elevation 2800 m a.s.l.), which is located 47 km west of the glacier terminus [60] (Figure 1). This weather station provided a gapless record of precipitation, air temperature, air humidity, air pressure and wind speed from 1951-1990. Unfortunately, these data were available in time steps of 10 days only. Based on this, daily temperature data were obtained by downscaling of ERA40 reanalysis data [61] using the regional climate model REMO ("Regional Model") [62] Version 2009. Precipitation from the ERA-40 reanalysis has not been taken into account because it is supposed to be of low quality, due to the low resolution of the underlying model in the topographically diverse terrain of Central Asia. Indeed, the comparison between 10-days precipitation time series from ERA-40 and available stations revealed large discrepancies in terms of mean, standard deviation and phase relationship (not shown). Instead, daily precipitation values were generated by a Markov Chain approach and a Gamma-Distribution [63]. The according 10-daily means were compared to the measured means and the daily data corrected accordingly using a Delta-T-approach under consideration of the generated dryness or wetness of the days. The so derived daily data were used for model calibration.
The virtual point precipitation at the weather station is extrapolated by the hydrological model to the catchment by vertical gradients. To account for the representativeness of the weather station, the model implies calibrated correction factors for liquid and solid precipitation. This is standard procedure to create basin precipitation in conceptual runoff models. However, the so derived precipitation is subject to large uncertainties. Especially in glacierized catchments, where runoff can be generated by both precipitation and glacier melt, the danger of error compensation occurs if the model is calibrated by runoff only. In such a case, runoff simulations can show high efficiencies despite incorrect reproductions of precipitation and melt. For this reason, it is very important to use snow-and glacier related criteria for calibration and validation in order to constrain model parameters and thus reduce parameter uncertainty [64,65].
To simulate melt and runoff under past (1970/1971-1999/2000) and future (2070/2071-2099/2100) climate conditions, we generated meteorological input data using two different downscaling approaches: One dynamical approach using a regional climate model, and a statistical-dynamical approach. To keep results comparable, the ECHAM5/MPI-OM global climate model [66,67] has been used as boundary condition for both downscaling approaches. The choice of this GCM has been motivated by the fact that, at the time of producing the meteorological input, it was the climate model with highest spatial and temporal resolution provided via the CERA database and it still is a well evaluated global climate model for the Central Asian target region [16,68].
For the dynamical downscaling approach, we use the regional climate model REMO [62]. For reasons of limited computational resources, we could realize only one long-term climate change projection with REMO using the moderate A1B emission scenario. It represents an intermediate pathway of future climate change within the bunch of more recent RCP scenarios and is still often used for climate change studies. The resulting regional climate model data was bias corrected based on station data, using a linear correction for temperature [69] and a distribution derived transformation using the Gamma-distribution for precipitation [70]. In the following, these data are referred to as "REMO".
For the statistical-dynamical methodology a circulation weather type (CWT) approach following Jones et al. [71] is utilized. This approach is applied to reanalysis (for calibration only) and GCM data in order to classify the daily lower-level atmospheric characteristics using 700 hPa geopotential height (GPH) fields as an input (for details, see Reference [17]). For future decades the three SRES scenarios B1 (optimistic, CO 2 concentration increase from 367 ppm in the year 2000 to 540 ppm in 2100), A1B (moderate, 367 ppm to 703 ppm), and A2 (pessimistic, 367 ppm to 836 ppm) are considered. For each scenario three realizations are available. The statistical moments (e.g., mean and standard deviation) of the temperature and precipitation distribution were calculated for each CWT and season independently. To obtain continuous meteorological time series a random number generator (RNG) was used to generate temperature and precipitation time series that correspond to the statistical moments that were calculated for each CWT and season. Finally, the generated time series were recombined and weighted with the frequency of occurrence of each CWT [17]. To generate time series of precipitation and temperature from global climate projections, only the frequency of occurrence of each CWT was considered. Thus, this approach provides an estimate of the lower bound of possible changes in temperature and precipitation over this region as it does not consider possible changes of rainfall and temperature within the CWTs. These datasets are referred to as "RG" (for random generator) with their SRES scenario (e.g., RGB1, RGA1B and RGA2) in the following.
The downscaled meteorological input for the hydrological model does not arise from a multi-GCM multi-RCM ensemble approach. However, uncertainty of past and future meteorological boundary conditions is imposed by the two different downscaling approaches (dynamical versus statistical-dynamical) and by the statistical assessment in the framework of model data post-processing (see above), replacing the ensemble technique to a certain extent. In addition, scenario uncertainty is introduced by considering three different emission scenarios within the RG methodology.

Calibration Data
Runoff was measured at the gauging station near Inylchek village, 50 km downstream of the glacier snout (Figure 1) in the summer months of 1962-1965 and 1980-1981. The low flow winter period was not recorded and was therefore estimated manually by linear interpolation [47]. For model calibration, both, meteorological and runoff data are necessary, limiting the calibration years to 1963/1964, 1964/1965, and 1980/1981. Snow cover information was extracted from Level 1b data of the Advanced Very High Resolution Radiometer (AVHRR) onboard the NOAA-9 satellite with a spatial resolution of 1100 m. If more than 50% of the model pixels inside one (non-) snow-covered AVHRR pixel were also (non-) snow-covered, this AVHRR pixel was counted as correctly attributed. The accuracy of snow cover classification from AVHRR is discussed in detail in Peters et al. [72]. They found a true positive rate of 87% between classified and directly observed snow in the Aksu basin. We chose seven images between April and September 1986 for our analyses, because 1986 was the only year with AVHRR and meteorological data coverage. The raw AVHRR images have to be pre-processed prior to snow cover extraction. This includes radiometric calibration and calculation of albedo and radiances, as well as the correction of sensor degradation and non-linearities in the measurements. Additionally, the images were geometrically rectified and referenced to UTM. The calibrated and referenced images were classified into the classes "snow", "no snow", and "clouds" using a dichotomous multi-channel classification scheme [73]. The model relies on surface temperature, albedo and NDVI to separate the classes. The thresholds of the subsequent steps were adapted to seasonal characteristics to improve accuracy. Cloud coverage has been reduced by employing a modified version of the Modsnow cloud removal algorithm [74].

Validation Values from Literature
Several glaciological studies have already been conducted at the Inylchek glaciers, their results could partly be used for model validation (Table 4). Aizen et al. [75] determined the mean mass balance at an elevation of 6147 m a.s.l and the equilibrium line altitude (ELA) for the years 1969-1989. Regional mass balances were calculated for the Inylchek glaciers based on surface elevation changes derived from multi-temporal digital terrain models (DTMs) [11] and ice transport rates [54]. For more details about the calculation see Mayr et al. [47].

Spatial Input
The model requires spatial information about glacier area, debris area, debris thickness, as well as the distribution of ice cliffs and supraglacial lakes. These elements may change over time, especially in a changing climate [42,76,77]. For the model calibration period and the past period 1971-2000, different sources of remote sensing data were used to generate the model input "Past Area" while for the future, a model input with possible changes was calculated, called "Best Guess" in the following. Additionally, three more glacier surface conditions were created to conduct a sensitivity analysis. Spatial input data for calibration period and past period ("Past Area") The glacier area and the debris-covered glacier part were manually delineated based on Hexagon KH-9 data from 1974 [11].
DCT was estimated by an empirical approach using thermal satellite imagery and local debris thickness measurements [31]. Measurements are available for 13 locations along a longitudinal profile (elevation from 3200 to 3300 m a.s.l.) in 2013. Land surface temperature (LST) was obtained following the approach by Pi et al. [78], using ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) Level 1B granule images of 12 May 2007. These data have a ground resolution of 30 m, they have been resampled to 180 m to match the model resolution. Mihalcea et al. [79] observed a strong correlation between ASTER-derived LST and debris thickness. Furthermore, in this study, debris thickness measurements in 25 pixels of the ASTER image showed a good correlation with surface temperature, thus LST could be used to estimate debris thickness distribution.
Ice cliffs and supraglacial lakes were mapped manually using a RapidEye image from 2011 [80]. According to the product specifications, RapidEye images have a positional accuracy of less than 10 m. This is more than enough for our purpose, where the absolute position of the features is not that important. Ice cliffs were differentiated into north, south and east/west facing cliffs. For melt calculation, the actual cliff area had to be calculated based on this horizontal data. Sakai et al. [41] found a range of ice cliff angles between 30 • and 60 • at Lirung Glacier. Therefore, a standard slope angle of 45 • was used here [33].
Debris distribution, as well as ice cliff mapping, was performed using more recent remote sensing data compared to the calibration period. A comparison of the two periods revealed only minor changes in debris area and ice cliff distribution. Lakes seem to be slightly more numerous in 2011. This variation has very probably only a small impact on the melt calculation, because the lakes cover only 0.24% of the total glacier area in 2011.
Spatial input data for future period ("Best Guess") Glaciers in the Tian Shan lost both area and mass during the last decades [10,[81][82][83][84]. The major reason for the ice loss can be attributed to increasing air temperatures and changes in precipitation portioning [81,84]. Because increasing temperatures in the future will further enhance ice and snow melt, glacier area and mass are also expected to further decrease. Additionally, intensified melt and reduced ice flux increases the debris-covered area extent [85,86]. For this reason, the available spatial data had to be adapted for the projection period from 2071 to 2100.
To adapt the debris cover, we followed the approach of Jouvet et al. [87] and compared the RapidEye image from 2011 with a Corona KH-4B scene from August 1967 [88]. The KH-4B image was orthorectified using Remote Sensing Software Package Graz Version 7.46 with an RMS-X of 4.67 px and a RMS-Y of 4.02 px. The observed rate of debris front propagation was used to calculate a possible debris cover in 2085, as the middle of the 30 years' future period. The influence of a possible increase of negative mass balances in the future was not considered. For this reason, the debris expansion might be slightly underestimated.
To adjust debris cover thickness, the following assumptions were made: Anderson [85] postulates that medial moraines show a relatively uniform debris thickness of only a few centimeters. The moraine widths increase downstream, but the thickness is relatively constant. Debris thickness starts to increase in areas where the moraines merge. Furthermore, at Inylchek, it is visible that debris at the moraine margins tends to be slightly thinner than in the middle. Based on these assumptions, new debris-covered areas were set to the mean debris thickness of the adjacent former marginal pixels. At these pixels, the debris thickness is supposed to increase. Their new debris thickness was calculated as the mean of their former debris thickness and the adjacent non-boundary pixels.
Volume changes of debris-covered glaciers are not reflected by variations of the glacier fronts, which commonly reveal a very stable behavior [89]. For this reason, the usual methods for glacier area adaption like the accumulation-area-ratio method [90] cannot be applied at debris-covered glaciers.
To estimate the potential future glacier area, we used a similar method as for debris cover evolution. Consequently, images from 1967 (Corona KH-4B) and 2011 (RapidEye) were compared to identify the area decrease for 100 m elevation bands. To take a future atmospheric warming into account, a linear equation was derived based on these measurements and the ELA (4400 m to 4500 m) [75] as level with no area loss. To calculate the future glacier area, the ELA was shifted upwards. Model simulations performed during this study using the original glacier area under future climate conditions revealed that the ELA will increase up to 5300 m. The mean between the old ELA of 4400 m and the new one was used as mean ELA of our calculations. Based on this, it was possible to calculate the future area losses based on the previously mentioned equation. This approach results in a glacier area reduction of 22.5%.
Spatial input data for sensitivity analysis As first modelling results revealed that runoff will still increase with reduced glacier area, we decided to conduct a sensitivity analysis and generated three more glacier inputs. First and second, the areal loss of the above explained "Best Guess" input was increased by 300% and 350% respectively, leading to area reductions of 40.9% and 43.0%. Third, in order to analyze the effect of the debris cover on ice melt Northern and Southern Inylchek, we generated a hypothetical glacier extent without debris cover, ice cliffs and supraglacial lakes based on the glacier area of the Best Guess.

Model Calibration and Validation
To calibrate conceptual hydrological models, automatic calibration algorithms are frequently used [91]. One very promising approach is the multi-algorithm, genetically adaptive method (AMALGAM) of Vrugt and Robinson [92]. This algorithm was already successfully applied for the calibration of hydrological models [47,93].
To obtain valid calibration results, the multi-objective calibration approach is established. Additional to runoff, other measured data, such as seasonal and annual mass balances, are used for calibration. The simulation quality is assessed by the difference between measured and simulated data, the so called "objective functions". The multi-objective approach usually leads to a reduced quality of runoff simulation, but better mass balance simulation [19,56,65,94,95]. If mass balance data are not available, snow cover from remote sensing imagery (e.g., MODIS, AVHRR) can be used as additional calibration dataset [96][97][98].

Objective Functions
The three years with all necessary data for model calibration were divided into differential split samples of two calibration years and one validation year [99,100]. In the following, these are called Sample 1, Sample 2 and Sample 3.
Runoff can be evaluated using the Nash-Sutcliffe coefficient (R 2 ) and the volume error ratio (VER). Only runoff data prior to the onset of the GLOF of Merzbacher Lake were used for the calculation: where n is the total number of calculated time steps, Q m [m 3 s −1 ] is the observed discharge and Q s [m 3 s −1 ] is the simulated discharge.
The objective function snow-covered area (SCA) is calculated as the percentage of correctly attributed snow-covered and snow-free pixels (px corr ) in the total number of pixels in the catchment (px tot ): SCA = px corr px tot × 100.
As only one year of SCA data was available, it was used in each of the split samples.

Calibration Procedure
The quite complex arrangement of data and model inputs during the different steps of model application in the past and the future is summarized graphically in Figure 3.  At the first calibration step, the number of free model parameters (Table 1) is reduced from 23 to 19 by fixing the less sensitive parameters to reasonable values [54,56]. The automatic calibration procedure considered runoff (R 2 ) and SCA (%). The most promising parameter set of each split sample was selected, called compromise solution in the following [101].
Although runoff and snow cover could be simulated very satisfactory (see Section 3.2. for details) in the first step, the modelled glacier mass balances were far too negative when compared to values from the literature (see Table 4). As this deviation could not be avoided by automatic calibration, we conducted a manual optimization run as a second calibration step. Here, we changed the parameters that control liquid precipitation and glacier mass balance in order to obtain more realistic mass losses. The other model parameters remained unchanged.
As a third calibration step, another automatic calibration run involving only the runoff parameters has been performed. In this step, the runoff simulation is optimized again under the changed precipitation and glacier mass regime. By fixing the other parameters, we force the model to find another optimum in a neighborhood with realistic glacier mass balances.

Climate Scenarios
The results of the climate models are visualized in Figure 4 (annual changes) and Figure 5 (monthly changes). At the first calibration step, the number of free model parameters (Table 1) is reduced from 23 to 19 by fixing the less sensitive parameters to reasonable values [54,56]. The automatic calibration procedure considered runoff (R 2 ) and SCA (%). The most promising parameter set of each split sample was selected, called compromise solution in the following [101].
Although runoff and snow cover could be simulated very satisfactory (see Section 3.2. for details) in the first step, the modelled glacier mass balances were far too negative when compared to values from the literature (see Table 4). As this deviation could not be avoided by automatic calibration, we conducted a manual optimization run as a second calibration step. Here, we changed the parameters that control liquid precipitation and glacier mass balance in order to obtain more realistic mass losses. The other model parameters remained unchanged.
As a third calibration step, another automatic calibration run involving only the runoff parameters has been performed. In this step, the runoff simulation is optimized again under the changed precipitation and glacier mass regime. By fixing the other parameters, we force the model to find another optimum in a neighborhood with realistic glacier mass balances.

Climate Scenarios
The results of the climate models are visualized in Figure 4 (annual changes) and Figure 5 (monthly changes).  Temperature and precipitation of the REMO and RG scenarios differ considerably ( Figure 4, Table 2). For the RG scenarios, the changes compared to past climate conditions are smaller than for REMO, as this approach only provides a lower bound of the possible changes of temperature and precipitation in this region. While precipitation decreases considerably in the REMO projections, the RG projections do not show a clear precipitation trend, though the majority of the RG projections simulate less precipitation for the future climate. Generally, differences between the three RG scenarios are rather small. For temperature, the magnitude of change is also highest for the REMO projections. However, both approaches reveal a strong warming for the future climate. The variations between the RG realizations partly covers the temperature differences between the scenarios, but   Temperature and precipitation of the REMO and RG scenarios differ considerably ( Figure 4, Table 2). For the RG scenarios, the changes compared to past climate conditions are smaller than for REMO, as this approach only provides a lower bound of the possible changes of temperature and precipitation in this region. While precipitation decreases considerably in the REMO projections, the RG projections do not show a clear precipitation trend, though the majority of the RG projections simulate less precipitation for the future climate. Generally, differences between the three RG scenarios are rather small. For temperature, the magnitude of change is also highest for the REMO projections. However, both approaches reveal a strong warming for the future climate. The variations between the RG realizations partly covers the temperature differences between the scenarios, but Temperature and precipitation of the REMO and RG scenarios differ considerably (Figure 4, Table 2). For the RG scenarios, the changes compared to past climate conditions are smaller than for REMO, as this approach only provides a lower bound of the possible changes of temperature and precipitation in this region. While precipitation decreases considerably in the REMO projections, the RG projections do not show a clear precipitation trend, though the majority of the RG projections simulate less precipitation for the future climate. Generally, differences between the three RG scenarios are rather small. For temperature, the magnitude of change is also highest for the REMO projections. However, both approaches reveal a strong warming for the future climate. The variations between the RG realizations partly covers the temperature differences between the scenarios, but they generally lie close together. Compared to the RG projections, the REMO scenario is warmer and dryer in the past period and this effect is further enhanced in the future period.

Calibration Results
During the first calibration step (automatic calibration) using the objective functions R 2 and SCA, the three split samples converged to nearly stable values within 20,000 model evaluations. All Pareto fronts show small trade-offs in comparable ranges for both objective functions. Nash-Sutcliffe coefficients range from 0.71 to 0.94 and SCA range from 76 to 85.5% of pixels with correct assignment to snow-covered and snow-free areas. Based on this calibration, a compromise solution for each sample was chosen (Table 3). Model runs over longer time periods using the selected compromise solution showed that local and overall mass balances are simulated far too negative in the ablation zone and not positive enough in the accumulation area (Table 4). This indicates that underestimated precipitation was compensated by overestimated melt rates. This error compensation has been reversed by the second calibration step, where the corresponding parameters were tuned manually. This second step improved glacier mass balances significantly, but reduced the quality of the runoff simulation. The third step should optimize runoff as much as possible, without losing the realistic glacier mass balance values. For this purpose, another automatic calibration run only considering the runoff routine was conducted. Again, the objective functions R 2 and VER were used. After 20,000 model evaluations, runoff is realistically simulated with Nash-Sutcliffe coefficients above 0.85 and VER of maximum 13% for the calibration year (Table 3), while the calculated mass balances remain within the range reported in the literature ( Table 4). The snow-covered area is still well simulated with more than 82% of pixels with correct assignment to snow-covered and snow-free areas.

Melt-and Runoff Scenarios
We used the Compromise Solution parameter sets of the three split samples to simulate runoff for the Past and Future Periods. The Past Period was simulated with the Past Area as a model input, the Future Period was simulated with Past Area and Best Guess. Additionally, the model inputs of the sensitivity analysis were used. The main results are summarized in Table 5. Compared to the Past Period, the runoff sums increase significantly in the future (Table 5). Since the three realizations of the RG scenarios differ only by standard deviations of 0 to 3%, only their means are presented here. If the past glacier area is assumed, the runoff nearly doubles, with an increase of 91% for REMO and up to 97% for the RG scenario A2. Using the Best Guess model input, runoff still increases, but in smaller ranges between 46% and 54% with a similar ranking of the scenarios. Using the model inputs of the sensitivity analysis, an area reduction of 40.9% still leads to a slightly increased runoff. If the glacier area is reduced by 43.0%, this leads to similar runoff as calculated in Period 1. This means that an area reduction of more than 43.0% is needed to compensate enhanced melt rates in the future climate. After this tipping point, further area recession will cause a reduction of annual discharge.
The mean monthly runoff and ice melt of all scenarios is presented in Figure 6. Due to higher deviations, especially in the spring months from April to June, all realizations of the RG scenarios are presented. In the Past Period, all samples and all scenarios show, more or less, the same runoff characteristics. The gradual decrease of runoff in the winter months stops with a small increase in May, followed by a more significant increase in June and runoff maxima in July and August. The individual runoff curves are controlled by both, the split sample and the climate scenario. Generally spoken, Sample 1 tends to produce runoff peaks in July, but more distinct for the REMO than for the RG scenarios. Sample 2 and 3 show a less distinct runoff peak in July, RG scenarios even show a double peak in July and August. For the Future Period (both, Past Area and Best Guess) and all scenarios and samples, the first runoff increase is now visible in April. This effect is strongest for the scenario A2 and smallest for REMO. In the REMO sets, the runoff peak in July is less distinct or even shifted to August, while for all RG scenarios the runoff peak in July is much more pronounced. In general, the runoff curves of the original and adapted area only differ in runoff amounts, but not in their distribution. The mean monthly runoff and ice melt of all scenarios is presented in Figure 6. Due to higher deviations, especially in the spring months from April to June, all realizations of the RG scenarios are presented. In the Past Period, all samples and all scenarios show, more or less, the same runoff characteristics. The gradual decrease of runoff in the winter months stops with a small increase in May, followed by a more significant increase in June and runoff maxima in July and August. The individual runoff curves are controlled by both, the split sample and the climate scenario. Generally spoken, Sample 1 tends to produce runoff peaks in July, but more distinct for the REMO than for the RG scenarios. Sample 2 and 3 show a less distinct runoff peak in July, RG scenarios even show a double peak in July and August. For the Future Period (both, Past Area and Best Guess) and all scenarios and samples, the first runoff increase is now visible in April. This effect is strongest for the scenario A2 and smallest for REMO. In the REMO sets, the runoff peak in July is less distinct or even shifted to August, while for all RG scenarios the runoff peak in July is much more pronounced. In general, the runoff curves of the original and adapted area only differ in runoff amounts, but not in their distribution.  The differentiation of runoff into its sources rain, snow and ice is shown in Figure 7. In the Past Period, the biggest part of runoff is provided by ice melt, closely followed by snow melt and rain. Highest ice melt was simulated in all Periods by REMO and lowest by RG B1. At Sample 3, more rain and less snow melt are simulated compared to the other samples. In the future periods, ice melt increases considerably while snowmelt decreases and rain increases slightly. As expected, due to the smaller glacier area, the future period with adapted area shows smaller shares of ice melt. The differentiation of runoff into its sources rain, snow and ice is shown in Figure 7. In the Past Period, the biggest part of runoff is provided by ice melt, closely followed by snow melt and rain. Highest ice melt was simulated in all Periods by REMO and lowest by RG B1. At Sample 3, more rain and less snow melt are simulated compared to the other samples. In the future periods, ice melt increases considerably while snowmelt decreases and rain increases slightly. As expected, due to the smaller glacier area, the future period with adapted area shows smaller shares of ice melt. The mean elevation of the ELA ranges from 5210 m a.s.l. in Sample 3 in RG Scenario A2 up to 5283m a.s.l. in Sample 2 at the REMO Scenario (Table 6). The meltwater supplies of clean-ice areas, debris-covered areas and ice cliffs were exemplarily calculated for the past period of the REMO scenario with the Sample 1 parameter set ( Table 7). The percental share of melt below debris is higher than its share of the total glacier area. Ice cliffs only The mean elevation of the ELA ranges from 5210 m a.s.l. in Sample 3 in RG Scenario A2 up to 5283m a.s.l. in Sample 2 at the REMO Scenario (Table 6). The meltwater supplies of clean-ice areas, debris-covered areas and ice cliffs were exemplarily calculated for the past period of the REMO scenario with the Sample 1 parameter set ( Table 7). The percental share of melt below debris is higher than its share of the total glacier area. Ice cliffs only contribute with slightly less than 5% to total ice melt, which is more than six times their areal share. The runoff contribution of ice melt at supraglacial lakes is even smaller than its areal share. As shown in Table 4, regional and total mass balances were simulated satisfactory after the third calibration step. Additional to this calculation, the mass balance calculation using the fictional glacier input without debris cover, ice cliffs and supraglacial lakes are shown in Table 4 (grey shaded mass balance values). For this fictional clean-ice glacier, significantly higher mass losses are simulated for both, regional and overall mass balances.

Discussion
As confirmed already by several authors [97,98], SCA can be used successfully as objective function for model calibration. Here, the SCA simulation is considerably improved and the runoff simulation only slightly deteriorated compared to a calibration using runoff only. Duethmann et al. [98] recommend the use of SCA to obtain higher internal model consistency. In our study, this statement was verified using the validation values from the literature (Table 4). In fact, SCA helps to obtain better results for the ELA, which is simulated more realistically compared to model calibrations with runoff only. However, contrary to expectations, this does not automatically lead to better mass balance simulations. The model still provided too much runoff originating from snow and ice melt. This was also observed in an earlier study [23], where the HBV model underestimated basin precipitation and compensated this by overestimating glacier melt.
Manual parameter adjustment considering values from the literature seems to be a good opportunity to improve automatic calibration results. This procedure combines the advantages of both, automatic and manual calibration and can also be used if only soft data or expert knowledge is available [102]. The resulting deterioration of the runoff simulation could almost be eliminated by the second automatic calibration run for runoff improvement, where only parameters concerning runoff formation were changed.
The final Nash-Sutcliffe coefficients (R 2 ) and volume error ratio (VER)m as well as snow-covered area (SCA) of the three split samples, lie generally within acceptable ranges compared to other studies with a multi-objective calibration approach (e.g., Reference [21]). While in Sample 1, the runoff curve of the validation year is simulated worse, this effect decreases in Samples 2 and 3. Similar results were found by Mayr et al. [54]. In sample two, the parameter set found in comparably wet years was capable to reproduce the hydrograph of a dry validation year. This is an indicator that the model calibration is also valid for simulations outside the rather short calibration period.
Results SCA simulations are slightly deteriorated during the manual refinement, but they are still much better as in model runs without consideration of SCA. Other studies [86][87][88] have calculated SCA not on a raster basis, but for elevation zones and therefore, their results cannot be compared directly.

How Will Runoff and Melt Change in Future?
Future runoff scenarios in glacierized catchments are generally hampered by high uncertainties [103]. The most important factors that determine the development of runoff characteristics over the next decades are future climate, future glacier size and future debris cover. All of them cannot be projected in a straightforward way, but have to be simulated or estimated by more or less complex models and assumptions. Another limitation is that we cannot be sure if the parametrization of the runoff model will still valid in a future climate. Although successful calibration and validation in separate years with different meteorological conditions gives some cause for optimism, our series of hydrometeorological data is too short to enable a real differential split-sample test after Klemes [99]. Being aware that our Best Guess is only one possible future, we are convinced that this future is a very likely one and that a more accurate projection is hardly possible with todays' data and methods. Even if the magnitude of the hydrological response cannot be secured, we have confidence that trends and signals point into the right direction.

Annual Runoff Changes
Even with a glacier area reduction by 22.5%, the runoff will still increase significantly in the future. This is caused by higher melt rates in the ablation area due to higher temperatures, but also by the enlargement of the ablation area, due to the elevation shift of the ELA. A glacier area reduction of 43% is necessary to reduce runoff to the level of the Past Period. For the larger Aksu subbasins of Sari-Djaz and Kaksaal, Duethmann et al. [104] predict an overall increase in discharge for the period 2010-2039 (reference 1971-2000), but a decrease in 2070-2099 caused by an area reduction of 28-89% in the Sari-Djaz and 47-94% in the Kakshaal basin. Huss and Hock [105] also expect the runoff peak for the entire Tarim basin between 2030 (RCP2.6) and 2058 (RCP8.5) and thus a decline afterwards. These deviations from our results are due to the consideration of larger catchments with smaller degrees of glacierization, smaller average glacier size, and to a different methodology of estimating of future glacier extent, resulting in higher area losses.

ELA Changes
In our simulation, the ELA rises by about 800 m to elevations ranging from 5210 to 5283 m a.s.l., depending on sample and scenario. This increases the areal share of the ablation area from 49.9% to 87.9%, if the glacier area remains unchanged, and to 86.5% in case of the adapted area.
Modelling of the mass balance of Tian Shan glaciers until the end of the 21st century by Aizen et al. [106] revealed an mean increase of the equilibrium line altitude by 570 m and a decrease of the glacier covered area by 31%. Temperature is assumed to increase by 4 • C, which is comparable to our data. The deviation of ELA elevation is probably caused by the different assumptions on the precipitation trend, which is assumed to increase by 10% [106]. But even with a lower ELA, a smaller glacier area is assumed [106]. This clearly shows the high uncertainty of future glacier area assumptions. For this reason, we conducted a sensitivity analysis to identify the glacier coverage necessary for decreasing runoff.

Monthly Runoff Changes
Generally, the discharge hydrographs of the Past Period show a more or less pronounced runoff maximum in July and ice melt maximum in August. The hydrographs of the mean monthly runoff sums deviate only slightly in the different samples and scenarios ( Figure 6). The most noticeable difference is the timing of the runoff peaks. Compared to RG, REMO has a lower June and higher July temperature, resulting in a melt-driven runoff peak in July. All samples of the scenarios show an ice melt maximum in August, because at this time, the ablation area is almost free of snow. Here, differences in the timing of runoff peaks can be explained by the used model parameters. For example, Sample 1 shows a higher runoff in July, combined with higher ice melt. Both are caused by a higher snow melt parameter R snow , which increases snow melt in July. The already snow-free glacier area increases the ice melt at the same time. Runoff in August is reduced, due to the already melted snow cover.
In the Future Period, the temperature increase in winter for the REMO scenario does not influence the melt rates because temperature remains below the melting point [8]. The earlier and more intense runoff in spring is caused by more intense snow melt, due to the overall higher spring temperatures, this effect was reported by many authors (e.g., by Hagg et al. [21] for the Rukhk catchment in the Pamir or by Akhtar et al. [107] for glaciers in the Hindukush-Karakorum-Himalaya region). In summer, the runoff curves changed differently for the REMO and RG scenarios. With REMO, the runoff peak in July is slightly lowered (Figure 6), because higher temperatures cause a shift of the snowmelt maximum from July to June. In contrast, the runoff peak in July is even more pronounced for the RG projections. The existing shift of snow melt is superimposed by an increase of ice melt in July, due to the temperature maximum in July in the RG scenarios. The deviations between the samples are similar to those of the Past Period. In general, the differences of runoff and ice melt between the different samples are relatively small. This supports the reliability of our modelling results.
Sorg et al. [83] suppose that continuing glacier mass losses in Tian Shan will finally lead to a change from glacial-nival to nival-pluvial runoff regimes, which will also increase the year-to-year variability of runoff. The timing of this change is highly dependent on the individual catchment and its glacier area development. This is not yet reached in our scenarios, as Northern and Southern Inylchek are comparatively large and therefore glacier melt still provides the main part of runoff.

What Is the Contribution of Snow, Ice and Precipitation to Runoff?
The share of ice melt in total runoff of the Past Period is 37-45%. Considering the high glacierization of more than 70%, this value seems small compared to the 36-38% calculated for the whole Aksu catchment (glacierization: 18%) by Aizen and Aizen [9]. In fact, these numbers cannot be compared properly, because Aizen and Aizen [9] consider all runoff from glaciers, including snow melt and liquid precipitation. It also should be noted that in very arid lowlands, the fraction of melt that is determined in the mountainous part of the basin, does not change anymore downstream. Comparing our findings to results from other catchments is even more difficult, because the fraction of ice melt strongly relies on climate. For example, the Rukhk catchment in the Pamir yields a significant share of ice melt of 32% with a glacier coverage of only 10%. This is caused by a very arid climate with an annual precipitation of only 295 mm [21]. In contrast, the Langtang Khola catchment in the Himalaya has glacier coverage of 43.5%, but the share of glacier runoff share was estimated to be only about 18%, due to its monsoon-influenced climate [108].
In the Future period, the share of ice melt increases significantly, due to the expected temperature increase. This warming will also cause more liquid precipitation, increasing the share of rain and decreasing snow melt. Enhanced ice melt will clearly overcompensate the decrease in glacier area and is expected to double until the end of this century. The timing of the tipping point, from which on area losses dominate and runoff decreases again, is mainly controlled by the initial glacier size and volume.

How Important Are Debris Cover, Ice Cliffs and Supraglacial Lakes for Melt Rates at
Debris-Covered Glaciers? Table 7 shows that the share of sub-debris melt in total runoff is higher than the fraction of debris covers in total glacier area. Since we assume that the presence of supraglacial debris causes an overall reduction of melt rates, this result may be surprising. It is caused by the fact that debris covers are situated at the lower parts of the glacier tongue, where air temperatures and resulting melt rates are higher than close to the equilibrium line, where bare ice dominates. The protecting effect of debris is confirmed by our test run using a fictional bare ice glacier (Table 4, grey shaded mass balance value). Here, simulated ice melt increases significantly compared to the original, debris-covered glacier. This supports the finding of Zhang et al. [32] at nearby Keqicar Baqi Glacier in the Chinese Tian Shan and Ragettli et al. [34] on debris-covered glaciers in the Langtang Himal, Nepal and conflicts with the results of other studies [36][37][38].
But even though thick debris covers have an insulating effect, considerable downwasting at heavily debris-covered glaciers has been observed, especially in parts which are interspersed with ice cliffs and supraglacial lakes [83,109,110]. The melt in supraglacial lakes is most likely underestimated in our results as lateral melt and calving could not be considered. The increased melt at ice cliffs, as stated by several authors [33,41,42,44,111] is confirmed by our findings (Table 7). Considering the small area covered by ice cliffs, their share in total ice melt is significant. However, according to our results, the enhanced ice melt at ice cliffs on Inylchek glaciers cannot compensate the overall isolating effect of the debris cover.

Conclusions
This study shows that additional objective functions, such as snow-covered area, improve the simulation quality of the respective variable, but do not guarantee valid model calibrations. Simulation results always need to be checked carefully on their plausibility, using available field data or expert knowledge. Considering these kinds of data during manual parameter adjustment is a good opportunity to improve automatic calibration results.
Using both, dynamical and statistical-dynamical modelling approaches to assess regional climate reveals significant projection differences. While the former estimates an upper bound for the possible changes in temperature and precipitation over the region, the latter provided a lower bound and considered multiple SRES scenarios, thus permitting an estimation of uncertainties. This emphasizes the need to meet uncertainties, due to climate modeling and the choice of the downscaling approach by including several downscaling realizations.
From 1970/1971 to 1999/2000, the share of ice melt in total runoff was 37-45%. By the end of this century, this already high share of ice melt, as well as total runoff, will still significantly increase, due to strongly enhanced ice melt in the warmer atmosphere and comparably small area losses of the heavily debris-covered Inylchek glaciers. As ice melt will still be the main source of runoff in this highly glacierized subbasin, a change from a glacial-nival to a nival-pluvial runoff regime will not be reached in this century according to our model results.
A major part of ice melt is provided by debris-covered glacier parts. This doesn't contradict the assumption of a protecting effect of debris covers, but is caused by the spatial distribution of debris-covered and clean-ice glacier areas. Simulations at a fictional clean-ice glacier confirm the protecting effect of a debris cover. Compared to the original debris-covered glacier, ice melt is increased significantly if a clean ice surface is assumed.
At ice cliffs, melt is significantly increased, but it is not high enough to compensate for the isolating effect of the debris cover. Melt in supraglacial lakes is most likely underestimated because lateral melt and calving could not be considered in the model approach. These processes should be implemented into the model in the future in order to further improve predictions of water availability in this region.
Although the projected increase in runoff will counteract man-made water shortages downstream and probably mitigate the negative ecological impacts, one should keep in mind that the additional water yield from glacier wastage is preliminary and not of sustainable nature. Acknowledgments: The German Research Centre for Geosciences (Geoforschungszentrum, GFZ) in Potsdam and the Central Asian Institute for Applied Geosciences (CAIAG) in Bishkek offered essential logistic help during the field campaigns. Gleb Glazirin, Tashkent, supplied important hydrometeorological data and Jasper Vrugt (University of California, Irvine) kindly provided his calibration algorithm. The comments of three anonymous reviewers are greatly acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.