Projections for headwater catchments of the Tarim River reveal glacier retreat and decreasing surface water availability but uncertainties are large

For model calibration and evaluation over 1957–2004, precipitation and temperature were interpolated from observed station measurements. Radiation and humidity data were retrieved from the Watch Forcing Data, based on ERA-40 and ERA-Interim (Uppala et al 2005, Dee et al 2011, Weedon et al 2011, 2014). The methods used to interpolate the meteorological data are reported by Duethmann et al (2015).

For the scenario analyses, we used daily climate projections from nine CMIP5 GCMs forced by three emission scenarios (RCP2.6, RCP4.5, and RCP8.5). The GCMs were selected to represent the range of changes projected by the CMIP5 ensemble in the study area (supplementary figure S1). With a spatial resolution of about 1.5°–3°, the GCMs cannot capture the small-scale topographic features. A better spatial resolution can be achieved through downscaling by RCMs. Data were available from two RCMs: the RCM REMO forced by the CMIP3 GCM ECHAM5 under the A1B emission scenario with a resolution of 0.166° (Mannig et al 2013), and the RCM CCLM forced by the CMIP5 GCM MPI-ESM-LR under the emission scenarios RCP2.6, RCP4.5, and RCP8.5 with a resolution of 0.44° (Wang et al 2013). The RCM based scenarios were used for impact assessment but to keep a consistent ensemble they were not included for uncertainty evaluation.

All climate model data were summarised to subcatchments and bias corrected using an empirical quantile mapping approach (Gudmundsson et al 2012), which allowed us to also consider changes in the variability. The bias correction was applied separately for each month.

Discharge and glacier geometry changes were simulated with the hydrological model WASA (Güntner and Bronstert 2004), which has also previously been applied in Central Asia (Duethmann et al 2013, 2014, 2015). The implementation, calibration, and evaluation of the WASA model to the Aksu headwater catchments is presented in detail in Duethmann et al (2015) and summarised in the following.

The model is applied at a daily time step, and the spatial discretization is based on hydrologic response units (HRUs) defined by subcatchments and elevation zones. Each HRU is characterised by its dominant land cover (based on MODIS land cover MOD12Q1 collection 5.1; Friedl et al 2002), and soil (Harmonized World Soils Database; FAO/IIASA/ISRIC/ISS-CAS/JRC 2012), its average elevation (derived from the SRTM DEM; Jarvis et al 2008), and its glacier fraction (based on the glacier inventory for ~1975 (~2008) for the calibration/control (future) simulations; Osmonov et al 2013, Pieczonka and Bolch 2015). The initial glacier ice thickness distribution of each glacier was based on estimates by a spatially distributed ice-thickness model (GlabTop2) (Linsbauer et al 2012, Frey et al 2014), with the SRTM DEM and the 1970s glacier inventory as inputs.

Snow and glacier melt is simulated using a temperature-index approach with seasonally varying melt factors. For each glacier, area and thickness is updated annually based on the simulated glacier mass balance using the Δh-approach (Huss et al 2010) (for details and model equations, refer to the supplement of Duethmann et al 2015). For model calibration, glacier area changes were based on the glacier inventories for ~1975 and ~2008 assuming linear decrease rates for each glacier.

A comprehensive multi-variable approach was used for model calibration and evaluation. Daily discharge variations were calibrated using an average of the Nash–Sutcliffe efficiency of linear and logarithmic discharge values to achieve a balanced evaluation of high and low discharges (equation (1))

$$\begin{eqnarray}&&{f}_{{Q}_{1}}=0.5\cdot \left(\left(1-\displaystyle \frac{\displaystyle {\sum }_{t=1}^{T}{({Q}_{{\rm{obs}}}(t)-{Q}_{{\rm{sim}}}(t))}^{2}}{\displaystyle {\sum }_{t=1}^{T}{({Q}_{{\rm{obs}}}(t)-\bar{{Q}_{{\rm{obs}}}})}^{2}}\right)\right.\\ &&\left.+\left(1-\displaystyle \frac{\displaystyle {\sum }_{t=1}^{T}{(\mathrm{log}({Q}_{{\rm{obs}}}(t))-\mathrm{log}({Q}_{{\rm{sim}}}(t)))}^{2}}{\displaystyle {\sum }_{t=1}^{T}{(\mathrm{log}({Q}_{{\rm{obs}}}(t))-\bar{\mathrm{log}({Q}_{{\rm{obs}}})})}^{2}}\right)\right).\end{eqnarray} \tag{ 1 }$$

${Q}_{{\rm{obs}}}(t)$ and ${Q}_{{\rm{sim}}}(t)$ are the observed and simulated daily discharge at time t, and T is the number of time steps. As high values of ${f}_{{Q}_{1}}$ may be achieved despite low performance with respect to interannual discharge variations, we also evaluated interannual variations of seasonal flow using the same criterion but applied to annual series of seasonal (DJF, MAM, JJA, SON) flow $({f}_{{Q}_{2}})$.

