Emptying Water Towers? Impacts of Future Climate and Glacier Change on River Discharge in the Northern Tien Shan, Central Asia

: Impacts of projected climate and glacier change on river discharge in ﬁve glacierized catchments in the northern Tien Shan, Kazakhstan are investigated using a conceptual hydrological model HBV-ETH. Regional climate model PRECIS driven by four di ﬀ erent GCM-scenario combinations (HadGEM2.6, HadGEM8.5, A1B using HadCM3Q0 and ECHAM5) is used to develop climate projections. Future changes in glaciation are assessed using the Blatter–Pattyn type higher-order 3D coupled ice ﬂow and mass balance model. All climate scenarios show statistically signiﬁcant warming in the 21st Century. Neither projects statistically signiﬁcant change in annual precipitation although HadGEM and HadCM3Q0-driven scenarios show 20–37% reduction in July–August precipitation in 2076–2095 in comparison with 1980–2005. Glaciers are projected to retreat rapidly until the 2050s and stabilize afterwards except under the HadGEM8.5 scenario where retreat continues. Glaciers are projected to lose 38–50% of their volume and 34–39% of their area. Total river discharge in July–August, is projected to decline in catchments with low (2–4%) glacierization by 20–37%. In catchments with high glacierization (16% and over), no signiﬁcant changes in summer discharge are expected while spring discharge is projected to increase. In catchments with medium glacierization (10–12%), summer discharge is expected to decline under the less aggressive scenarios while ﬂow is sustained under the most aggressive HadGEM8.5 scenario, which generates stronger melt.


Introduction
High mountains are often referred to as water towers of the world and this is particularly true in Central Asia where snow and glacier melt nourishes rivers providing water to the arid plains and sustaining agricultural production on the irrigated land [1].Various components of the mountain cryosphere contribute to the formation of discharge including melting snow pack, glacier ice, permafrost and ice contained in rock glaciers [2].All these components are affected by climate change.Observations suggest a decrease in the extent, depth, and duration of mountain snow pack [3].There is ample evidence for declining glacier mass balance [4] and wastage of glaciers across the Tien Shan mountains [5][6][7][8][9][10]. Results of both, observations and modelling, show that the depth of active layer of mountain permafrost is increasing while its spatial extent is declining [11].These changes affect river discharge and it is important to investigate both, the observed and projected changes in discharge in order to quantify the effects of changes in the cryosphere on water supply and establish the timing of peak flow after which discharge may decline in response to the diminishing snow and ice.
To date, relatively few studies assessed the observed and projected changes in discharge in glacierized catchments of Central Asia largely because discharge data in natural catchments, located at higher elevations, are not widely available following the closure of many gauging sites in the 1990s [2,8].Those studies that analyzed the long-term (over 50 years) discharge records in the headwater catchments show that to date, there has been no decline in the melt season discharge since the 1950s-1960s in the central [12], eastern [13], and northern [2] Tien Shan.Over shorter time periods, e.g., after the 1970s, an increase in the melt season discharge was registered in the catchments with glacierization exceeding approximately 10% of the total catchment area [2].
A number of modelling studies, aiming at the assessment of sensitivity of discharge to glacier shrinkage or at the development of future hydrological projections were undertaken for the Tien Shan either as part of a global assessment [14] or on a basin scale.The outcomes of these studies, their limitations and innovations are comprehensively reviewed by [15] and a further review is provided by [16].The limitations can be summarized as follows.Firstly, in a range of models, glacier melt is not represented although the models are applied to glacierized catchments.Secondly, streamflow alone is often used for model validation although the use of other objective functions (e.g., mass balance, snow water equivalent (SWE)) helps to constrain uncertainty [17].Thirdly, availability of the observational input data is restricted and data quality is often low.The lack of high-elevation meteorological stations and inhomogeneity of both meteorological and hydrological records are a well-known problem in the Tien Shan [8] restricting modelling efforts and potentially affecting the output.Temperature gradients are not constant throughout the year and precipitation gradients are not uniform in the mountain regions [18] and, therefore, interpolating data from the lowland stations generates additional uncertainty.Fourthly, future climate projections are usually derived from the CMIP3 and CMIP5 ensembles whereby spatial resolution of modelled climate data varies between 0.75 • and 3.25 • resulting in a considerable bias in input data in the complex terrain [19].The use of regional climate models (RCM) has been so far restricted in the Tien Shan although several studies produced dynamically downscaled data [20][21][22].Fifthly, there is a shortage of transient projections of glacier change [15] which are mostly limited to the large-scale regional assessments [23,24].Few hydrological studies employ transient scenarios of glacial change in Central Asia in contrast to the Hindukush-Himalaya region [25,26] except an assessment of future changes in the extent of glaciers and its impacts on water resources in the headwaters of the Amu Darya and Syr Darya [25] and global assessment by [14].Otherwise, sensitivity analyses are conducted instead whereby no change, complete disappearance or arbitrarily set glacier area reduction are used [27][28][29] or changes in glacier area, derived from parametrizations based on changes in mass balance and equilibrium line altitude (ELA) are projected for a specific time slice [30].This approach is useful in assessment of future water resources, however, it does not address the issue of timing of peak flow.
This paper presents future hydrological scenarios for five natural, glacierized catchments located in the northern Tien Shan (Table 1; Figure 1).All catchments are located in proximity to each other and in a relatively data-rich region where high-elevation meteorological, glacier mass balance and SWE data are available.The size, elevation span, and glacierization are different predetermining different responses of streamflow to climate change [2].The widely used HBV-ETH model [30][31][32] is applied and validated using additional objective functions.Hydrological modelling is performed for a range of climate scenarios generated using high-resolution regional climate model PRECIS and deglacierization scenarios, developed using the combined ice flow and mass balance models and projecting changes in area and volume of glaciers.and SWE data are available.The size, elevation span, and glacierization are different predetermining different responses of streamflow to climate change [2].The widely used HBV-ETH model [30][31][32] is applied and validated using additional objective functions.Hydrological modelling is performed for a range of climate scenarios generated using high-resolution regional climate model PRECIS and deglacierization scenarios, developed using the combined ice flow and mass balance models and projecting changes in area and volume of glaciers.
Figure 1.Study area.Letters refer to the location of meteorological stations: A-Chimbulak, B-Mynzhilki, C-Tuyuksu, D-Bolshoe Almatinskoe Lake.Flow gauging sites are numbered as in Table 1.

Study Area
The five modelled catchments of the rivers Kishi Almaty, Ulken Almaty, Turgen, Talgar, and Kaskelen are located on the northern slope of the Ile (Zailiiskiy) Alatau Range (Figure 1).The region is characterized by strong seasonal variations in temperature and precipitation.The westerly flow dominates in autumn and spring resulting in the precipitation maxima in April-May on the plains shifting towards June-July in the middle and high mountains, where snow accumulation peaks in spring and early summer.In winter, the western extension of the Siberian anticyclone predetermines sub-zero temperatures and small amounts of solid precipitation.In summer, the thermal Asiatic depression dominates driving advection from the south which results in hot and dry weather on the plains and high demand for irrigation water [33].
The Ile Alatau is extensively glacierized and in 2008, glaciers occupied 565 km 2 in the Kungei-Ile Alatau system [34].Between 1990 and 2008, glaciers were losing their area at an average rate of 0.9%  1.

Study Area
The five modelled catchments of the rivers Kishi Almaty, Ulken Almaty, Turgen, Talgar, and Kaskelen are located on the northern slope of the Ile (Zailiiskiy) Alatau Range (Figure 1).The region is characterized by strong seasonal variations in temperature and precipitation.The westerly flow dominates in autumn and spring resulting in the precipitation maxima in April-May on the plains shifting towards June-July in the middle and high mountains, where snow accumulation peaks in spring and early summer.In winter, the western extension of the Siberian anticyclone predetermines sub-zero temperatures and small amounts of solid precipitation.In summer, the thermal Asiatic depression dominates driving advection from the south which results in hot and dry weather on the plains and high demand for irrigation water [33].
The Ile Alatau is extensively glacierized and in 2008, glaciers occupied 565 km 2 in the Kungei-Ile Alatau system [34].Between 1990 and 2008, glaciers were losing their area at an average rate of 0.9% year −1 although the rates vary between the catchments depending on their elevation and size of glaciers ([34]; Table 1).
The snow and glacier melt period is limited to June-July-August (JJA) extending to September in individual years.The seasonal flow cycle is driven by snow melt in June-July and glacier melt in August [2].
Characteristics of the five study catchments are given in Table 1 while a more detailed description of the catchments is given in [2].The degree of glacierization varies from 16% in the Talgar and Ulken Almaty to 2.6% in the Kaskelen.

