Comparison of climate change signals in CMIP3 and CMIP5 multi-model ensembles and implications for Central Asian glaciers

. Central Asian water resources largely depend on melt water generated in the Pamir and Tien Shan mountain ranges. To estimate future water availability in this region, it is necessary to use climate projections to estimate the future glacier extent and volume. In this study, we evaluate the impact of uncertainty in climate change projections on the future glacier extent in the Amu and Syr Darya river basins. To this end we use the latest climate change projections generated for the upcoming IPCC report (CMIP5) and, for comparison, projections used in the fourth IPCC assessment (CMIP3). With these projections we force a regional-ized glacier mass balance model, and estimate changes in the basins’ glacier extent as a function of the glacier size distribution in the basins and projected temperature and precipitation. This


Introduction
The fate of Asian glaciers under climate change has been the topic of a heated scientific debate (Cogley et al., 2010;Immerzeel et al., 2010;Kargel et al., 2011;Bolch et al., 2012;Sorg et al., 2012). A main reason for this is the lack of systematic cryospheric observations and the absence of robust methods that can assess glacier evolution under climate change at the large river basin scale (Unger-Shayesteh et al., 2013). Downstream water availability in several large Asian rivers is highly sensitive to changes in snow and glacier extent (Immerzeel and Bierkens, 2012), and large populations depend on the water generated upstream. This dependence is likely to increase as irrigated areas further expand under population growth (Wada et al., 2011).
To assess future changes in high mountain hydrology, glacio-hydrological models forced by climate scenarios are used. Future glacier extent is a combined result of the glacier mass balance and ice-flow dynamics. While mass balance modeling is rather straightforward to implement and approaches of different complexity can be used (from simple degree-day to energy-balance models for calculation of ablation), changes in glacier geometry due to ice flow are more A. F. Lutz et al.: Climate change implications for Central Asian glaciers complex to include. At the same time, changes in glacier geometry have to be considered in regions where glacier melt makes a significant contribution to total runoff. Ideally, these should be simulated with mass balance models combined with two-or three-dimensional ice flow dynamics (e.g. Huss et al., 2007;Jouvet et al., 2008), but these are computationally demanding and require detailed knowledge of glacier bed geometry and ice thickness distribution. Other approaches have been developed in which ice is transported from the accumulation zone to the ablation zone through basal sliding or creep (e.g. Immerzeel et al., 2011Immerzeel et al., , 2013, but, like models of ice flow dynamics, this approach is only applicable for small catchments as it requires modeling at high spatial resolution. In several hydrological models, glaciers are treated as static entities that generate melt water and the glacier extent is modified for the future by making crude assumptions on the ice mass balance (e.g. Immerzeel et al., 2010) or by imposing hypothetical glacier scenarios (e.g. Singh and Bengtsson, 2004;Rees and Collins, 2006;Singh et al., 2006;Finger et al., 2012). A commonly used alternative method is to use volume-area scaling relationships (e.g. van de Wal and Wild, 2001;Möller and Schneider, 2010;Radić and Hock, 2011).
A parameterization of future glacier evolution has been developed for individual glacier systems (Huss et al., 2010). Although this approach can be applied to any area, it requires recalibration based on repeated digital elevation models (DEMs) for different glacier types. Several global scale models that simulate glacier mass balances have been developed (e.g. Hirabayashi et al., 2010;Radić and Hock, 2011), but limited approaches to assess glacier evolution at the large river basin scale are available. To our knowledge only few studies of glacier changes at basin scale have been conducted (Prasch, 2010;Weber et al., 2010;Prasch et al., 2013), all using the same modeling approach. This approach uses an energy-balance model for the calculation of melt and therefore requires additional atmospheric input besides air temperature. Thus, there is a strong need for an approach that can be applied at the large river basin scale, requires a minimum of data inputs which are readily available and which generalises changes in glacier extent over large areas without the need to model individual glaciers. At the same time this approach has to yield a reliable estimate of future glacier extent at the large river basin scale.
Models to estimate future ice areas and volumes are commonly forced by air temperature and precipitation provided by general circulation models (GCMs) which are downscaled to the study region. However, there is a large spread in the GCM projections Sutton, 2009, 2010;Radić and Clarke, 2011). This large spread is especially true for precipitation in Asia (Immerzeel et al., 2010). There is growing agreement that impact studies should be forced by an ensemble of GCMs outputs Sutton, 2009, 2010). While this has been done for North America (e.g. Radić and Clarke, 2011;Zhang et al., 2011), river basins originating in the European Alps (e.g. Huss, 2011;Farinotti et al., 2012), for river basins worldwide (e.g. Nohara et al., 2006), or for selected glaciers (e.g. Giesen and Oerlemans, 2010), no detailed assessments are available for Central Asia. Sutton (2009, 2010) identified three main sources of uncertainty in future climate projections: (i) model uncertainty due to the structural differences among GCMs, by which different models produce different projections for the same radiative forcing; (ii) scenarios uncertainty due to different radiative forcing; and (iii) uncertainty due to the natural climate variability. They showed that the first source of uncertainty is the larger throughout the century for both temperature and precipitation. It seems therefore imperative to take it into account in impact studies of glacier changes.
The aim of this study is to quantify the impact of uncertainty in climate change projections on the future glacier extent in the Amu Darya and Syr Darya river basins; two meltwater influenced rivers which provide the most important water sources in the Central Asian region. Therefore we analyse the differences in uncertainty range between the latest climate change projections provided by the fifth Coupled Model Intercomparison Project (CMIP5) generated for the upcoming fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC), and climate change projections used for the fourth IPCC assessment report (CMIP3). These projections for the climate from 2008 to 2050 are analysed at a monthly scale, and the results are used to force a glacier model simulating the future response of glaciers and changes in glacier geometry at the basin scale. We quantify the uncertainty in glacier projections as a result of the range in the climate change projections, and show how this uncertainty differs between the CMIP3 and CMIP5 ensembles. Moreover, the sensitivity of the presented approach to the model parameters is separately addressed, and the approach is validated.