For the glacier mass balance, the following criteria were applied. The cumulative glacier mass change over the period 1976–1999 was constrained to −0.52 ± 0.33 m w.e. a⁻¹ based on two geodetic mass balance estimates (Surazakov and Aizen 2006, Pieczonka et al 2013) (figure 1). The criterion was represented by a function that reaches an optimal value of one for the average mass balance estimate, a value of zero outside the uncertainty range, and linearly increasing/decreasing values in between. The temporal variation of the simulated glacier mass balances was calibrated by maximising the correlation to an in situ glacier mass balance series of Karabatkak Glacier (Dyurgerov and Meier 2005, WGMS 2012), located close to the study area (figure 1).

The model was automatically calibrated to the two discharge and the two glacier mass balance criteria with a multiobjective calibration algorithm (ε-NSGAII; Kollat and Reed 2006), using 24 years (1976–1999) for model calibration, and 24 years (1957–1975 and 2000–2004) for model validation. From the results of the multiobjective optimisation, solutions were selected if they showed a good performance for daily and interannual streamflow variations (${f}_{{Q}_{1}}\gt 0.75\;,$ and ${f}_{{Q}_{2}}\gt -2),$ were inside the uncertainty range of the geodetic glacier mass balance estimate, and represented observed discharge trends in the past. Since the calibration period was regarded as too short for the calculation of trends, we calculated trends in seasonal discharge time series over the entire historical period 1957–2004. Trend evaluation was constricted to trend direction and trend significance. As a result, 6 and 28 parameter sets were selected for the Sari-Djaz and Kakshaal catchment, respectively.