Observational Data
Several observational data sets were used to calibrate and validate the models (Table 2).Glaciological and meteorological data were obtained from the Tuyuksu group of glaciers (source of the Kishi Almaty) and Tuyuksu meteorological station (Figure 1).The Tuyuksu data were used to calibrate HBV-ETH for all catchments as this is the only high-elevation station in the region for which continuous and uninterrupted daily measurements were available.While climatic conditions and particularly precipitation vary strongly in the mountains [18], the spatial variability is reduced at higher elevations (e.g., ~3500 m a.s.l or 700 hPa) and we consider Tuyuksu representative of the high-elevation climate of all catchments [2].HBV-ETH Tuyuksu (Figure 1) is one of the world reference glaciers whereby mass balance and other glaciological measurements have been conducted by Kazakhstan Institute of Geography and reported to the World Glacier Monitoring Service (WGMS) since 1957.It was previously shown that the observed changes in area of the Tuyuksu glacier correlated strongly with changes in glacier area in the Ile Alatau as a whole [34].Accumulation and ablation are measured by the glaciological method using a network of approximately 100 stakes providing data for the estimation of summer and winter mass balance, ELA, and accumulation area ratio (AAR).Glacier flow velocities are measured using the same set of stakes, whose positions are recorded using GPS at the start and the end of winter and summer mass balance seasons.
Tuyuksu meteorological station is located in proximity to the Tuyuksu glacier (Figure 1; Table 2) and provides daily mean air temperature (derived from 3-hourly measurements) and daily precipitation (derived from precipitation measured twice a day).Changes in the summer and winter mass balance, ELA, AAR as well as seasonal temperature and precipitation at the Tuyuksu glacier and station respectively are reviewed in [2].
Data on snow water equivalent (SWE), measured daily, were obtained from the Chimbulak avalanche monitoring station located in the Kishi Almaty valley (Figure 1; Table 2).
Daily streamflow values were obtained for the gauging sites listed in Table 1.A detailed description of flow measurements in the selected catchments including measurement practices, compilation of the data sets and analysis of data quality are provided in [2].We note here that all streamflow data are for the natural sections of the catchments where there was no channel modification or water abstraction and where homogeneous flow records were available.The Talgar catchment is the only one where there was a lengthy interruption in flow measurements between 1994 and 2005, however, continuous observations were resumed at the same location using the same methods in 2005.
Glacier areas were measured in the study catchments using Landsat 5 imagery from 1990 and 2006 (Table 2) to provide input data for HBV-ETH.A semi-automated classification using bands 4 and 5, followed by manual corrections (where required), was applied.We note that changes in glacierized area in the Ile Alatau were extensively reviewed in [6,34] and will be addressed here solely in the context of changing streamflow.
Annual data on the area of the Tuyuksu glacier, required for glaciological modelling, was derived thorough the in situ GPS measurements which are a part of the Tuyuksu monitoring programme.

Climate Modelling
The baseline and future climate change scenarios were generated using Regional Climate Model (RCM) PRECIS (PREdicting Climate for Impact Studies) developed by the UK Met Office [35].PRECIS is a hydrostatic model with a horizontal resolution of 0.22 • (~25 km).Factional grid box land cover is used to improve the detail of surface characterization with two land cover types representing over 50% and over 25% of each grid box.PRECIS derives lateral boundary conditions from a range of GCM.
In this study, four integrations were used forced by three GCM including (i) HadGEM for Representative Concentration Pathway (RCP) 2.6 and 8.5 scenarios from CMIP5 and (ii) HadCM3Q0 and ECHAM5 for A1B scenario from CMIP3 ensembles.RCP 2.6 and 8.5 refer to the assumed radiative forcing of 2.6 W m −2 (least climatic warming) and 8.5 W m −2 (most aggressive climatic warming) at equilibrium.The A1B scenario represents the world in which energy production is balanced across all sources resulting in medium-to-strong warming.These combinations of scenarios and GCM were selected to represent the 'warm and dry' (HadCM3Q0) and 'cold and wet' (ECHAM5) models as well as a weak (RCP 2.6) and strong (RCP 8.5) warming.We note that model description as 'warm' or 'cold' refers to model sensitivity and not to over-or underestimation of regional or local temperatures.
The baseline climate was simulated for the 1981-2005 period and future climate for the 2006-2095 period.The model domain extended between 32-55 • N and 46-86 • E and encompassed the plains of Kazakhstan and Central Asia, the Tien Shan and the Pamir mountains.The present-day land cover was assumed in future climate simulations.The boundaries of the Aral Sea were adjusted manually on the seasonal basis.

Validation of Climate Models
Daily mean air temperature values and precipitation rates, simulated by PRECIS, were obtained for the grid points positioned over the study catchments and averaged to produce a single temperature and a single precipitation time series for all study catchments.These data were validated against the observed Tuyuksu data.
The modelled temperature data showed a good agreement with observations in reproducing the annual cycle.The differences in the observed and modelled values were largely due to the difference in Water 2020, 12, 627 6 of 31 elevation between the Tuyuksu station and PRECIS grid points although larger bias characterized the ECHAM5-driven simulation.The modelled precipitation data exhibited stronger bias overestimating the cold season and underestimating summer precipitation (Figure 2c).The study catchments are located on the northern slope of the Tien Shan and extend from the maximum altitude of 4400 m to the plain.There is a progressive shift in the timing of the precipitation maximum from spring to mid-summer with elevation and above 3000 m a.s.l.precipitation peaks in June-July while on the plains it peaks in April-May.In the uncorrected PRECIS simulations, seasonal distribution of precipitation is consistent with precipitation regime observed on the plain.

Glaciological Scenarios
Evolution of glacier area and volume were modelled using a Blatter-Pattyn type of a higherorder (HO) 3-D ice dynamics model whose structure, equations and application are discussed in [44][45][46][47].This is a modified version of the original Blatter-Pattyn model [44] with the modifications outlined in [45].The HO 3-D model is coupled to a surface mass balance model based on the energy- Overestimation of cold season precipitation is evident in all simulations and particularly in the ECHAM5-driven simulation while the underestimation of summer precipitation is evident in the HadCMQ0-and HadGEM-driven simulations.
The QM methods correct distributions of the modelled data in such a way that they match the distributions of the observed data and the difference between the methods is in distributions and thresholds they use.An advantage of these methods is robust correction of extreme values which makes them particularly suitable for correcting precipitation data.Thus QM methods usually eliminate the 'drizzle effect' which results in a low number of dry days in the modelled data and enable a better correction of extremely high precipitation whose frequency and value are often exaggerated in modelled data [40].LOCI corrects means as well as the wet-day frequencies and intensities by setting wet-day precipitation thresholds in such a way that their exceedance in the modelled data matches the wet-day frequency in observations [39,41].PTP, which is often applied to LOCI-corrected data, further adjusts variance statistics in precipitation time series [42].
Figure 3 presents a comparison of mean monthly temperature values the simulation driven by ECHAM 5 model.The EQM and VARI methods performed best in the correction of temperature projections in this study in line with previous comparisons for western China [43].The EQM-corrected data were selected to force the hydrological model as this method showed a marginal advantage over the VARI method.

Glaciological Scenarios
Evolution of glacier area and volume were modelled using a Blatter-Pattyn type of a higherorder (HO) 3-D ice dynamics model whose structure, equations and application are discussed in [44][45][46][47].This is a modified version of the original Blatter-Pattyn model [44] with the modifications outlined in [45].The HO 3-D model is coupled to a surface mass balance model based on the energybalance approach [47][48][49], which simulates spatial and temporal distribution of accumulation and ablation.The models operate on a grid with a spatial resolution of 25 m.The HO 3-D model has a temporal time step of approximately 1 week (0.02 year) while mass balance is updated annually at The published comparisons of the performance of these methods [39,43] applied to precipitation suggest that QM and PT methods perform best in correcting frequency distributions in the modelled data while LOCI results in a closer match between the modelled and observed records with higher coefficients of determination and Nash-Sutcliffe Efficiency (NSE) index.These methods, however, are mostly aimed at the correction of bias associated with frequency and intensity of wet days.The main deficiency of the modelled precipitation in the study region is a rapid decline of precipitation following the well-reproduced spring maximum and underestimation of precipitation during the dry summer season characteristic of all PRECIS-GCM combinations (Figure 2c).Delta method was the only method that corrected this bias (Figures 2d and 3) and was applied.
DownscaleR, an R-based package for bias correction (www.meteo.unican.es/en/node/73440), was used for correcting both temperature and precipitation.