Study area
The sources of the Amu Darya and Syr Darya rivers are located in the Pamir and Tien Shan mountains, respectively ( Fig. 1), and both rivers drain into the Aral Sea. Water allocation is a highly sensitive topic in the region. The upstream countries (Kyrgyzstan and Tajikistan) use water mainly for hydropower production during winter, whereas the downstream countries (Uzbekistan, Turkmenistan and Kazakhstan) utilise water for irrigation during summer where around 22 million people depend on irrigated agriculture (Siegfried et al., 2012). Glacier melt provides an important source of water in both basins, given the dry and warm climate downstream Sorg et al., 2012). The total glacierized area is 10 289 km 2 (1.3 % of total 799 261 km 2 basin area) in the Amu Darya basin and 1596 km 2 (0.14 % of total 1 117 625 km 2 basin area) in the Syr Darya basin, as calculated from the Randolph Glacier Inventory version 2.0 (Arendt et al., 2012), which for Central Upstream parts of the Amu and Syr Darya river basins (in green and pale blue, respectively), the main river system (blue lines), the initial glacierized fraction per 1 km grid cell (red shades) and political boundaries (black lines).

Digital elevation models
In this study two DEMs are used. Both are based on the Shuttle Radar Topographic Mission (SRTM) DEM at a nominal resolution of 90 m. For the downscaling of GCMs, this DEM is resampled to 1 km resolution. From here on, this DEM will be referred to as the 1 km DEM, and 1 km will also be the spatial resolution of the glacier model. For sub-grid calculations, the SRTM DEM at 90 m resolution is used. This DEM is referred to as the 90 m DEM.

Climate data
A dataset of precipitation and temperature spanning thirty years  is used as reference for the climate change assessment. For this period, we use the Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE, Yatagai et al., 2012) dataset for precipitation and Princeton's Global Meteorological Forcing Dataset (PGMFD, Sheffield et al., 2006) for temperature. APHRODITE is a long-term continental-scale daily precipitation product based on a dense network of rain gauges, with spatial resolution of 0.25 • (≈18-30 km in the studied area). The PGMFD was constructed by combining a suite of global observationbased datasets with the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis and it has a daily resolution and a spatial scale of 0.5 • (≈36-60 km in the studied area). Daily precipitation data are bilinearly interpolated to 1 km resolution from the APHRODITE 0.25 • gridded precipitation dataset grid cell centers. Gridded daily average near-surface air temperature data at 1 km resolution are obtained by bilinear interpolation from grid cell centers in the PGMFD 0.5 • gridded temperature dataset, which are subsequently corrected for elevation using the 1 km DEM and a vertical temperature lapse rate (Table 1).