Nash-Sutcliffe efficiencies for daily (monthly) discharge are within 0.77–0.87 (0.81–0.94) in the calibration and 0.67–0.84 (0.80–0.92) in the validation period (ranges over both catchments and the selected parameter sets). The model shows an acceptable performance with respect to the interannual variations of seasonal discharge and represents the observed discharge trends. Trends over 1957–2004 are 19–21 mm decade⁻¹ (simulated) and 21 mm decade⁻¹ (observed) for the Sari-Djaz catchment (mean annual runoff 382 mm a⁻¹), and 10–11 mm decade⁻¹ (simulated) and 10 mm decade⁻¹ (observed) for the Kakshaal catchment (mean annual runoff 151 mm a⁻¹) (trends evaluated with Sen's slope estimator Sen 1968). In agreement with the observed geodetic glacier mass balance estimates, the simulated glacier mass balance estimates over 1976–1999 are in a range of −0.83 to −0.51 m w.e. a⁻¹ (−0.85 to −0.25 m w.e. a⁻¹) for the Sari-Djaz (Kakshaal) catchment.

In this study, we analysed uncertainties from three sources, including three emission scenarios, nine GCMs, and six (28) parameter sets of the hydrological model for the Sari-Djaz (Kakshaal) catchment. This led to 162 ensemble members for the Sari-Djaz catchment and 756 for the Kakshaal catchment. Uncertainties are characterised by the 5–95 percentile range of the ensemble. For consistency, only GCM-based simulations were included in the ensemble. Contributions of individual uncertainty sources were analysed by their contribution to the total variance (see supplement).

The analysis of changes in climate and glacio-hydrologic response was carried out by comparing the near future 2010–2039 ('2020s'), mid future 2040–2069 ('2050s'), and far future period 2070–2099 ('2080s') to the control period 1971–2000. Increased air temperatures are projected by all climate models under all emission scenarios (figure 2). Differences in the projected temperature increase between the different RCPs are negligible for the '2020s' but substantial for the '2080s'. Temperature increases at a similar rate over the whole 21st century under RCP8.5, while for RCP4.5, the rate of temperature increase is projected to slow down toward the end of the 21st century, and in RCP2.6, temperature is projected to stabilise or decrease about the end of the century. The resulting temperature increase in the '2080s' compared to 1971–2000 is 3.2 °C–7.2 °C for RCP8.5, 1.6 °C–4.5 °C for RCP4.5 and 0.7 °C–2.8 °C for RCP2.6. Seasonally, the strongest temperature increases are expected for winter (supplementary figure S2).

Zoom In Zoom Out Reset image size

Precipitation is projected to be about 10% higher for the '2020s' as compared to 1971–2000 for all RCPs, but the variability between GCMs is high. Over the 21st century, most GCMs project further small increases in precipitation. In contrast to temperature, differences between the RCPs remain small.

Compared to the median of the GCMs, the two RCMs show a tendency toward lower precipitation and temperature increases, but are mostly within the interquartile range of the GCMs. From the ensemble, we depicted three GCMs representing warm-wet (MIROC-ESM), warm-dry (IPSL-CM5A-LR), and cold (GFDL-ESM2M) conditions to better comprehend the change mechanisms of the glacier and hydrological variables described in the next sections.

Changes in humidity and radiation are mostly small and within ±5%. Larger decreases in humidity are linked to the warm-dry climate projections (supplementary figure S4).

The projected changes in temperature and precipitation are simulated to cause predominantly negative mass balances and significant losses in glacier extent (figure 3) and volume (supplementary figure S5). By the end of the century, projections suggest a reduction in glacier area by −66(−89 to −28)% for the Sari-Djaz and −78(−94 to −47)% for the Kakshaal catchment (uncertainties refer to the 5–95 percentile range). The variability of the GCM ensemble contributes most to the overall uncertainty (figure 4). As expected, the strongest glacier area reduction is projected for the warm-dry GCM (figure 3). The influence of the emission scenario increases from the 2050s, and uncertainties from the GCMs and emission scenarios are similar by the end of the century. The pattern at the beginning of the scenario period in figure 4 relates to very small total uncertainties (figure 3) and is therefore not analysed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The glacier area reduction by the CCLM simulations is slightly lower than the ensemble median, in accordance with the lower than average temperature increase (figure 3). The REMO simulations suggest a strong glacier reduction from the 2050s onwards, which is in line with the continuous temperature increase under the A1B scenario.

In the '2020s', discharge is generally projected to be higher than during 1971–2000 in all seasons (figure 5). In both catchments and over the whole 21st century, discharge in spring and early summer (April–June) is higher than during the control period for nearly all ensemble members. The largest changes in relative terms are observed during spring, where discharge during the reference period is low (supplementary figure S6). Discharge changes in late summer are more heterogeneous. In the Sari-Djaz catchment, August discharge is increased in the '2020s' but decreased in the '2080s'. For the Kakshaal catchment, both increases and decreases are observed during all scenario periods, and the direction of change strongly depends on the choice of GCM, which contributes most to overall uncertainties (figure 6). By the end of the century, emission scenarios are an important uncertainty source during spring and early summer, while the contribution of the hydrological model parameters generally remains small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Snowmelt is shifted toward earlier in the season leading to higher snowmelt in spring and lower snowmelt in summer as shown for the '2080s' (figure 7). Increase in spring rainfall further contributes to higher discharge during spring. During the '2080s', lower summer discharge is caused by the shift in snowmelt, decreases in ice melt, and increases in actual evapotranspiration (AET). Summer rainfall shows a large variability in the GCM projections. A decrease in summer rainfall as in the warm-dry GCM may further exacerbate the discharge decreases, while rainfall increases as projected under the warm-wet GCM partly compensate the changes in snowmelt, ice melt, and AET.

Zoom In Zoom Out Reset image size

Glacier melt (defined as ice melt from glaciers) is simulated to change substantially over the 21st century in the Sari-Djaz catchment (figure 8(a)). After increases at the beginning of the 21st century, glacier melt drastically decreases during 2040–2060. Peak glacier melt is projected to occur roughly around 2030 (figure 8(a)). In the Kakshaal catchment, glacier melt is generally much lower and peak melt is passed earlier (figure 8(b)).

Zoom In Zoom Out Reset image size

The changes in glacier melt have a significant impact on the evolution of total annual runoff of the Sari-Djaz River, resulting in the reduction of discharge by −22(−51 to +3)% in the '2080s' compared to the control period (figure 9(a), table 1, supplementary figure S7). In contrast, the impact of changes in glacier melt is smaller for the Kakshaal River due to its lower glacier melt contribution to discharge (figure 9(b), table 1, supplementary figure S8).

Zoom In Zoom Out Reset image size

The projected evolution of glacier melt for the RCM-based simulations is similar to those of the GCM-based simulations, while the projected discharge toward the end of the century is slightly lower than the median of the GCM-based simulations. This is in line with a lower precipitation increase than projected by the median of the GCMs.

Changes in the interannual variability of discharge were evaluated using the coefficient of variation (CV). The CV was calculated after applying a 20 year moving average filter to exclude influences of systematic discharge changes. In both catchments, a tendency toward increased interannual variability of discharge can be seen; however, considering the large variability for different ensemble members, the changes are only small (figures 10(a) and (d)). A moderate glacier cover of a catchment is linked with low interannual discharge variability: in years with low precipitation, glacier melt begins earlier and can therefore (partly) compensate lower snowmelt (Jansson et al 2003, Hock et al 2005). The increase of interannual variability may therefore be caused by decreasing glacier area and glacier melt. In addition, changes in the interannual variability of climate variables might also play a role; there is a slight increase in the interannual variability of temperature (figures 10(c) and (f)), but changes in interannual variability of precipitation are small (figures 10(b) and (e)).

Zoom In Zoom Out Reset image size

Our results for glacier area reduction of −36(−69 to −12)% by 2050 and −69(−90 to −32)% by 2099 compare well with other studies for Central Asia (supplementary table S1). The projected percentage glacier area loss is larger in the Kakshaal than in the Sari-Djaz catchment, where glaciers are generally larger and have a higher thickness.

Changes in the runoff regime with lower runoff in August and higher runoff in May/June are typical for mountain catchments influenced by snow and glacier melt and have been projected for many regions including other catchments in Central Asia (Hagg et al 2007, Sorg et al 2014, Gan et al 2015, Ma et al 2015).

While this study projects decreases in summer discharge in the '2080s' (figure 5), an earlier study for headwater catchments of the Aksu, Hotan, and Yarkand projected only little changes (±5%) for summer discharge in the period 2081–2100 compared to 1981–2000 (Liu et al 2013). These differences are likely due to the static (and thus overestimated) glacier areas in the model of Liu et al (2013), which lead to an overestimation of glacier melt.

Estimates of the timing of peak melt vary by catchment and climate projection. Averaged over entire Central Asia, Bliss et al (2014) found a steady decline in glacier runoff over the 21st century, indicating that many glaciers have already passed peak melt. For the Chon-Kemin catchment in the northern Tien Shan, peak melt is projected to occur in the 2020s under warm climate scenarios but glacier runoff may also remain constant until 2100 under cooler climate scenarios (Sorg et al 2014). Only small decreases in glacier runoff until 2050 were projected for the Tanimas Basin in the Pamir (Hagg et al 2013).

For estimating water availability in the mainstream of the Tarim River, one needs to consider changes in land and water management as well as changes in other tributaries, which provide only a small contribution to the Tarim today but possibly a higher contribution in the future. Due to lacking meteorological and glacier mass balance data in the Hotan and Yarkand Basins, glacio-hydrologic simulations and climate change impact analyses in these basins are, however, highly uncertain.

In agreement with many other impact studies, our evaluation of uncertainties showed that the largest uncertainties stem from the climate models, while emission scenarios become a further important uncertainty source toward the end of the 21st century (Wilby and Harris 2006, Kay et al 2009, Prudhomme and Davies 2009, Ott et al 2013, Addor et al 2014). However, in data scarce mountain regions, uncertainties of hydrological models are large, and parameter uncertainties of hydrological models can become the largest source of uncertainties (Finger et al 2012, Ragettli et al 2013). The availability of glacier mass balance data in our study likely contributed to relatively small uncertainties from the hydrological parameters. In our case, parameter uncertainties were derived from different solutions of the multiobjective calibration. One would expect more diverse parameter sets if they had been derived by independent optimisation runs or Monte Carlo simulations and further parameter sets with a lower model performance would have been accepted. However, the simulated glacier mass changes (−0.83 to −0.51 m w.e. a⁻¹ and −0.85 to −0.25 m w.e. a⁻¹) cover large parts of the range given by the estimate from the observations (−0.85 to −0.19 m w.e. a⁻¹) and thus indicate that the parameterizations cover a good range of uncertainties. It is to note that uncertainties of the hydrological model are larger than the analysed parameter uncertainties due to uncertainties in model input and structure, as discussed below.

Further uncertainties in the climate change signal result from the coarse resolution of the climate models, which only partly resolve the complex topography of this region, and from the assumption of the stability of the bias correction under a changed climate. With respect to the glacio-hydrological modelling, further uncertainties stem from uncertainties in the calibration and input data, as well as the model structure. For example, while for the Sari-Djaz Basin, the glacier mass balance estimate used for model calibration largely overlaps with a more recent geodetic estimate derived specifically for this catchment (Pieczonka and Bolch 2015), an estimate specifically for the Kakshaal catchment is still missing. The applied ice-thickness estimates were based on a distributed ice-thickness model. While this type of thickness estimation method is seen as more reliable than approaches based on volume-area scaling, which tend to overestimate ice volume for large and steep glaciers (Frey et al 2014), uncertainties are still estimated to be in the range of 30% and contribute to uncertainties in the glacier and runoff evolution (Gabbi et al 2012, Huss et al 2014). Uncertainties in the model structure result from simplifications of the described processes and the assumption of stable model parameters under a changed climate. For example, the application of a temperature-index approach implicates that influences on the energy balance by changes in radiation or the wind field cannot be considered. However, comprehensive model testing showed that the model represents observed discharge trends and glacier mass loss in the past, which increases our confidence in the model. While the average influence of debris cover on glacier melt is implicitly taken into account, influences on glacier retreat by possible changes of the debris-cover extent cannot be considered by the current model. Our simulations did also not consider feedback effects of the vegetation, such as changes in the length of the growing period, the vegetation types and the water use efficiency. These aspects should be addressed by future research.