Glaciological Scenarios
Evolution of glacier area and volume were modelled using a Blatter-Pattyn type of a higher-order (HO) 3-D ice dynamics model whose structure, equations and application are discussed in [44][45][46][47].This is a modified version of the original Blatter-Pattyn model [44] with the modifications outlined in [45].The HO 3-D model is coupled to a surface mass balance model based on the energy-balance approach [47][48][49], which simulates spatial and temporal distribution of accumulation and ablation.The models operate on a grid with a spatial resolution of 25 m.The HO 3-D model has a temporal time step of approximately 1 week (0.02 year) while mass balance is updated annually at the end of each hydrological (mass balance) year.The continuity equation links ice dynamics and the surface mass balance enabling calculation of ice thickness [47,50].
The HO 3-D model was applied to the observed glacier geometry (including ice thickness) to simulate changes in area and volume of the Tuyuksu glacier.The input data for the ice flow model included (i) ice thickness derived from the measurements conducted in 2013 using ground-penetrating Water 2020, 12, 627 8 of 31 radar (GPR) [51] and (ii) DEM derived from the void-filled 30 m SRTM DEM and resampled to 25 m resolution.Meteorological input data, used in the mass balance model, included daily temperature and precipitation from the Tuyuksu meteorological station (Section 3.1; Table 2).In the mass balance model, air temperature controls snow and ice melt (mass discharge from the surface of a glacier) and the ratio between liquid and solid precipitation.Solid precipitation comprises mass accumulation on the surface of a glacier.
The ice flow model was calibrated against the observed values of ice flow velocities derived from the annual displacement of ablation stakes on the Tuyuksu glacier (Section 3.1; Table 2).The key parameters of the model were selected to achieve the best fit between simulated and observed ice flow velocities in the coinciding or closest grid points for fixed glacier geometry.To constrain surface velocity we used an approach described in [46] which is based on the least square fit of the simulated annual mean velocities to the observed annual displacement of ablation stakes.The outcomes of the model calibration were specific mass balance and surface flow velocity field simulated with a 25 m spatial resolution.The modelled surface flow velocities are shown in Figure 4.As expected, the simulated pattern shows maxima on the steepest slopes in the accumulation zone of the Tuyuksu and its tributary and low velocities in the ablation zone.
Water 2020, 12, 627 8 of 32 the end of each hydrological (mass balance) year.The continuity equation links ice dynamics and the surface mass balance enabling calculation of ice thickness [47,50].
The HO 3-D model was applied to the observed glacier geometry (including ice thickness) to simulate changes in area and volume of the Tuyuksu glacier.The input data for the ice flow model included (i) ice thickness derived from the measurements conducted in 2013 using groundpenetrating radar (GPR) [51] and (ii) DEM derived from the void-filled 30 m SRTM DEM and resampled to 25 m resolution.Meteorological input data, used in the mass balance model, included daily temperature and precipitation from the Tuyuksu meteorological station (Section 3.1; Table 2).
In the mass balance model, air temperature controls snow and ice melt (mass discharge from the surface of a glacier) and the ratio between liquid and solid precipitation.Solid precipitation comprises mass accumulation on the surface of a glacier.
The ice flow model was calibrated against the observed values of ice flow velocities derived from the annual displacement of ablation stakes on the Tuyuksu glacier (Section 3.1; Table 2).The key parameters of the model were selected to achieve the best fit between simulated and observed ice flow velocities in the coinciding or closest grid points for fixed glacier geometry.To constrain surface velocity we used an approach described in [46] which is based on the least square fit of the simulated annual mean velocities to the observed annual displacement of ablation stakes.The outcomes of the model calibration were specific mass balance and surface flow velocity field simulated with a 25 m spatial resolution.The modelled surface flow velocities are shown in Figure 4.As expected, the simulated pattern shows maxima on the steepest slopes in the accumulation zone of the Tuyuksu and its tributary and low velocities in the ablation zone.Following model calibration, the bias-corrected temperature and precipitation scenarios (Section 3.2) were assimilated by the coupled ice flow-mass balance model to simulate future mass balance and geometry of the glacier.The assimilation procedure is described in [47].In these experiments, potential changes in net radiation were neglected because prognostic skills of radiation components projected by RCM are relatively low at the temporal and spatial scales required by the coupled glacier model [50].Following model calibration, the bias-corrected temperature and precipitation scenarios (Section 3.2) were assimilated by the coupled ice flow-mass balance model to simulate future mass balance and geometry of the glacier.The assimilation procedure is described in [47].In these experiments, potential changes in net radiation were neglected because prognostic skills of radiation components projected by RCM are relatively low at the temporal and spatial scales required by the coupled glacier model [50].
The integration started on 1 January 2013 (when ice thickness data were obtained by the GPR survey [51]) and ended on 30 December 2095.The integration did not start from the equilibrium state of the glacier which is a preferable practice [50].This is because Tuyuksu is far away from its equilibrium state and both annual and cumulative mass balance of the glacier has been persistently negative [2].Therefore, prognostic changes in area and volume intrinsically include response of the glacier not only to the current and future but also to the past climate fluctuations [52].This approach is Water 2020, 12, 627 9 of 31 justified for two reasons.Firstly, the currently observed changes in the extent of glaciers involve their response to the past climate variations.Secondly, the purpose of this study is to compare evolution of the glacier under different climate scenarios all of which have very similar baseline climates and, therefore, past evolution of the glacier geometry is the same or close under each climate scenario.
The output data were aggregated by the model over four month periods.Annual mean values, derived from the four month data, were used as input to the hydrological model.

Streamflow Modelling: Model Description, Input Data, and Experiment Design
The HBV model is a classic conceptual runoff [53] model which exists in several versions.Already in the 1990s, a glacier routine has been implemented into the HBV structure, now referred to as "HBV-ETH" model [31].Since then, further developments, especially in melt calculation and calibration, have been implemented and the model was applied to glacierized catchments in the Alps [54], the Rocky Mountains [55], the Hindukush-Karakoram-Himalayas [29,31], the Tien Shan [16,27,28], the Pamirs [30], and the Caucasus [56].
A full description of the HBV-ETH model version used in this study is given in [17].Figure 5 illustrates the model structure.The snow and glacier melt are calculated using temperature index method and meltwater refreezing in snowpack is accounted for by the use of a refreezing coefficient which is applied when temperature falls below the melt threshold.Solid and liquid precipitation are distinguished using vertical temperature gradients and a threshold temperature while basin precipitation is calculated using a seasonally stable precipitation gradient.Rainwater and meltwater provide input for the soil moisture routine while actual evapotranspiration is calculated from potential evaporation dependent on the soil moisture content.

Calibration and Validation of the Hydrological Model
The HBV-ETH model was calibrated in a semi-automated mode using the NSE criterion [57], 2 The input data, used for model calibration, consist of topographical data, glacierized area and daily temperature and precipitation.Daily streamflow data are used for model calibration and validation.SWE data are used as an objective control function to constrain the model.The model output consists of daily streamflow, basin precipitation, evapotranspiration and storage change for glaciers, snow and groundwater.
Details of the hydrological and meteorological data are provided in Section 3.1 and in Table 2.The topographical input data included total and glacierized area and aspect presented separately for three orientation classes (north, south, and east-west-horizontal) for 200 m elevation bands.The latter is required to correct for the exposition-dependent differences in solar energy input.Terrain with slopes smaller than 5 • was classified as horizontal.DEMs were derived from the void-filled SRTM DEM data for each study catchment.Glacier masks for 1990 and 2006 were used (Section 3.2).
The simulations include model calibration and validation using observational data followed by four simulations of streamflow for each catchment using the bias-corrected climate scenarios and deglacierization scenarios.The baseline period was set to 1980-2005 in line with the baseline period used in the climate models and future flow projections were developed for 2006-2095.

Calibration and Validation of the Hydrological Model
The HBV-ETH model was calibrated in a semi-automated mode using the NSE criterion [57], coefficient of determination R 2 and mean model error (ME) [58] to characterize the agreement between the measured and simulated daily streamflow.These metrics were calculated for mean flow as well as low and high flow.We note that we use Q90 (flow which was equaled or exceeded for 90% of the specified term) and Q10 (flow which was equaled or exceeded for 10% of the specified term) as indicators of low and high flow, respectively.In addition, SWE was used as objective function comparing both daily and mean monthly values of SWE from Chimbulak station (Table 2) with those simulated by HBV-ETH.
Model calibration was performed using glacierization data from 1990 and 2006 for the two sets of hydrological years, 1985-1988 and 2001-2005, respectively except the Talgar catchment where data for 2007-2010 hydrological years were used because of the gap in measurements.Model calibration was performed in two stages.In order to increase the chance to find a global optimum, a Monte Carlo simulation with broad parameter limits and a high number of model runs, was performed as the first stage.Here, 120,000 parameter sets were created randomly.These model runs partly show high Nash-Sutcliffe efficiencies (>0.9), but values of SWE were not always realistic.
At the second stage, parameter values were restricted to narrower ranges which were in line with previous modelling experiments in Central Asia [9,27,30] and observations and analysis of model uncertainty (Section 3.6).Restricting model parameters and introducing a realistic simulation of SWE as additional calibration criterion, reduced the ability of the model to reproduce streamflow.At this stage, the model parameters were re-tuned in order to achieve a satisfying representation of both runoff and SWE values.The best sets of parameters achieved at this stage were used in all subsequent model runs.