Climate change projections
We use the set of global climate change simulations which is used as basis for the upcoming fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC), the CMIP5 multi-model ensemble (Taylor et al., 2012). All simulations which were available online in the PCMDI database (http://cmip-pcmdi.llnl.gov/cmip5/) earlier than 15 December 2011 are included in the analysis. In order to compare the CMIP5 multi-model ensemble to the previous generation of global climate change simulations, the CMIP3 multi-model ensemble (Meehl et al., 2007), which is the basis of the fourth IPCC assessment report, is also analysed.

Glaciers
Glacier covered areas in the Amu and Syr Darya river basins are extracted from the Randolph Glacier Inventory version 2.0 (RGI 2.0) dataset (Arendt et al., 2012). We updated the RGI 2.0 with more recently mapped glacier outlines provided by T. Bolch. The updates include outlines for the large glacier systems in the Fedchenko glacier region, which are not available in RGI 2.0 as well as more accurate outlines for numerous other glaciers in the Pamir and Tien Shan mountain ranges. We assume this compiled dataset of glacier extent to represent the glacier extent at the end of the reference period (2007), and to form the starting point for the future simulations of glacier extent.   (Immerzeel et al., 2012a), MB OBS is taken from (WGMS, 2011 From this dataset, the size distribution of glaciers is extracted (Fig. 2). In the Amu Darya and Syr Darya river basins combined, 50 % of the total glacier area consists of glaciers with a surface area smaller than 25 km 2 and 11 % of the glacier area consists of glaciers smaller than 1 km 2 . The median glacier size in the basin is 0.21 km 2 . From this distribution 26 different glacier size classes are defined and used for further analysis (Fig. 2).
The dataset with glacier extents is also used for the calculation of an initial fractional glacier cover per 1 km grid cell, to be used as starting point for the glacier model simulations. Each 1 km grid cell of the 1 km DEM is assigned a fractional glacier cover varying from 0 (no glacier cover) to 1 (entirely covered with glaciers) (Fig. 1).
For model calibration, the average of the observed annual mass balance in the region's mountains is used, which is approximately −0.47 m water equivalent (w.e.) per year between 1978 and 2007, based on five glaciers with mass balance records in the region (WGMS, 2011) ( Table 2).

Downscaling of GCM output
Downscaling of the GCMs outputs is necessary due to the large scale discrepancy between the climate models (operated on grids of 100 km grid distance or more) and the glacier model (operating on the 1 km scale). Since in our study, the major focus is on uncertainty stemming from the climate simulations, we include as many climate simulations as possible. We consider the CMIP3 and CMIP5 simulations based on all available emission scenarios: SRES B1, A1B, and A2 (Nakicenovic et al., 2000) in the case of CMIP3, and rcp2.6, rcp4.5, rcp6.0, and rcp8.5 (Meinshausen et al., 2011) in the case of CMIP5. Since it is difficult to associate probabilities to the emission scenarios, we do not use any prior assumption and give the same weight to all scenarios. We therefore calculate percentiles for all GCM realizations according to the inverse number of GCM realizations per scenario. We extract the grid cells of the climate models over the study region and analyse projected annual and monthly temperature and precipitation averaged over the period 2021-2050 The left panel shows model runs used for the fourth assessment report of the IPCC (AR4) for 5 three different emission scenarios (A1B (53 runs), A2 (36 runs), B1 (44 runs)). The right 6 panel shows model runs that will be used for the fifth assessment report (AR5, all simulations 7 available before 15 December 2011 are included) for four representative concentration 8 pathways (RCP2.6 (26 runs), RCP4.5 (32 runs), RCP6.0 (17 runs), RCP8.5 (29 runs)). The 9 plotted values are means over the values assigned to the grid cells of the climate models over 10 the study region. 11 12 Fig. 3. Range of projected changes (2021-2050 relative to 1961-1990) in yearly average temperature and precipitation in the upstream areas of the Amu and Syr Darya river basins. The left panel shows model runs used for the fourth assessment report of the IPCC (AR4) for three different emission scenarios (A1B (53 runs), A2 (36 runs), B1 (44 runs)). The right panel shows model runs that will be used for the fifth assessment report (AR5, all simulations available before 15 December 2011 are included) for four representative concentration pathways (RCP2.6 (26 runs), RCP4.5 (32 runs), RCP6.0 (17 runs), RCP8.5 (29 runs)). The plotted values are means over the values assigned to the grid cells of the climate models over the study region. and compare it to the period 1961-1990. Hence, the climate change signals refer to the changes during 60 yr. We do this for the simulations in both ensembles. The range of temperature and precipitation projections is shown in Fig. 3 for both ensembles.
We derive the 10th (Q10), 25th (Q25), 50th (Q50), 75th (Q75) and 90th (Q90) percentile values of the changes in precipitation and temperature for each month for the entire CMIP3 and CMIP5 ensemble. We compute a transient "delta change" value for 1961-2050 by linearly interpolating the changes between 1961-1990 and 2021-2050. This is done for every percentile and every month. For each simulated year in 2008-2050, we select a random year from the 1 km × 1 km reference period climate dataset  and we superimpose the basin-averaged monthly temperature and precipitation change values to construct a transient time series from 2008 to 2050. These time series are then used as meteorological forcing for the glacier model, which is run with all the combinations of the percentile values of changes in precipitation and temperature. This well established "delta change" approach (Arnell, 1999;Kay et al., 2008) removes large parts of climate models biases, which cancel out in the climate change signals. We have selected the delta change method as it allows us to include a large number of climate scenarios.

Glacier model
The method used in this study to estimate the glacier evolution is an approach with minimum data requirements. We use a mass balance model with parameterization of glacier area changes and subsequent aggregation of regional glacier characteristics. The model estimates the fractional glacier cover (G F ) for each 1 km grid cell at a monthly time step from 2008 until 2050. The model requires monthly average temperature and monthly precipitation sums, terrain elevation data, the initial fractional glacier cover for each 1 km grid cell and the distribution of glaciers over terrain elevation and glacier size classes as data input (Sect. 3). Figure 4 provides a schematic representation of the modeling steps. First, one basin scale hypsometric curve is derived for the study area, which describes the distribution of glacierized area over the terrain elevation. Subsequently, we calculate a monthly basin scale specific glacier mass balance. We do this by specifying the accumulation area and ablation area using a monthly basin scale 0 • C isotherm and the basin scale hypsometric curve. The model is calibrated against the average of the observed mass balance in the basins during the climatic reference period. For the future, the basin scale mass balance is used to derive an annually updated area for the glaciers in each glacier size class by volume-area scaling (Bahr et al., 1997).
The changes in area are aggregated for all glaciers in all size classes to obtain the basin scale changes in glacier area and construct a basin scale area depletion curve. Finally, the basin scale area depletion curve can be used to calculate an updated fractional glacier cover per 1 km grid cell from 2008 until 2050.

Basin scale hypsometric curve
To generalise the hypsometry of the glaciers in the basins, we construct a basin scale hypsometric curve from the initial fractional glacier cover in the 1 km grid cells. To this end we need to derive the median elevation of the fractional glacier cover (H GLAC ) in a 1 km grid cell. First we use the 90 m DEM to calculate the average terrain altitude (H AVG ), standard deviation of the terrain altitude (H SD ), and maximum terrain altitude (H MAX ) within each 1 km grid cell at the 90 m subgrid. We then derive H GLAC for each grid cell based on the distribution of terrain elevation and G F , assuming that within a 1 km grid cell the distribution of ice follows the terrain elevation distribution and glaciers occupy the highest (coldest) end of the terrain elevation distribution. Figure 5 shows schematically how H GLAC can be determined from H AVG , H SD and G F . It shows the terrain elevation distribution within a 1 km grid cell and the part of the terrain elevation distribution occupied by glacier ice. If we assume the terrain elevation distribution to be approximately normal, then we can estimate the median elevation of the fractional glacier cover as is the 1 − G F 2 quantile of the standard normal distribution and H MAX is the maximum terrain elevation within the 1 km grid cell. H GLAC is limited by H MAX because the median elevation of the fractional glacier cover cannot be higher than the maximum terrain elevation in the 1 km grid cell. If for example H AVG = 4000 m a.s.l., H SD = 200 m and G F = 0.4, then H GLAC = 4168 m a.s.l. When G F = 1, the entire cell is covered with ice and thus Eq. (1) yields H GLAC = H AVG .
We sort the data for H GLAC from low to high values for all grid cells with G F > 0, with each value assigned a weight according to its fractional glacier cover as part of the total glacier area (i.e. the sum of G F for all grid cells) in order to derive one basin scale hypsometric curve (  (1). Using the grid cell's mean terrain elevation (H AVG ) in combination with the standard deviation of terrain elevation within the grid cell (H SD ) and the fractional glacier cover of the grid cell (G F ), the median elevation of the part of the grid cell that is covered by ice can be determined (H GLAC ). Basin scale averaged temperature and elevation for grid cells with glaciers (T AVG and H AVG ) are calculated (2). Values of H GLAC for all grid cells from step 1 are used to construct a basin scale hypsometric curve. Basin scale mass balance calculations are done for all glaciers in 26 glacier size classes with a monthly time step (3). Using H AVG , T AVG and a temperature lapse rate (T lapse ) the basin scale 0 • C isotherm can be determined (H 0 ). By combining H 0 with the hypsometric curve the accumulation area ratio (AAR) can be calculated. With the AAR the amount of ablation (A) and accumulation (C) can be derived. A representative temperature for the ablation zone (T ABL ) is calculated at the mean elevation of the ablation zone (H ABL ). A degree day factor (DDF) is used to calculate the actual ablation. The accumulation consists of the precipitation (P ) over the accumulation zone. Using A and C a monthly mass balance ( M) is calculated. Applying volume-area scaling in October each year an updated glacier area is calculated for the glaciers in each size class and the change in area can be tracked ( S). With the result from step 3 a basin scale area depletion curve is constructed to derive an updated basin scale median elevation of the glacierized part of the basins (H GLAC ) for each month (4). With H GLAC and the elevation distribution within a grid cell (mean terrain elevation (H AVG ) and standard deviation of elevation (H SD )), the basin scale model output is downscaled to the grid cell scale for each month, to provide an updated fractional glacier cover (G F ) per grid cell (5). which represents the average glacier altitude distribution in the study area. We construct the hypsometric curve using the initial fractional glacier cover and distribution of terrain elevation in a 1 km grid cell instead of computing it directly from the glacier outlines and 90 m DEM for consistency with the calculation of the updated fractional glacier cover at the end of each time step during the simulation (Sect. 4.2.5).

Basin scale 0 • C isotherm and accumulation area ratio
Once the basin scale hypsometric curve is obtained, we want to use it to calculate a basin scale monthly mass balance. The idea is to determine the basin scale 0 • C isotherm for each month and combine it with the basin scale hypsometric curve to determine the basin scale accumulation area ratio, which in turn can be used to calculate the ablation and accumulation for each month and for the glaciers in each glacier size class.
To determine the altitude of the basin scale 0 • C isotherm, we calculate the basin scale mean elevation (H AVG ) and the monthly basin scale average temperature (T AVG ). Then, using H AVG and T AVG , we derive the altitude of the basin scale 0 • C isotherm (H 0 ) for each month: where T lapse is a temperature lapse rate ( • C m −1 ), which is the mean of the saturated and dry adiabatic lapse rates (Table 1). for AAR is looked up in the upper horizontal axis of Fig. 6 for the corresponding value of H 0 on the vertical axis. For example, in Fig. 6 H 0 = 4800 m a.s.l. and the associated AAR = 43 % as derived from the basin scale hypsometric curve. The next step is to use the monthly AAR to scale the ablation area and accumulation area for each month, to calculate month specific accumulation and ablation.

Basin scale mass balance
For each month, a specific mass balance ( M [m w.e. yr −1 ]) is determined at basin scale: where C (m) is the monthly accumulation and A (m) is the monthly ablation. The monthly accumulation at basin scale is calculated as where P is the monthly precipitation sum over the glacierized area in the basins (m) and T ACC is the basin scale average temperature representative for the accumulation zone. T ACC can be derived from the median elevation of the accumulation zone at basin scale (H ACC ), which is derived from the hypsometric curve (Fig. 6). For example in Fig. 6, H 0 = 4800 m a.s.l. and AAR = 43 %. Thus the upper 43 % of the glacier area is located in the accumulation zone. The median elevation of this zone (H ACC ) is 0.5 · 43 % = 21.5 % on the AAR-axis. Deriving H ACC from the hypsometric curve yields H ACC = 5133 m a.s.l. We calculate the temperature for the accumulation zone (T ACC ) according to Accumulation occurs when T ACC is below 2 • C as stated in Eq. (4), in which case all precipitation over the accumulation zone is assumed to be solid.

The monthly ablation (A [m]) is calculated as
where T + ABL is the positive (set to zero when negative) basin scale monthly average temperature representative for the ablation zone (see derivation below), DDF is a composite degree day factor (mm w.e. • C −1 day −1 ) calculated as the weighted mean of two distinct values referring to debris free and debris covered ice (Table 1). Weighting is performed according to the fraction of debris free glaciers (85 %) and debris covered glaciers (15 %). This ratio is based on observations in the western Tien Shan . The number of days in the month is d, and AAR is the accumulation area ratio. The degree day factors for debris free glaciers and debris covered glaciers were calibrated in a related hydrological study for the same river basins (Immerzeel et al., 2012a).
To calculate T ABL we derive the median elevation of the ablation zone at basin scale (H ABL ) using the hypsometric curve and the AAR. For example in Fig. 6, H 0 = 4800 m a.s.l. and AAR = 43 %. Thus the lower 57 % of the glacier area is part of the ablation zone. The median elevation of this zone (H ABL ) is 100 %-0.5 × (100 %-AAR) = 71.5 % on the AAR-axis. Deriving H ABL in the hypsometric curve yields H ABL = 4261 m a.s.l.
We calculate the temperature for the ablation zone (T ABL ) according to For each month a specific mass balance is calculated at basin scale as specified in Eq. (3).

Updating glacier area for glaciers in each size class
An intitial mean ice thickness is determined for the glaciers in each size class using volume-area scaling (Bahr et al., 1997). Volume-area scaling is based on physical arguments (Bahr et al., 1997) and has been extensively used (e.g. Farinotti et al., 2009;Radić and Hock, 2010;Grinsted, 2013). The volume-area scaling can be expressed as a relation between the mean glacier thickness (h [m]) and glacier area (A [m 2 ]) (Radić and Hock, 2010;Huss and Farinotti, 2012): where c and γ are scaling parameters. We use the same scaling parameters as Radić and Hock (2010) use for mountain glaciers (c = 0.2055, γ = 1.375). With this relation we derive an initial mean ice thickness for the glaciers in each size class. This thickness is updated every month (t) for the glaciers in each size class (i) with the basin scale specific mass balance (Sect. 4.2.3): To simulate future glacier extent, we force the model with the downscaled temperature and precipitation projections described in Sect. 3.3 for 2008 until 2050 at a monthly time step. Each year at the beginning of a new glaciological year (in this study on October 1st), we use the inverse of Eq. (8) to calculate the new glacier area for each size class from the updated ice thickness. By aggregating the results for all glaciers in all size classes, the percentile change in total glacierized area in the basins from 2008 to 2050 with respect to 2007 is determined (Fig. 7). An area depletion curve can be fitted through the time series of percentile changes in glacier area (Fig. 7). By looking up the H values on the vertical axis in Fig. 6 that correspond to the values of the area depletion curve for each time step on the upper horizontal axis, a time series of updated basin scale median elevation of the glacierized part of the basins (H GLAC ) (Fig. 7) is constructed, which can later be used to downscale the basin scale averaged changes in glacier area to monthly updated fractional glacier cover for each 1 km grid cell. For example in Fig. 7, on 1 January 2040 the glacierized area is 52.0 % of the glacierized area in 2007 as can be derived from the area depletion curve. Using the fractional glacier cover value 0.520 (=52.0 %) in the lower horizontal axis in Fig. 6 yields 4586 m a.s.l. for H GLAC from the hypsometric curve.

Updating fractional glacier cover per grid cell
To create monthly maps of glacier extent, we update the fractional glacier cover (G F ) for each grid cell for each month from 2008 until 2050 using H GLAC and the distribution of terrain elevation within a 1 km grid cell. Assuming that the glacier distribution follows the distribution of terrain elevation, and that the latter can be described by a normal distribution, we calculate G F for a 1 km grid cell using the cumulative standard normal curve function: For example in Fig. 5, when H GLAC = 4270 m a.s.l. and a given grid cell has H AVG = 4000 m a.s.l. and H SD = 200 m, then G F = 0.18. If H GLAC moves up, G F decreases. G F has an upper limit of 1, as the fractional glacier cover cannot exceed this value. Thus, when H GLAC ≤ H AVG , G F = 1.

Calibration
We calibrate the model for the reference period .
Based on the average of the observed mass balance in the region during the reference period (Sect. 3.4, Table 2) the model is calibrated by correcting the monthly mean temperature for the reference period with a temperature correction ("CorT") (Table 1), which is added to the temperature forcing. With the calibrated CorT, the model produces the same mass balance for the reference period as the average of the observed mass balance in the basins (MB OBS , Table 1). The CorT parameter accounts for a combined effect of errors in the forcing data, temperature differences within a 1 km grid cell, vertical and horizontal errors from interpolation in the reference period climate dataset (Sect. 3.2) and errors from averaging over the two basins. The degree day factors for debris free glaciers and debris covered glaciers where calibrated for a related hydrological study for the same river basins (Table 1) (Immerzeel et al., 2012a). The degree day factors are within the range of other studies reported in the region (Mihalcea et al., 2006;Zhang et al., 2006;Hagg et al., 2008;Immerzeel et al., 2010Immerzeel et al., , 2012b. In addition we take into account variation in degree day factors in the uncertainty analysis described in Sect. 5.3.

Validation
Since data scarcity in Central Asia makes it difficult to validate the model performance, we validate the method for the Austrian Alps, where multiple glacier inventories and glacier mass balance time series for twelve glaciers are available. We use two glacier inventories, marking the starting point and endpoint of the simulation. A glacier inventory representative for the year 1969 (Patzelt, 1978)  Kuhn, 2007). We assume this inventory to be representative for 1997, since 81 % of the glacier area was mapped in 1997 and 1998. Thus, 1997 is the last year of the simulation. We force the model with daily air temperature and daily precipitation from the PGMFD (Sheffield et al., 2006). We use the same DEMs, the same degree day factors and the same volume-area scaling coefficients as used for the application in Central Asia. The average of the observed mass balance in the Austrian Alps is −0.37 m w.e. yr −1 between 1969 and 1997 based on mass balance records from twelve individual glaciers (WGMS, 2011). We calibrate the CorT parameter to this average of the observed mass balance yielding CorT = 0.76 • C and simulate the changes in glacier area until 1997. The simulated decrease in glacier area between 1969 and 1997 is 24.5 %. Figure 8 shows the complete simulation of changes in glacier area from 1969 until 1997. The black error bars indicate simulation results when calibrated for the average of the observed mass balance plus one standard deviation (positive error) and the average of the observed mass balance minus one standard deviation (negative error). The observed decrease in glacier area according to the two glacier inventories equals 19.0 % (Fig. 8). The estimated error in the glacier inventories (Lambrecht and Kuhn, 2007) is displayed with the blue error bars. Considering the fact that our approach is a first order estimate of basin scale glacier area the change over 60 years (1961-1990 to 2021-2050). The boxes represent the range from Q25 10 to Q75, divided by the median value (Q50). The whiskers represent the range between Q10 11 and Q25 (at the lower end of the distributions) and the range between Q75 and Q90 (at the 12 higher end of the distributions). 13 14 Fig. 9. Box-whisker plots for projected changes in temperature (left) and precipitation (right) for three AR4 SRES emission scenarios and four AR5 representative concentration pathways extracted from the CMIP3 (SRES) and CMIP5 (RCP) databases. The A1B (53 GCM runs), A2 (36 runs) and B1 (44 runs) AR4 scenarios are used and the RCP2.6 (26 runs), RCP4.5 (32 runs), RCP6.0 (17 runs) and RCP8.5 (29 runs) AR5 scenarios are used. The values are mean delta change values for GCM grid cells covering the study area and represent the change over 60 years (1961-1990 to 2021-2050). The boxes represent the range from Q25 to Q75, divided by the median value (Q50). The whiskers represent the range between Q10 and Q25 (at the lower end of the distributions) and the range between Q75 and Q90 (at the higher end of the distributions). changes, the uncertainties in the methodology (as discussed in Sect. 5.4) and uncertainties in the glacier outlines in the inventories, we conclude that the model performs satisfactory.

Future climate
All results stated are for the Amu Darya and Syr Darya basins combined and the climate change signals refer to the changes during 60 yr (change between 1961-1990 and 2021-2050). Both the CMIP3 and CMIP5 ensembles show large variation in temperature and precipitation changes between models and between emission scenarios (Fig. 9). On average, temperature is expected to rise by about 2 • C and precipitation to remain nearly constant. The uncertainty in temperature projections ( T ), expressed as the 90th and 10th percentiles, is estimated to range from 1.3 to 2.4 • C in the CMIP3 ensemble and from 1.7 to 2.9 • C in the CMIP5 ensemble (Fig. 9, left panel). For precipitation projections ( P ) the 90th and 10th percentiles range from −6 to +7 % in the CMIP3 ensemble and from −8 to +15 % in the CMIP5 ensemble (Fig. 9, right panel). Though the climate projections of both ensembles mainly cluster around the same values (about 2 • C and 0 %, for temperature and precipitation, respectively), the new CMIP5 ensemble includes the possibility of more extreme climate change. There are several "warmer" simulations (up to +3.5 • C) and many of those are also extreme The definition of the boxplots is as in Figure 9. 4 5 Fig. 10. Box-whisker plots for projected changed per month in temperature (upper panel) and precipitation (lower panel) for the CMIP3 ensemble (red) and CMIP5 ensemble (blue). The definition of the boxplots is as in Fig. 9.
in precipitation change (Fig. 9). Note that this observation not only holds across scenarios, but also between GCM runs within a given scenario, e.g. RCP 2.6, 6.0 and 8.5 show similar extremes in temperature and precipitation. The CMIP5 ensemble also shows a larger average warming than CMIP3 (Fig. 9, left panel). In addition, the variation between scenarios is also larger for CMIP5 for both precipitation and temperature (Fig. 9). Looking at the projections on a monthly scale (Fig. 10), mean projections for temperature ( T , Q50) in July to September do not differ much between the two ensembles, although the range in temperature projections is higher for the CMIP5 ensemble compared to the CMIP3 ensemble (Fig. 10, upper panel). However, mean temperature projections for October to May are higher for the CMIP5 ensemble compared to the CMIP3 ensemble. The spread in precipitation projections is generally larger for the CMIP5 ensemble compared to the CMIP3 ensemble (Fig. 10, lower  panel). Especially for March to September the mean projections for precipitation ( P , Q50) are higher for the CMIP5 ensemble compared to the CMIP3 ensemble, while little differences in mean projections for precipitation are observed for October to February.
As we choose to include as many climate projections as possible in our study, we do not use particular GCMs but apply the quantile approach as described in Sect. 4.1. A disadvantage of this approach is that systematic changes in the daily variability are not included. However, since our glacier model is forced with monthly data, we accept this for the benefit of including as many climate projections as possible.

Implications of climate change for Central Asian glaciers
We force the glacier model with all quantile combinations of the downscaled temperature and precipitation as analysed on a monthly basis (monthly delta change values) and obtain the basin scale cumulative mass balance for the simulated period . Figure 11 shows the cumulative mass balance for the average projection ( T Q50, P Q50), for the very warm and very dry ( T Q90, P Q10) case, for the very cold and very wet ( T Q10, P Q90) case, for the very warm and very wet case ( T Q90, P Q90) and for the very cold and very dry ( T Q10, P Q10) case. The range of projections is higher for the CMIP5 (Fig. 11, right panel) ensemble compared to CMIP3 (Fig. 11, left panel). When forced with the CMIP3 ensemble the cumulative mass balance for 2007-2050 ranges from −32.3 m w.e. for the very cold, very wet case to −44.9 m w.e. in the very warm, very dry case. Forcing with the average projection yields −38.9 m w.e. When forced with the CMIP5 ensemble the range is from −32.2 m w.e. to −47.7 m w.e. for the very cold, very wet case and the very warm, very dry case, respectively. Forcing with the average projection yields a cumulative mass balance for 2008-2050 of −38.6 m w.e. Figure 12, spanning the frequency space between the 10th and 90th percentiles for both temperature and precipitation, shows the percentile glacier retreat in 2050 for the CMIP3 and the CMIP5 case. Both cases show variability in future glacier extent. For the CMIP3 projections, a reduction in glacier area varying between 54.5 % in 2050 when the model is forced by the T Q10 and the P Q90, and a reduction of 63.5 % in 2050 when forced by the T Q90 and P Q10 is observed. By keeping T constant at the Q50 level a 0.8 % range in potential glacier retreat is found (from 59.0 to 59.8 % decrease) over the full P range for the CMIP3 case and a range of 6.7 % is found (from 56.0 to 62.7 % decrease) when P is kept constant at the Q50 level. For the CMIP5 case this range is larger with a 1.1 % range (from 59.1 to 60.2 % decrease) when T is kept constant at the Q50 level, and a 7.8 % range (from 55.7 to 63.5 % decrease) when P is kept constant at the Q50 level. So, the range in temperature projections has a much larger impact on the predicted glacier extent as the range in precipitation projections.
The range for the CMIP5 based projection for glacier extent is slightly wider than for CMIP3. The T Q10 and the P Q90 combination results in a projected decrease of 54.4 %, while the T Q90 and the P Q10 combination leads to a decrease of 65.1 % (Fig. 12). Although the mean temperature projection on an annual basis is higher for the CMIP5 ensemble compared to the CMIP3 ensemble and the mean precipitation projections are almost similar, the projected decrease in glacier extent is practically the same (even 0.1 % more decrease for the CMIP3 case). This can be explained by the fact that mean temperature projections ( T Q50) for July to September, when most of the melting Hydrol. Earth Syst. Sci., 17, 3661-3677 takes place, are similar for CMIP3 and CMIP5 (Fig. 10), and mean precipitation projections ( P Q50) are also similar for CMIP3 and CMIP5 during October to February when most accumulation takes place. From these results it is evident that it is important to assess climate change projections at the seasonal level rather than at the annual level, when making projections for future glacier extent. Figure 13 shows the decrease in total glacier area in the Amu Darya and Syr Darya basins for the entire simulated period based on the CMIP3 (Fig. 13, left panel) and CMIP5 (Fig. 13, right panel) model runs. The range of glacier extent projections for the CMIP5 ensemble and the CMIP3 ensemble are very similar. The fact that the very cold, very dry projection is closer to the very cold, very wet projection than to the average projection for both ensembles again shows that the uncertainty in temperature projections has a much larger impact on the uncertainty in glacier extent than uncertainty in precipitation projections, and change in temper-ature is the main driver for future decrease in glacier extent in these areas. Figure 14 shows, for the CMIP5 case, the projected glacier extent for 2050 for a selected area covering the large glacier systems in the central Pamir (Fig. 14b) as compared with the initial glacier extent (Fig. 14a). The three lower left panels (Fig. 14c, e, g) show the projected fractional glacier cover per 1 km grid cell in 2050 for the average projection and the two most extreme projections (very cold, very wet and very warm, very dry). The right panels (Fig. 14d, f, h) show the change in fractional glacier cover per 1 km grid cell with respect to the initial situation for these three cases. It can be clearly seen that the fractional glacier cover decreases strongest in the lowest glacierized parts, and that mainly the tongues in the valleys are affected. A similar figure shows a selected area in the Tien Shan (Fig. 15). In the Tien Shan mountains the glaciers are smaller than in the central Pamir and many are located at lower elevations. As a result, in the Tien Shan the impact of climate change will lead to a more rapid decrease in glacier extent than in the Pamir.

Parametric uncertainty analysis
Besides uncertainty in glacier extent as a result of the uncertainty in the climate change projections, the projected glacier changes are subject to other uncertainties. These include parametric uncertainty, uncertainty in present glacier extent and volume, uncertainty in the volume-area scaling, uncertainty in climate evolution, uncertainty in climatic forcing for the reference period, uncertainty in mass-balance time series and uncertainties stemming from simplifications and assumptions applied to the model. Since the mass balance model is based on an empirical approach requiring calibration we also evaluate, besides uncertainty in the climate change projections, how the uncertainties in the model parameters as well as uncertainty in the observed historical glacier mass balance translate in uncertainty in the future glacier extent by running the model for different sets of parameters and observed glacier mass balance. We assume the three critical model parameters (vertical temperature lapse rate (T lapse ), degree day factor for clean ice glaciers (DDF CI ), degree day factor for debris covered glaciers (DDF DC )) to be three independent normally distributed (random) variables. The temperature correction (CorT) is recalibrated for each set of parameters. We use a mean DDF DC = 3.97 mm • C −1 d −1 and DDF CI = 7.94 • C −1 d −1 and both with σ = 1 • C −1 d −1 . For T lapse we use a mean −0.0068 • C m −1 and assume a standard deviation of 0.0012 • C m −1 , which is based on the difference between the dry and saturated adiabatic lapse rate. The average of the observed glacier mass balance (MB OBS ) used is −0.47 m yr −1 with a standard deviation of 0.082 m yr −1 (Sect. 3.4, Table 2). For the observed mass balance we use an uncertainty range of two standard deviations. Based on these assumptions we sample 50 parameter sets and mass balance values. We then run a full simulation until 2050 with each of these 50 parameter-mass balance combinations (i.e. of T lapse , DDF CI , DDF DC , MB OBS and associated CorT, which is separately calibrated for each combination) and we estimate uncertainty by taking the standard deviation of the 50 simulations (Ragettli and Pellicciotti, 2012). This analysis allows to estimate the propagation of parameter uncertainty to uncertainty in the glacier model simulations.
The uncertainty resulting from model parameters is displayed for the very cold, very wet and the very warm, very dry cases in Fig. 13. The effect of model parameter uncertainty leads to an additional uncertainty of ±8.6 % in total glacier extent in 2050 for both the CMIP5 and the CMIP3 case, showing that parameter uncertainty has roughly the same importance as uncertainty in the climate projections.

Limitations in the methodology
The advantage of low data requirements associated with the approach described in this paper of course comes with its limitations. We use volume-area scaling to estimate the initial ice volume based on the initial glacierized area and to translate new ice volumes to areas (Bahr et al., 1997). Approaches that use volume-area scaling are sensitive to the scaling parameters used (Grinsted, 2013), but have been largely used for large areas. Other methods based on ice physics and fluxbalance principles have been suggested to estimate the initial ice volume (Farinotti et al., 2009;Huss and Farinotti, 2012;Paul and Linsbauer, 2012), which could yield different results when applied in our modeling study.
We are interested in simulating the behavior of the glaciers as a result of climate perturbations at the basin scale. We do not model individual glaciers, and therefore we use an average of the observed mass balance for the five glaciers in calibration. This regionalization is justifiable over a longer period, but not at smaller time steps.
In our model set-up, we construct one average hypsometric curve for the two river basins. This simplification constitutes represents the area enlarged in the other panels. Panel (a) shows the initial fractional glacier cover per 1 km grid cell. The three lower left panels show the simulated fractional glacier cover per 1 km grid cell in 2050 for the CMIP5 runs. Panel (c) shows the fractional glacier cover for the run with the 50th percentile (Q50) values of temperature and precipitation change. Panel (e) shows the projection with the strongest decrease in glacier cover, when the model is forced with the 90th percentile (Q90) for temperature change and 10th percentile (Q10) for precipitation change. Panel (g) shows the projection with the least decrease in glacier cover, when the model is forced with the 10th percentile (Q10) for temperature and 90th percentile (Q90) for precipitation change. The three lower right panels (d, f, h) show the change in fractional glacier cover per grid cell for the 2050 projections in the three lower left panels (c, e, g) with respect to the initial glacier cover (a). a drawback as regional differences are neglected. To retain more regional differences a more accurate glacier modeling could be done by constructing different hypsometric curves for different (sub)basins, or theoretically for every grid cell. The same holds for basin scale averaged temperature and precipitation. As we use the initial hypsometric curve during the entire simulation, another improvement could be inclusion of regular recalculation of the hypsometric curve during the 46 1 Fig. 15. Projected fractional glacier cover in 2050 for a slected area in the Tien Shan mountains. The square in the top right panel (b) represents the area enlarged in the other panels. Panel (a) shows the initial fractional glacier cover per 1 km grid cell. The three lower left panels show the simulated fractional glacier cover per 1 km grid cell in 2050 for the CMIP5 runs. Panel (c) shows the fractional glacier cover for the run with the 50th percentile (Q50) values of temperature and precipitation change. Panel (e) shows the projection with the strongest decrease in glacier cover, when the model is forced with the 90th percentile (Q90) for temperature change and 10th percentile (Q10) for precipitation change. Panel (g) shows the projection with the least decrease in glacier cover, when the model is forced with the 10th percentile (Q10) for temperature and 90th percentile (Q90) for precipitation change. The three lower right panels (d, f, h) show the change in fractional glacier cover per grid cell for the 2050 projections in the three lower left panels (c, e, g) with respect to the initial glacier cover (a). simulation based on the updated fractional glacier cover per grid cell.
Another area for improvement is the melt modeling. We now use a combined degree day factor for debris free and debris covered glaciers, which reflects the different behaviour of the two surfaces with melt decreasing under a thick debris cover (Nicholson and Benn, 2006;Brock et al., 2010).
If the exact extent of both types of glaciers is available it would be recommendable to model the two types separately. However, melt modeling under debris covered glaciers is not trivial as it crucially depends on debris thickness, which is not commonly available. Strong spatial variation is observed in the Alps as a result of the type and thickness of the debris layer. Improved models for melt under debris should be used that account for the effect of debris thickness (Reid et al., 2012), provided that the thickness and characteristics of the debris layer are known. Apart from modeling melt under debris cover, melt modeling can be improved by including incoming solar radiation (e.g. Pellicciotti et al., 2005), and considering other components of the energy balance. A general limitation of degree day melt models is the necessity to calibrate the parameters for each case as the parameters are not transferable in time and space.
Given the limitations discussed above, we are aware that the glacier model used in this study is too coarse to reproduce the response of single glaciers and the complexity of processes involved. The model choice is imposed by the limited amount of data available and the large scale of our application. However, the model is suitable for our aim, i.e. to translate downscaled future climate scenarios into glacier response at the basin scale, and to assess how the spread and differences in the future climate scenarios transform into differences in glacier response. Despite the fact that the simulated glacier response is subject to the uncertainties we discuss, the simulated trends are apparent. This study shows that parameter uncertainty and differences between GCMs should be taken into account and that the impact of climate change signals should take account of seasonal variation.

Conclusions
Both CMIP3 and CMIP5 climate change projections point towards a decline of glacier extent in Central Asia. Our results show that uncertainty about the range of this decline remains large. The range of projections for temperature and precipitation in the Central Asian region until 2050 for the CMIP5 ensemble is larger than for the CMIP3 ensemble and the median projection for CMIP5 models shows greater warming than for CMIP3 models. The CMIP5 ensemble shows higher projections for winter temperatures compared to CMIP3 while summer temperature projections are similar. On the other hand, the CMIP5 ensemble shows higher precipitation projections for the summer months compared to CMIP3 ensemble, while precipitation projections for the winter months are similar for both ensembles. As a result, the CMIP5 ensemble leads to a slightly wider range in projected glacier extent. For temperature and precipitation projections, the median projection shows a decrease in glacier extent between 2007 and 2050 of 59.4 % for the CMIP5 ensemble compared to 59.6 % in the CMIP3 case. The projected decrease in glacier extent ranges from 54.4 to 65.1 % for the CMIP5 ensemble compared to 54.5 to 63.5 % for the CMIP3 ensemble. Large spread is evident among models within both ensembles, in agreement with recent studies that have indicated that the differences among GCMs due to their structure and characteristics is the main source of uncertainty in future climate. Parametric uncertainty leads to additional uncertainty in the projections of future glacier extent, and has roughly the same importance as uncertainty in the climate projections. The mentioned ranges in projected glacier extent decrease demonstrate substantial uncertainty in climate change projections and associated glacier response for Central Asia. Furthermore, it shows that it is imperative to use a representative selection of climate models and emission scenarios that span the entire range of possible future climates in climate change impact studies, to provide a complete picture of possible climate change impact. At the same time it shows that climate change signals should be analysed at a seasonal scale, when used to assess the response of glaciers to the changes in climate. The wide range in the projections implies an uncertain future for Central Asian glaciers.