Model Sensitivity and Uncertainty Analysis
The model sensitivity and uncertainty analyses were performed using the model outputs resulting from (i) changes in individual model parameters and (ii) the use of various combinations of model parameters.In order to reduce the number of simulations, sensitivity analysis was conducted by running HBV-ETH for a shorter period between 2007-2008 and 2012-2013 hydrological years.Variations in temperature gradient (TGRAD), rainfall correction factor (RCF), storage discharge constants (k 0 ; k 1 ), field capacity (FC) and parameter describing percolation from upper to lower storage of water (CPERC) affected model output to a larger extent than changes in other parameters.This information was used at the Stage 2 of model calibration (Section 3.5) and in uncertainty analysis using the Generalized Likelihood Uncertainty Estimation (GLUE) method [59,60].GLUE is a common method to assess uncertainty in hydrological model predictions and provides insight into model sensitivity to combinations of model parameters enabling improvements in model calibration.
GLUE approach considers outputs of hydrological modelling for each time step not as a single numerical value but as an interval defined by the prediction bounds reflecting the uncertainties due to the input data.In this study, different sets of model parameters were used to simulate streamflow using NSE as a likelihood measure.Based on the previous studies [60,61], the NSE threshold was set to 0.5 for the selection of simulations to be used in GLUE analysis discarding simulations with lower NSE values.The upper and lower prediction bounds for the simulated streamflow were set as 5th and 95th likelihood-weighted quantiles.
We used containing (containment) ratio (CR), the most commonly applied metric in GLUE analysis defined as a percentage of the observed streamflow values enveloped by the prediction bounds.A high CR for the estimated prediction bounds represents a good model fit with 100% CR being an ideal result [62].

Hydrological Simulations Using PRECIS Data for the Baseline Period and Future Scenarios
Discharge was simulated for the baseline period and future hydrological scenarios were generated using the bias-corrected temperature and precipitation data from four PRECIS-GCM simulations (Section 3.2).Streamflow during the baseline period (1980-2005) was simulated using the sets of parameters obtained with the 1990 and 2006 glacier masks for the 1980-2000 and 2001-2005 periods, respectively.We note that all assessments of future changes in discharge were based on comparisons with discharge in the baseline period simulated using PRECIS data.This approach ensures consistency of climate data and helps to avoid uncertainty due to differences between the observed and modelled climate.
Future data on catchment glacierization were derived from the projections of glacier area change developed for four climate scenarios (Section 3.3).Computational limitations and a lack of ice thickness and glacier flow velocities measurements outside the Tuyuksu group of glaciers did not allow us to model glacier area change in all study catchments.We, therefore, applied elevation-dependent rates of area change, projected for Tuyuksu (source of the Kishi Almaty), to other catchments.This approach is a source of uncertainty in the development of hydrological projections because glaciers of different types positioned at different elevations have a different response time to climate change [5] although it was previously shown that the observed rates of change of Tuyuksu glacier are strongly correlated with rates of glacier area change in the whole of the Ile Alatau [34].
Prior to running future simulations, two approaches to changing glacier area were compared: (i) shrinking glaciers from below and (ii) elevation-dependent glacier area change.For the Ulken Almaty catchment, which together with the Talgar catchment, has the highest glacierization (Table 1), there was no difference in the mean July streamflow values simulated using two sets of the projected glacier area data in the 2006-2020 period while after 2040, the differences were 5-10% of the mean July flow with the altitude-dependent changes generating higher flow.Subsequently, glacier area changes were considered as a function of altitude whereby the rates of relative (%) change simulated for Tuyuksu were applied to the same elevation bands in the modelled catchments.3 show future air temperature and precipitation scenarios for the study region.All scenarios show a statistically significant increase in air temperature.The HadCM3Q0-driven A1B scenario and HadGEM RCP 8.5 show the strongest warming, the former in summer and the latter across all seasons with summer temperatures increasing by 6-7 • C and 7-8 • C respectively (Table 3).The HadGEM-driven RCP 2.6 scenario shows a significant warming in the first quarter of the 21st Century in comparison to the baseline period in all seasons followed by a smaller increase in summer temperature.A brief analysis of the non-downscaled CMIP5 simulations showed that most simulations project temperature stabilization around 2040.Under this scenario, the July and August monthly temperatures increase by approximately 2.5 • C in the 2006-2025 period in comparison with the 1980-2005 period with a further increase of approximately 1 • C towards the end of the 21st Century.The ECHAM5-driven A1B scenario shows the smallest increase in summer temperature but a stronger warming during the cold season between November and March.statistically significant changes in cold season precipitation with an increase of 30-50% by the end of the 21st Century (Table 3).Two scenarios, driven by HadGEM2.6 and HadCM3Q0 project a decline in August precipitation (Figure 7; Table 3).In all simulations, changes in annual total precipitation at the end of the century do not exceed ±5%.

The Projected Climate and Glacier Change
It should be noted that the ECHAM5-driven simulation projects a decline in summer precipitation in the 2006-2025 period in comparison with the baseline period: 22% in July and 14% in August.This is followed by an increase in precipitation towards the end of the century (Figure 7).
Figure 8 shows future projections of area and volume of the Tuyuksu glacier based on the developed climate scenarios (Figures 6 and 7).In all simulations, glacier area and volume decline by the middle of the 21st Century.Projections based on the HadGEM 2.6 and ECHAM5-driven simulations exhibit a faster change, particularly in glacier volume, possibly because of the projected reduction in precipitation in the first part of the 21st (Figure 7).Under these scenarios, however, changes in area and     Future changes in precipitation rates, averaged over the 25-year time slices, remain within the limits of interannual variability in precipitation measured during the baseline period and in the modelled baseline precipitation.Only two scenarios, driven by HadGEM8.5 and ECHAM5, projecting the strongest warming overall and a strong winter warming, respectively, project statistically significant changes in cold season precipitation with an increase of 30-50% by the end of the 21st Century (Table 3).Two scenarios, driven by HadGEM2.6 and HadCM3Q0 project a decline in August precipitation (Figure 7; Table 3).In all simulations, changes in annual total precipitation at the end of the century do not exceed ±5%.
It should be noted that the ECHAM5-driven simulation projects a decline in summer precipitation in the 2006-2025 period in comparison with the baseline period: 22% in July and 14% in August.This is followed by an increase in precipitation towards the end of the century (Figure 7).
Figure 8 shows future projections of area and volume of the Tuyuksu glacier based on the developed climate scenarios (Figures 6 and 7).In all simulations, glacier area and volume decline by the middle of the 21st Century.Projections based on the HadGEM 2.6 and ECHAM5-driven simulations exhibit a faster change, particularly in glacier volume, possibly because of the projected reduction in precipitation in the first part of the 21st (Figure 7).Under these scenarios, however, changes in area and volume of the glacier stabilizes by the middle of the century in line with a slow increase in summer temperatures after the 2050s.The HadGEM2.6-driven scenario projects the lowest volume and area loss of 38% and 34% respectively by the end of the century.Under the HadCM3Q0 and particularly HadGEM8.5-drivenscenarios, wastage of the Tuyuksu glacier continues until the end of the century despite a statistically significant increase in cold-season precipitation projected by the HadGEM8.5-drivenscenario (Table 3).This scenario projects the strongest loss of volume and area of 50% and 39% respectively by 2095.Volume of the glacier stabilizes by the middle of the century in line with a slow increase in summer temperatures after the 2050s.The HadGEM2.6-driven scenario projects the lowest volume and area loss of 38% and 34% respectively by the end of the century.Under the HadCM3Q0 and particularly HadGEM8.5-drivenscenarios, wastage of the Tuyuksu glacier continues until the end of the century despite a statistically significant increase in cold-season precipitation projected by the HadGEM8.5-drivenscenario (Table 3).This scenario projects the strongest loss of volume and area of 50% and 39% respectively by 2095.
Neither scenario predicts the complete disappearance of the Tuyuksu glacier by the end of the 21 st Century.

Streamflow Model Calibration and Validation: Model Efficiency
The HBV-ETH model showed robust performance in both calibration and validation periods as shown by the model efficiency statistics calculated using daily streamflow data; comparison of mean, high (Q10) and low (Q90) flow values (Table 4) and comparison of the observed and simulated hydrographs (Figure 9).In both periods, NSE and R 2 values, characterizing daily variations in the observed and modelled flow, were within 0.84-0.93 and 0.84-0.95ranges, respectively, and mean error was around 5% indicating good performance of the model in comparison with similar models applied across a range of high-elevation catchments worldwide and in Central Asia [17,56,63].The highest model efficiency has been achieved for the catchments with the highest glacierization, i.e., the Kishi Almaty, Ulken Almaty and Talgar (Tables 1 and 4).In the Kaskelen catchment, which has the lowest glacierization (Table 1), the model efficiency was the lowest with NSE values of 0.82 and 0.73-0.75 in the in the calibration and validation periods, respectively, and the mean error of 10-11% (Table 4).During the 2000-2013 period, the model showed a slightly better fit in comparison with the 1985-1998 period in both, calibration and validation years.We attribute this improvement to a better quality of the observational data in the 2000s when restoration of the gauging network took place in Kazakhstan [2] and availability of longer sequences of uninterrupted data for model calibration.
Comparison of the simulated and observed mean monthly flow for the validation years shows the best model performance in the Kishi Almaty catchment, where both Tuyuksu glacier and station are located, and in the most heavily glacierized Talgar catchment.In the Ulken Almaty catchment, the observed flow exceeds the modelled flow in August on average by 12%.This difference is substantially lower in the 1980s-1990s, when it is limited to 5-7% in the validation years, increasing Neither scenario predicts the complete disappearance of the Tuyuksu glacier by the end of the 21st Century.

Streamflow Model Calibration and Validation: Model Efficiency
The HBV-ETH model showed robust performance in both calibration and validation periods as shown by the model efficiency statistics calculated using daily streamflow data; comparison of mean, high (Q10) and low (Q90) flow values (Table 4) and comparison of the observed and simulated hydrographs (Figure 9).In both periods, NSE and R 2 values, characterizing daily variations in the observed and modelled flow, were within 0.84-0.93 and 0.84-0.95ranges, respectively, and mean error was around 5% indicating good performance of the model in comparison with similar models applied across a range of high-elevation catchments worldwide and in Central Asia [17,56,63].The highest model efficiency has been achieved for the catchments with the highest glacierization, i.e., the Kishi Almaty, Ulken Almaty and Talgar (Tables 1 and 4).In the Kaskelen catchment, which has the lowest glacierization (Table 1), the model efficiency was the lowest with NSE values of 0.82 and 0.73-0.75 in the in the calibration and validation periods, respectively, and the mean error of 10-11% (Table 4).During the 2000-2013 period, the model showed a slightly better fit in comparison with the 1985-1998 period in both, calibration and validation years.We attribute this improvement to a better quality of the observational data in the 2000s when restoration of the gauging network took place in Kazakhstan [2] and availability of longer sequences of uninterrupted data for model calibration.registered in the Ulken Almaty [2].In these two anomalous years, used for validation, model efficiency declined with NSE and R 2 values decreasing to 0.74-0.77.However, in other post-2008 validation years, model efficiency remained high with NSE and R 2 values varying between 0.84 and 0.92.A comparison of the daily records of simulated and measured streamflow showed that in all these years, the model simulated the occurrence of individual peaks in streamflow while systematically underestimating streamflow values in the late July-August.
Figure 9. Example of the observed and simulated hydrographs using temperature and precipitation data from the Tuyuksu station for the validation years.Note that different scales are used for (i) Kishi Almaty, Ulken Almaty, Kaskelen and (ii) Talgar, Turgen.

Model Parameters and Additional Objective Functions
The achieved model parameters are shown in Table 5.The values are in good agreement with the values obtained for other glacierized catchments of Central Asia [9,27,30,32].Several meteorological stations, located at different elevations between 2500 m a.s.l. and 3500 m a.s.l.provide data for the Ulken Almaty and Kishi Almaty catchments (Table 2, Figure 1) enabling evaluation of model parameters, e.g., precipitation gradients.Thus precipitation gradients between the Bolshoe Almatinskoe Lake station, located at the bottom of the Ulken Almaty catchment at 2500 m a.s.l. and the Tuyuksu station, were 2.3% per 100 m and 1.9% per 100 m in the 1985-1999 and 2000-2013 periods, respectively, while modelled gradients were within 3-2.4% per 100 m range (Table 5  Comparison of the simulated and observed mean monthly flow for the validation years shows the best model performance in the Kishi Almaty catchment, where both Tuyuksu glacier and station are located, and in the most heavily glacierized Talgar catchment.In the Ulken Almaty catchment, the observed flow exceeds the modelled flow in August on average by 12%.This difference is substantially lower in the 1980s-1990s, when it is limited to 5-7% in the validation years, increasing after the 2007-2008 hydrological year when the highest on record JJA air temperature was registered in the region and the seasonal temperature anomaly exceeded the 1974-2015 mean, derived from the Tuyuksu record, by 1.5 • C. In the summer of 2010, the highest on record streamflow values were registered in the Ulken Almaty [2].In these two anomalous years, used for validation, model efficiency declined with NSE and R 2 values decreasing to 0.74-0.77.However, in other post-2008 validation years, model efficiency remained high with NSE and R 2 values varying between 0.84 and 0.92.A comparison of the daily records of simulated and measured streamflow showed that in all these years, the model simulated the occurrence of individual peaks in streamflow while systematically underestimating streamflow values in the late July-August.

Model Parameters and Additional Objective Functions
The achieved model parameters are shown in Table 5.The values are in good agreement with the values obtained for other glacierized catchments of Central Asia [9,27,30,32].Several meteorological stations, located at different elevations between 2500 m a.s.l. and 3500 m a.s.l.provide data for the Ulken Almaty and Kishi Almaty catchments (Table 2, Figure 1) enabling evaluation of model parameters, e.g., precipitation gradients.Thus precipitation gradients between the Bolshoe Almatinskoe Lake station, located at the bottom of the Ulken Almaty catchment at 2500 m a.s.l. and the Tuyuksu station, were 2.3% per 100 m and 1.9% per 100 m in the 1985-1999 and 2000-2013 periods, respectively, while modelled gradients were within 3-2.4% per 100 m range (Table 5).
Precipitation correction factor values were above 1.0 in the Ulken Almaty, Kishi Almaty, and Turgen catchments, indicating that basin precipitation is higher than precipitation received by the Tuyuksu station where measurements might experience a systematic underestimation of precipitation mainly due to wind and errors related to evaporation.In the Talgar and Kaskelen catchments, precipitation correction factor values were close to 1 (Table 1).
To test the model performance further, a comparison between the observed and simulated snow properties was performed.The model successfully reproduced individual snowfalls and timing of snow melt in the study catchments and generated realistic values of SWE (Figure 10).Comparison of the monthly means of SWE simulated for the Kishi Almaty catchment and SWE measured at the Chimbulak station showed a particularly strong correlation with R 2 values, reaching 0.63 in the 2002-2012 time period that includes both calibration and validation years (Figure 10) and 82% for the 5-day SWE means for calibration years only.In other catchments, the R 2 values were lower within the 0.30-0.40range for both calibration and validation years.This disparity, however, can result from strong spatial variability in SWE and it can be expected that catchments, located further away from Chimbulak, exhibit different SWE values.Precipitation correction factor values were above 1.0 in the Ulken Almaty, Kishi Almaty, and Turgen catchments, indicating that basin precipitation is higher than precipitation received by the Tuyuksu station where measurements might experience a systematic underestimation of precipitation mainly due to wind and errors related to evaporation.In the Talgar and Kaskelen catchments, precipitation correction factor values were close to 1 (Table 1).
To test the model performance further, a comparison between the observed and simulated snow properties was performed.The model successfully reproduced individual snowfalls and timing of snow melt in the study catchments and generated realistic values of SWE (Figure 10).Comparison of the monthly means of SWE simulated for the Kishi Almaty catchment and SWE measured at the Chimbulak station showed a particularly strong correlation with R 2 values, reaching 0.63 in the 2002-2012 time period that includes both calibration and validation years (Figure 10) and 82% for the 5day SWE means for calibration years only.In other catchments, the R 2 values were lower within the 0.30-0.40range for both calibration and validation years.This disparity, however, can result from strong spatial variability in SWE and it can be expected that catchments, located further away from Chimbulak, exhibit different SWE values.

GLUE Analysis
GLUE analysis was applied to a range of randomly selected model runs which generated NSE values over 0.50 (Section 3.6) with varying model parameters for validation years between 2007 and 2013.The CR values were above 75% indicating that in the overall majority of simulations, monthly values of the observed flow were contained within the 5th to 95th percentile envelope of monthly flow values simulated by the model.Figure 11 shows results for the Talgar catchment as an example.In the Talgar and Kaskelen catchments, the highest uncertainty characterizes winter flow with monthly values of the observed flow being closer to the 5th (high flow) percentile of the simulated flow or outside the predicted envelope.This was a persistent overestimation which became evident when averaged over five hydrological years as well as in individual months in the Talgar (Figure 11).

GLUE Analysis
GLUE analysis was applied to a range of randomly selected model runs which generated NSE values over 0.50 (Section 3.6) with varying model parameters for validation years between 2007 and 2013.The CR values were above 75% indicating that in the overall majority of simulations, monthly values of the observed flow were contained within the 5th to 95th percentile envelope of monthly flow values simulated by the model.Figure 11 shows results for the Talgar catchment as an example.In the Talgar and Kaskelen catchments, the highest uncertainty characterizes winter flow with monthly values of the observed flow being closer to the 5th (high flow) percentile of the simulated flow or outside the predicted envelope.This was a persistent overestimation which became evident when averaged over five hydrological years as well as in individual months in the Talgar (Figure 11).During the warm season, the observed flow was close to the 50th percentile of the simulated flow values in the Talgar (Figure 11) as well as other catchments except the Ulken Almaty and Kaskelen where simulated flow was below the 50th percentile.In the Ulken Almaty, this underestimation probably reflects non-stationarity in the observed summer flow which was increasing in comparison with the pre-2006 calibration years [2].

Streamflow Scenarios
Streamflow was modelled using PRECIS data for the baseline period of 1986-2005 and scenarios were generated for the 25-year time slices with future projections starting in 2006 using four different climate and deglacierization scenarios (Section 4.1).
Future changes in streamflow were estimated relative to the streamflow during the baseline period simulated using the modelled climate data from PRECIS experiments.This is a standard practice of analyzing climate and hydrological projections [27,30,56] ensuring consistency between the baseline period and future scenarios and avoiding uncertainty due to the discrepancy between modelled and observed data.
In the context of water management, not only magnitude of the projected change is important but also how large it is compared to natural variability.Systems are adapted to climate variability and it is important to focus on changes that exceed the observed variability as they can trigger unprecedented impacts [64].Here, we compare the projected changes in discharge to variability (characterized by standard deviation) in monthly streamflow measured during the baseline period (Figures 12-16).The projected changes in discharge, which exceed standard deviations in the observed flow, are classed as significant.We acknowledge that variability in the observed data is different from internal variability in the modelled data and stress that for practical water management, comparison with observational data is important.

Streamflow Scenarios
Streamflow was modelled using PRECIS data for the baseline period of 1986-2005 and scenarios were generated for the 25-year time slices with future projections starting in 2006 using four different climate and deglacierization scenarios (Section 4.1).
Future changes in streamflow were estimated relative to the streamflow during the baseline period simulated using the modelled climate data from PRECIS experiments.This is a standard practice of analyzing climate and hydrological projections [27,30,56] ensuring consistency between the baseline period and future scenarios and avoiding uncertainty due to the discrepancy between modelled and observed data.
In the context of water management, not only magnitude of the projected change is important but also how large it is compared to natural variability.Systems are adapted to climate variability and it is important to focus on changes that exceed the observed variability as they can trigger unprecedented impacts [64].Here, we compare the projected changes in discharge to variability (characterized by standard deviation) in monthly streamflow measured during the baseline period (Figures 12-16).The projected changes in discharge, which exceed standard deviations in the observed flow, are classed as significant.We acknowledge that variability in the observed data is different from internal variability in the modelled data and stress that for practical water management, comparison with observational data is important.), HBV-ETH reliably simulates the annual cycle of streamflow in all catchments and for all climate scenarios.Glacier areas derived from the satellite imagery were used in this period and, therefore, uncertainties of the baseline hydrological simulations are mainly associated with model calibration and, more importantly, the quality of climate model output which serves as input for HBV-ETH.Thus in the Kishi Almaty (Figure 12), streamflow simulations based on RCM output is very close to the observed flow and this is consistent with high model efficiency during the calibration and validation years.In the Ulken Almaty and the Talgar catchments, an underestimation of flow in August, a month dominated by the glacier and ground-ice melt, is evident (Figures 12 and 13) similarly to the calibration and validation years (Figure 9).), HBV-ETH reliably simulates the annual cycle of streamflow in all catchments and for all climate scenarios.Glacier areas derived from the satellite imagery were used in this period and, therefore, uncertainties of the baseline hydrological simulations are mainly associated with model calibration and, more importantly, the quality of climate model output which serves as input for HBV-ETH.Thus in the Kishi Almaty (Figure 12), streamflow simulations based on RCM output is very close to the observed flow and this is consistent with high model efficiency during the calibration and validation years.In the Ulken Almaty and the Talgar catchments, an underestimation of flow in August, a month dominated by the glacier and ground-ice melt, is evident (Figures 12 and 13) similarly to the calibration and validation years (Figure 9).Despite bias correction, the use of PRECIS outputs contributes to the overall uncertainty.For example, application of the data from the PRECIS runs, driven by the 'cold and wet' ECHAM5 model, result in a stronger underestimation of streamflow in August, particularly in the catchments with higher elevations and higher glacierization, e.g., the Ulken Almaty (Figure 12d) and the Talgar (Figure 14d).In this case, the use of the 'cold and wet' model enhances uncertainty of calibration (Section 4.2.1;Section 4.2.3).By contrast, simulations using the RCM outputs driven by the 'warm and dry' HadCM3Q0 model result in higher flows in July-August, which better matches observed flows, but overestimates runoff in June (Figures 12c and 14c).
The largest discrepancies between the observed flow and flow simulated for the baseline period using PRECIS data are found in the Kaskelen and Turgen catchments.Both, model efficiency and performance of the objective functions (Section 4.2.2) were worse than in these catchments during the calibration and validation years when the model was forced with meteorological observations (Table 4) than in those with higher glacierization.In the Turgen, flow simulations based on all PRECIS outputs (ii) In the Turgen and Kaskelen catchments, the strongest reduction is projected under the most aggressive HadGEM 8.5-driven scenario.Under the HadGEM-driven scenarios (both RCP 2.6 and RCP 8.5), in both catchments summer flow declines already in the 2006-2025 period with small changes in the later years projected under the RCP 2.6 scenario.Under the A1B scenario, significant reduction in summer discharge is projected for the second half of the century.(iii) In the Kishi Almty catchment (the smallest of all by area but with the highest specific discharge), the strongest reduction in August-September and July-September flow is simulated under the scenarios projecting the weakest warming, i.e., HadGEM2.6 and ECHAM5-driven simulations, respectively (Table 5).Reduction in the August flow, dominated by the glacier melt, is projected already for the 2005-2025 period.The long-term observational discharge record extending to 2013 confirms this trend [2].(iv) In the catchments with higher glacierization, the Ulken Almaty and the Talgar, no significant change in flow is expected in summer except a reduction in August flow under the HadCMQ0-driven scenario.An increase in summer flow between 1980-2005 and 2006-2025 is projected under the HadGEM2.6 and HadCM3-driven scenarios and this is in agreement with the long-term observations extending to 2018 [2].(v) In catchments with higher glacierization, e.g., the Kishi Almaty, Ulken Almaty and Talgar, the smallest reduction in flow in summer and the largest increase in the cold season are projected under the most aggressive HadGEM 8.5 scenario (Table 5).Under this scenario, glacier melt is more intense and glacier retreat continues to the end of the 21st Century (Figure 8) and, therefore, greater specific discharge is generated from smaller glacier area; (vi) Displacement of hydrographs towards higher flow in spring-early summer (particularly May-June), consistent with earlier and more intensive snow melt, is evident in a number of simulations for the Ulken Almaty, Talgar, Turgen, and Kaskelen catchments.(vii) Increase in the cold season flow is projected under most catchment-scenario combinations except the least aggressive ECHAM5-driven scenario for the Kishi Almaty, Turgen, and Kaskelen.The projected increase is consistent with the observed long-term trends [2].(viii) Peak flow is projected already for the first quarter of the 21st Century in catchments with lower glacierization, i.e., the Kaskelen and Turgen, as well as for the Kishi Almaty.

Discussion
This paper presented projections of river discharge for five catchments in the northern Tien Shan with different specific discharge and glacierization varying between 2.6% in the Kaskelen catchment to 16% in the Ulken Almaty and Talgar catchments for four climate and glacier scenario-model combinations.Investigating a variety of catchments is important to assess the difference in response of discharge to the projected climate change and deglacierization.We stress that our assessment refers to the natural headwater catchments with glacio-nival nourishment and characterizes changes in input of water into the river system before any water abstraction or catchment management occurs.
The HBV-ETH model showed a robust performance in the calibration and validation years.Model efficiency was higher in the catchments with higher glacierization probably because melt is driven by temperature, which is spatially more homogenous than precipitation and therefore can be simulated more accurately by a simple degree-day approach.A strong elevation-dependent change in environmental conditions characterizes larger catchments, extending into the foothills (e.g., the Turgen and Kaskelen) and distributed models (e.g., [65]) may be more suitable than a lumped model.However, in all catchments simulated mean, high (Q10) and low (Q90) flows were close to the observed values (Table 4).Model parameters (Table 5) were within realistic ranges previously defined for the glacierized alpine catchments [17] and in agreement with observational data from the Kishi Almaty catchment.The SWE measurements were used as additional criteria to reduce parameter uncertainty.A good agreement was achieved between the observed and modelled data in the Kishi Almaty catchment, where measurements were made (R 2 of 0.63 for monthly mean in calibration and validation years and 0.82 for 5-day means in calibration years), but only a moderate agreement (R 2 of 0.3-0.4 for monthly means in calibration and validation years) in the Turgen and Kaskelen catchments.In the Ile Alatau, SWE is known to vary strongly with elevation, aspect and prevailing vegetation and the achieved moderate correlation with the modelled data can be due to the observed spatial variability in SWE [66].
GLUE analysis generated CR values in excess of 75% indicating a good agreement between the observed and modelled data [61,62].Observed flow values fall outside the modelled Q05 and Q95 range, being close to or exceeding Q05, predominantly in winter months (Figure 11).During this season, when flow values are low, uncertainty in flow measurements increases [2] which may be one reason for the discrepancy between the observed and modelled data.
In summer, the observed flow values were well within the simulated Q05-Q95 range.However, uncertainty in the simulation of late July-August flow in the Ulken Almaty catchment increases in the late 2000s in comparison with the 1980s-1990s in contrast to other catchments.We attribute this discrepancy to the increasing melt of ground ice, i.e., permafrost and rock glaciers which are particularly widespread in this catchment.The area occupied by rock glaciers constitutes about 30% of the glacierized area [67].Their melt intensified recently and our filed observations confirmed considerable discharge from rock glaciers in 2015-2019.Temperature of permafrost in this catchment increased by 0.3-0.6 • C between 1974 and 2009 [68] and there is widespread evidence for an increase in active layer thickness [67].Ground ice melt is not simulated by HBV-ETH because its input in streamflow was previously considered insignificant [27].We suggest that under the rapidly warming climate, it may be a source of systematic error in the simulation of the late summer flow as well of the underestimation of winter flow.An increase in the observed winter flow has been registered in all catchments except the Turgen since the 1990s [2] despite the sub-zero temperatures potentially due to the longer period of ground-ice melt and later freezing of soil.
For the catchments with high glacierization and significant contribution of glacier melt to total flow, the future evolution of the glacier cover is crucial.In contrast to previous studies, the developed hydrological projections were driven not only by three climate scenarios-model combinations and assumed glacier change, but also by deglacierization scenarios based on the combination of glacier dynamics and mass balance modelling.This approach provides data on transient evolution of glacier volume and area and enables an assessment of the timing of peak flow.
Previous projections of glacier change in High Asia mostly used parametrizations based on mass balance and did not explicitly include ice flow components (e.g., [9,[23][24][25][26]30,65]).In these studies, the projected trends showed a near linear decline in glacier volume and area to the end of the 21st Century.By contrast, our results generated by the combined HO 3-D ice dynamics and mass balance model show that the observed retreat of the Tuyuksu glacier group will continue until the middle of the 21st Century.By this time, a loss of about 40% of glacier volume and 35% of area is projected, followed by stabilization.These results, however, agree with projections developed by [24] for glaciers of northern Asia (e.g., the Altai) which is consistent with the location of Tuyuksu on the northernmost flank of Central Asia.Under the most aggressive HadGEM8.5-drivenscenario, glaciers continue retreating further, albeit at a slower rate losing 50% of their volume and 39% of their area by 2095 (Figure 8).The difference between changes in area and volume projected for different scenarios for the end of the century is relatively small.For the scenarios projecting the highest (HadGEM8.5)and the lowest (HadGEM2.6)loss, the difference is 18% for volume and 11% for area.The difference between HadGEM2.6,HadCMQ0 and ECHAM5-driven scenarios are 4-6% for area and 5-10% for volume.
A limitation of this approach is high data requirements which means that only changes in the Tuyuksu group of glaciers, which nourish the Kishi Almaty, could be modelled.Elevation-dependent rates of area change derived for the Tuyuksu group were applied to other glaciers.It was previously shown that changes in glacier area in the Ile and Kungey Alatau correlate closely with those of the Tuyuksu group [34].Over 90% of all glaciers in the study area have the same (northern) aspect.However, types and size vary and, therefore, on a scale of catchment glaciers can respond differently to climate change.Thus while in the Kishi Almaty catchment glacier area, represented largely by the Tuyuksu group, declined by 10.6% between 1990 and 2006, in the Kaskelen catchment it declined by 23% although in other catchments closer area reduction rates were observed (Table 1).In the absence of wider glaciological measurements required by HO 3-D model, this limitation has to be accepted as a source of uncertainty.
Another important limitation is that glaciological and climate models are not coupled and, therefore, elevation feedbacks on the surface mass balance have not been taken into account in future scenarios.Using degree-day factors, derived for the present climate in future simulations by HBV-ETH, introduces additional uncertainty to future scenarios.Degree-day factors are not time independent and can change substantially on both seasonal and decadal timescales due to varying feedbacks from energy balance components which are not determined by the temperature, e.g., albedo [69,70].Therefore the factors, derived for the present conditions, may break down in a future climate.
Our hydrological scenarios show that in catchments with low glacerization (Kaskelen, Turgen), total river discharge in July-August declines under all scenarios by the last quarter of the 21st Century in comparison with the baseline period.A particularly strong decrease of 34-37% occurs under the HadGEM8.5-drivenscenario (Table 6; Figures 15 and 16).Under other scenarios, the July-August discharge decreases by 9-30% and, with exception of the HadCM3Q0-driven scenario for the Kaskelen, these changes exceed variability in the observed and modelled discharge during the baseline period.By contrast, in the Ulken Almaty and Talgar catchments with high glacierization (16% as in 2006; Table 1), the July-August flow is not projected to change beyond the currently observed interannual variability (except HadCM3Q0-driven scenario for the Talgar) (Table 6; Figures 12 and 14).In the Kishi Almaty catchment, characterized by lower glacierzation of 12% and high specific discharge, the least aggressive scenarios (driven by HadGEM2.6 and ECHAM5 models) generate flow reduction of about 30% in August.Under the more aggressive HadGEM8.5-drivenscenarios, which project glacier retreat beyond the middle of the century, flow reduction is smaller at 2-19% and does not exceed variability in the observed flow during the baseline period (Table 6; Figure 13) because more the intense melt from a smaller glacier area compensates for the projected decline in summer precipitation (Figure 7) and increase in evaporation.These results are different from the earlier sensitivity studies (e.g., [27]) or discharge projections which used glacier parametrization based on mass balance trends and excluding ice flow component (e.g., [9,30]).All of these studies projected a stronger decline in river flow for more aggressive scenarios in catchments with high glacierization in Central Asia.Flow increase was projected in the Hindukush-Himalaya region.However, the latter is driven by the projected increase in precipitation rather than glacier melt [26,65].By contrast, changes in annual precipitation in the northern Tien Shan at the end of the century are limited to ±5% in comparison with the baseline period and are not statistically significant under any scenario.Therefore, in the northern Tien Shan, in catchments with high glacierization, the most aggressive climate scenarios are expected to cause least stress to water resources in the 21st Century.This, however, applies to the upstream natural catchments without water abstraction and with relatively low evaporation.
An increase in May-June flow, dominated by snow melt, and/or displacement of seasonal high flow from July-August to May-June are projected for most catchment-scenario combinations (Table 6; Figures 12-16).In particular, in the high-elevation Talgar catchment, spring flow is projected to increase by the factor of 1.5-2.5 under the HadGEM8.5-drivenscenario whereby snow melt occurs at earlier stages and at higher elevations.
The timing of peak flow is important with regard to water resources and development of adaptation strategies and practice.A decline in July-August flow is projected already for the 2005-2025 time slice in the Kaskelen and Turgen under both HadGEM scenarios.In the Kishi Almaty, August flow is expected to decline under all scenarios although the projected reduction is within the observed interannual variability.Flow measurements in the Kishi Almaty and Kaskelen catchments show that July-August flow was marginally lower in 2001-2013 in comparison with two previous decades while there has been no change in summer flow in the Turgen catchment between 1990s and 2000s [2].In the Ulken Almaty catchment, an increase in late spring and/or summer flow is projected for 2006-2025 under the Hadley Centre model scenarios with subsequent return to the baseline values in 2026-2050.The projected increase is in agreement with observations extending to 2018.Overall, changes in the observed flow [2] and simulations of future changes in river discharge suggest that in the northern Tien Shan summer flow peaks in the first quarter of the 21st Century.This conclusion is important for the development of adaptation strategies in both water and agricultural sectors.To substantiate it further, a wider range of scenario-model combinations is required given a large spread of climate projections which impact glaciological and hydrological scenarios [24].

Conclusions
The Tien Shan Mountains are water towers of Central Asia and the state of runoff generated in the mountains is of upmost importance to national economies and international relations in the region.South-eastern Kazakhstan-a region accommodating large cities, industrial centres, and agriculturerelies on discharge from the melting snow and glaciers of the Ile Alatau.Our results suggest that peak flow in the region occurs in the first quarter of the 21st Century.Will the regional water towers run dry in the future?
Our study showed that under the warming climate and insignificant changes in precipitation, glaciers will lose half of their volume and third of their area by the end of the 21st Century but will not disappear completely.Impact of the projected climate change on river discharge will depend on glacierization of catchments.Those with low (currently 2-4%) glacierization (Turgen, Kaskelen) will be vulnerable as 20-37% reduction in summer flow is expected particularly under the more aggressive scenarios.The projected decline may have significant impact on water availability for summer irrigation which dominates water consumption in these catchments.In catchments with high glacierization (currently 16% and over; Ulken Almaty, Talgar), no statistically significant changes in summer flow are expected while spring flow increases under the aggressive scenarios.The Ulken Almaty supplies water for municipal use in Almaty city and is regulated via a reservoir.In the Talgar catchment, the projected increase in spring flow may require development of flood prevention measures.In catchments with medium glacierization (10-12%, Kishi Almaty), reduction in summer flow is projected for the less aggressive scenarios while under the strongest warming, more intensive glacier and snow melt will sustain water supply.

Figure 1 .
Figure 1.Study area.Letters refer to the location of meteorological stations: A-Chimbulak, B-Mynzhilki, C-Tuyuksu, D-Bolshoe Almatinskoe Lake.Flow gauging sites are numbered as in Table1.

Figure 2 .
Figure 2. Observed and simulated temperature (a, b) and precipitation rates (c, d) before (a, c) and after (b, d) bias correction using the Empirical Quantile Mapping (EQM) and Delta methods for temperature and precipitation respectively for the baseline period.

Figure 3 .
Figure 3.Comparison of bias correction methods applied to the simulations of air temperature for the baseline period.PRECIS simulations driven by ECHAM5 are used as example.

Figure 2 .
Figure 2. Observed and simulated temperature (a,b) and precipitation rates (c,d) before (a,c) and after (b,d) bias correction using the Empirical Quantile Mapping (EQM) and Delta methods for temperature and precipitation respectively for the baseline period.

Figure 2 .
Figure 2. Observed and simulated temperature (a, b) and precipitation rates (c, d) before (a, c) and after (b, d) bias correction using the Empirical Quantile Mapping (EQM) and Delta methods for temperature and precipitation respectively for the baseline period.

Figure 3 .
Figure 3.Comparison of bias correction methods applied to the simulations of air temperature for the baseline period.PRECIS simulations driven by ECHAM5 are used as example.

Figure 3 .
Figure 3.Comparison of bias correction methods applied to the simulations of air temperature for the baseline period.PRECIS simulations driven by ECHAM5 are used as example.

Water 2020, 12 , 627 10 of 32 Figure 5 .
Figure 5. Schematic structure of HBV-ETH model (after[17]).Input data include daily temperature and precipitation as well as catchment topography and glacier cover statistics.The model has snow and glacier modules and uses seasonally-varying degree day factors.The upper panel shows the main model parameters.

Figure 5 .
Figure 5. Schematic structure of HBV-ETH model (after [17]).Input data include daily temperature and precipitation as well as catchment topography and glacier cover statistics.The model has snow and glacier modules and uses seasonally-varying degree day factors.The upper panel shows the main model parameters.

Figures 6
Figures 6 and 7 and Table3show future air temperature and precipitation scenarios for the study region.All scenarios show a statistically significant increase in air temperature.The HadCM3Q0-driven A1B scenario and HadGEM RCP 8.5 show the strongest warming, the former in summer and the latter across all seasons with summer temperatures increasing by 6-7 • C and 7-8 • C respectively (Table3).The HadGEM-driven RCP 2.6 scenario shows a significant warming in the first quarter of the 21st Century in comparison to the baseline period in all seasons followed by a smaller increase in summer temperature.A brief analysis of the non-downscaled CMIP5 simulations showed that most

Figure 6 .
Figure 6.Air temperature scenarios for the study region.Note an increase in air temperature in the first quarter of the century under HadGEM 2.6 scenario and a continuing increase in temperature under HadGEM 8.5 scenario.

Figure 6 . 32 Figure 6 .
Figure 6.Air temperature scenarios for the study region.Note an increase in air temperature in the first quarter of the century under HadGEM 2.6 scenario and a continuing increase in temperature under HadGEM 8.5 scenario.

Figure 7 .Figure 7 .
Figure 7. Precipitation scenarios for the study region.Note that neither model projects significant changes in precipitation in spring-early summer when precipitation maximum is observed.A shifts in seasonal maximum towards spring is projected under the HadGEM8.5 scenario.

Water 2020, 12 , 627 15 of 32 Figure 8 .
Figure 8. Projected changes in volume and area of Tuyuksu glacier.Glacier stabilization is projected at the middle of the century.Under the most aggressive HadGEM 8.5 scenario glacier wastage continues until the end of the century.

Figure 8 .
Figure 8. Projected changes in volume and area of Tuyuksu glacier.Glacier stabilization is projected at the middle of the century.Under the most aggressive HadGEM 8.5 scenario glacier wastage continues until the end of the century.

Figure 9 .
Figure 9. Example of the observed and simulated hydrographs using temperature and precipitation data from the Tuyuksu station for the validation years.Note that different scales are used for (i) Kishi Almaty, Ulken Almaty, Kaskelen and (ii) Talgar, Turgen.

Figure 10 .
Figure 10.Comparison of monthly mean values of snow water equivalent measured at the Chimbulak meteorological station and simulated by HBV-ETH for the Kishi Almaty catchment for both calibration and validation years in 2002-2012.

Figure 10 .
Figure 10.Comparison of monthly mean values of snow water equivalent measured at the Chimbulak meteorological station and simulated by HBV-ETH for the Kishi Almaty catchment for both calibration and validation years in 2002-2012.

32 Figure 11 .
Figure 11.Envelope of simulated values ranging between 5th, 50th, and 95th percentiles containing the observed and simulated mean monthly flow values for the Talgar catchment.Data are averaged over the 2007-2013 period.Validation years are used.

Figure 11 .
Figure 11.Envelope of simulated values ranging between 5th, 50th, and 95th percentiles containing the observed and simulated mean monthly flow values for the Talgar catchment.Data are averaged over the 2007-2013 period.Validation years are used.

Figure 12 .
Figure 12.Streamflow scenarios for the Ulken Almaty based on the PRECIS outputs for the baseline period (1981-2005; glacierization for 1990 and 2006) and future 25-year time slices.Grey bars show observed flow values and their standard deviations (±1 σ shown by vertical bars).

Figure 12 . 32 Figure 12 .
Figure 12.Streamflow scenarios for the Ulken Almaty based on the PRECIS outputs for the baseline period (1981-2005; glacierization for 1990 and 2006) and future 25-year time slices.Grey bars show observed flow values and their standard deviations (±1 σ shown by vertical bars).

Figure 13 .
Figure 13.As in Figure 12 but for the Kishi Almaty.

Figure 13 .
Figure 13.As in Figure12but for the Kishi Almaty.

Figure 14 .
Figure 14.As in Figure 12 but for the Talgar.Note that observational data refer to the 1983-1995 period.

Figure 15 .
Figure 15.As in Figure 12 but for the Kaskelen.

Figure 16 .
Figure 16.As in Figure 12 but for the Turgen.4.3.1.Simulation of Flow Using Biased-Corrected PRECIS Data for the Baseline Period In the baseline period (1981-2005, grey bars and black line in Figures12-16), HBV-ETH reliably simulates the annual cycle of streamflow in all catchments and for all climate scenarios.Glacier areas derived from the satellite imagery were used in this period and, therefore, uncertainties of the baseline hydrological simulations are mainly associated with model calibration and, more importantly, the quality of climate model output which serves as input for HBV-ETH.Thus in the Kishi Almaty (Figure12), streamflow simulations based on RCM output is very close to the observed flow and this is consistent with high model efficiency during the calibration and validation years.In the Ulken Almaty and the Talgar catchments, an underestimation of flow in August, a month dominated by the glacier and ground-ice melt, is evident (Figures12 and 13) similarly to the calibration and validation years (Figure9).

Figure 16 .
Figure 16.As in Figure 12 but for the Turgen.4.3.1.Simulation of Flow Using Biased-Corrected PRECIS Data for the Baseline PeriodIn the baseline period, grey bars and black line in Figures12-16), HBV-ETH reliably simulates the annual cycle of streamflow in all catchments and for all climate scenarios.Glacier areas derived from the satellite imagery were used in this period and, therefore, uncertainties of the baseline hydrological simulations are mainly associated with model calibration and, more importantly, the quality of climate model output which serves as input for HBV-ETH.Thus in the Kishi Almaty (Figure12), streamflow simulations based on RCM output is very close to the observed flow and this is consistent with high model efficiency during the calibration and validation years.In the Ulken Almaty and the Talgar catchments, an underestimation of flow in August, a month dominated by the glacier and ground-ice melt, is evident (Figures12 and 13) similarly to the calibration and validation years (Figure9).Despite bias correction, the use of PRECIS outputs contributes to the overall uncertainty.For example, application of the data from the PRECIS runs, driven by the 'cold and wet' ECHAM5 model, result in a stronger underestimation of streamflow in August, particularly in the catchments with higher elevations and higher glacierization, e.g., the Ulken Almaty (Figure12d) and the Talgar (Figure14d).In this case, the use of the 'cold and wet' model enhances uncertainty of calibration (Section 4.2.1;Section 4.2.3).By contrast, simulations using the RCM outputs driven by the 'warm and dry' HadCM3Q0 model result in higher flows in July-August, which better matches observed flows, but overestimates runoff in June (Figures12c and 14c).The largest discrepancies between the observed flow and flow simulated for the baseline period using PRECIS data are found in the Kaskelen and Turgen catchments.Both, model efficiency and performance of the objective functions (Section 4.2.2) were worse than in these catchments during the calibration and validation years when the model was forced with meteorological observations (Table4) than in those with higher glacierization.In the Turgen, flow simulations based on all PRECIS outputs

Table 1 .
Characteristics of the study catchments.Minimum elevation corresponds to the elevation of the flow gauge which provided data used in this study.

Table 3 .
Projected changes in monthly air temperature and precipitation in 2076-2095 in comparison with the baseline period(1985/1986-2004/2005).Absolute changes are shown for temperature and relative for precipitation.Changes in precipitation significant at 5% confidence level are highlighted in bold.

Table 4 .
Statistics of river runoff and model efficiency during model calibration and validation with 1990 and 2006 glacier masks.'Year' refers to the year for which glacier mask was derived.'O' stands for observed; 'S' for simulated values.The 1985-1995 and 2001-2005 (2007-2010 for the Talgar) periods are denoted as '1' and 2' respectively.NSE-Nash-Sutcliffe efficiency criterion; ME is mean error. ).
As in Figure12but for the Talgar.Note that observational data refer to the 1983-1995 period.Figure 13.As in Figure 12 but for the Kishi Almaty.Figure 14.As in Figure 12 but for the Talgar.Note that observational data refer to the 1983-1995 period.