Meteorological data pre-processing

Based on the DEM, the glacierized region of the YRSR is divided, at intervals of 100 m, into a set of elevation bands. Each elevation band is assumed to exhibit different homogeneous hydrological characteristics, and the temperature and precipitation time series are interpolated according to its mean elevation. The interpolation is based on an altitude-dependent regression of the observations at stations located in or near the studied glacierized catchments. For the temperature time series, a constant lapse rate (0.50°C (100 m)⁻¹) is applied to the temperature series measured at the closest meteorological station (Reference Ding, Li, Liu and YuDing and others, 1992).

The precipitation increases linearly with altitude in the glacierized regions of the YRSR (Reference Ding, Li, Liu and YuDing and others, 1992; Reference YaoYao, 2002). The precipitation gradient is set to 11.2% (100 m)⁻¹ (Reference ZhangZhang and others, 1997; Reference YaoYao, 2002). Note that the spatial variability of precipitation is high in the glacierized regions of high mountainous catchments, and good spatial interpolation of precipitation data series is therefore difficult to obtain and the prevailing extreme conditions imply significant measurement uncertainty. In addition, the value of the precipitation gradient is difficult to estimate in the glacierized regions of the YRSR where precipitation data are scarce. This difficulty in spatial interpolation of precipitation is a significant source of modelling uncertainty for glacier runoff simulation. Nevertheless, glaciers represent the most important water storage reservoir in the glacierized regions of the YRSR, and for glacier runoff simulation any under- or overestimation of the area-average precipitation can be compensated by simulated ice melt.

Glacier accumulation and melt

For each elevation band, the aggregation state of precipitation is determined according to a simple temperature threshold:

where P _tot is the total precipitation on a given day, and P _snow and P _liq are solid and liquid precipitation. T is the mean daily air temperature and T ₀ is the threshold temperature.

Correct estimation of the aggregation state of precipitation is essential for the modelling of hydrological processes in high mountain catchments. Although a modelling approach based on a simple threshold function (Equation (1)) does not reflect the observed phenomenon, the results of previous studies have shown that the use of a simple temperature threshold is acceptable for the modelling of hydrological processes in the glacierized regions of high mountain catchments (WMO, 1986; Reference Ye, Yang, Ding, Han and ToshioYe and others, 2004; Reference Zhang, Shi, Li, Wang and DaiZhang and others, 2006c; Reference Hagg, Braun, Kuhn and NesgaardHagg and others, 2007). Fortunately, detailed data of precipitation and simultaneous meteorological conditions are available for the four national meteorological stations in the YRSR for the period 1961–79 (Fig. 1). Based on these data, a relationship between the probability of solid precipitation and temperature is obtained (Fig. 2). As indicated in Figure 2, precipitation is almost always solid when the air temperature is below 2.0°C. Above 2.0°C the probability of liquid precipitation significantly increases. This finding is in agreement with previous studies (e.g. Reference Ye, Yang, Ding, Han and ToshioYe and others, 2004; Reference Yang, Yang, Ding, Chen and LiuYang, 2006; Reference Zhang, Shi, Li, Wang and DaiZhang and others, 2006c). We therefore set the value of the threshold parameter to 2.0°C. A mixture of snow and rain is assumed for a transition zone ranging from 1 K above to 1 K below the threshold temperature. Within this temperature range, the snow and rain percentages of total precipitation are obtained from linear interpolation. Redistribution of the original snowfall by wind transport or avalanches is not considered.

Fig. 2. Relationship between the probability of solid precipitation and temperature in the YRSR.

Glacier melt, M, in the ith elevation band is computed using a degree-day approach, which rests upon a claimed relationship between snow- or ice melt and air temperature, usually expressed in the form of accumulated positive temperatures. Recent work shows that the degree-day approach has been justified more on physical grounds than previously assumed (Reference OhmuraOhmura, 2001). Although involving a simplification of complex processes that are more properly evaluated by the energy balance of the glacier surface, the degree-day approach often matches the performance of energy-balance models on the catchment scale (WMO, 1986). Reference Braithwaite and RaperBraithwaite and Raper (2002) show that there is excellent consistency between the results of energy-balance models and the degree-day model. The degree-day approach is given by

where i is an index for the elevation band, DDF is the degree-day factor, different for snow and ice surfaces, and T _m is the threshold temperature for melting, which is set to 0°C.

Glacier runoff

For the entire YRSR, glacier runoff, Q, on a given day is computed as

where s(i) is the area of the ith elevation band, derived from the DEM, and f is the refreezing rate.

Note that since significant amounts of meltwater have been captured as superimposed ice (Reference Fujita, Seko, Ageta, Pu and YaoFujita and others, 1996), refreezing of meltwater must be taken into account in simulating glacier runoff in the YRSR where cold-type glaciers are dominant (Reference HuangHuang, 1990). Recent work has suggested that approximately 20% of meltwater is preserved at the snow–ice interface due to refreezing processes in the mid- to northern Tibetan plateau (Reference Fujita, Ohta and AgetaFujita and others, 2007). Therefore, herein the refreezing rate f is set to 20%.

On the other hand, evaporation is of limited importance in the water balance of the glacierized regions in high mountainous catchments (Reference BraunBraun and others, 1994; Reference ZhangZhang and others, 1997). Previous work has suggested that the estimated annual evaporation from the glacier surface accounts for <6.0% of the total river runoff in the Dongkemadi glacier catchment of the YRSR (Reference ZhangZhang and others, 1997; Reference YaoYao, 2002). In addition, because the variation of glacier runoff primarily depends on changes of glacier area and climatic conditions, the influence of evaporation from the glacier surface on glacier runoff is small (Reference Ye, Chen and ShiYe and others, 1997). It is difficult to estimate evaporation from the glacier surface because of the sparsity of meteorological data in the glacierized regions of the YRSR, and it is not taken into account in the present analyses.

Model performance assessment

The model has two parameters (DDF_snow and DDF_ice) to be determined in addition to the above-mentioned fixed model parameters (Table 3). Computations of DDF_ice and DDF_snow are based on the available ablation datasets on the north slope of the Tanggula mountains over different periods (Reference Fujita and AgetaFujita and Ageta, 2000; Reference Yang, Yang, Ding, Chen and LiuYang, 2006). The spatial variability of degree-day factors in western China is analyzed from snow- and ice-melt data and meteorological data from different glaciers (Reference Zhang, Shi, Li, Wang and DaiZhang and others, 2006b). It is found that regional patterns of degree-day factor are detectable and that they vary very little across the Tibetan Plateau (Reference Zhang, Shi, Li, Wang and DaiZhang and others, 2006b). We therefore assume fixed values (Table 3) in the model for these parameters.

Glacier runoff in the glacierized regions of the YRSR is computed from meltwater and rainwater over the period 1961–2000, on the basis of the parameters shown in Table 3. Because few observational data exist in the glacierized regions of the YRSR, it is difficult to assess the performance of the model. The modelling performance is therefore assessed using the following processes.

We use the observed mass balance and ELA of Xiao Dongkemadi glacier to assess the performance of the model. Xiao Dongkemadi glacier is located in the Dongkemadi glacier catchment (Fig. 1), where the mass balance has been measured at a representative set of points in the accumulation area and the ablation area by the stake method once or twice a year since 1988 (Reference Fujita and AgetaFujita and Ageta, 2000). The annual mass balance at a given point on a glacier is defined as the sum of water accumulation in the form of snow and ice minus the corresponding ablation over the whole year (Reference PatersonPaterson, 1994). The above modelling approach enables the estimation of the annual mass balance based on the hydrological simulation outputs. For each elevation band, the mean annual mass balance is calculated based on the simulated snow accumulation and the simulated snow- and ice melt, which is obtained as

where b _a,i is the annual mass balance of the ith elevation band, t ₀ is the start date of the measurement year (here 1 October) and t ₁ is the end of the measurement year (30 September the following year). The annual mass balance of the entire glacier is estimated as the area-weighted sum of the mass balance of all elevation bands, which is given as

where B _a is the simulated total annual mass balance of the glacier and s _g is the area of the entire glacier.

The annual mass balance and ELA of Xiao Dongkemadi glacier were calculated for the period 1988/89–1997/98 using the above approach. As indicated in Figure 3, simulations are in good agreement with observations. Correlation of simulations with observations gives coefficients of 0.93 and 0.88 for the mass balance and ELA of Xiao Dongkemadi glacier respectively, significant at the <1% level. We submit that the processes of glacier ablation and accumulation are sufficiently well simulated by the above modelling approach, which considers only precipitation and temperature as underlying driving forces, and that the applied spatial interpolation of the driving forces can be assumed representative of the glacierized region.

Fig. 3. (a) Measured (squares) and modelled (circles) mean annual mass balance, and (b) the ELA of Xiao Dongkemadi glacier.

The total ice-volume change in the YRSR for the period 1968–2002 is used to assess the overall performance of the above modelling approach. The volume is estimated using a modified equation suggested by Reference Liu, Sun, Shen and LiLiu and others (2003):

where V is the ice volume and S is the total glacier area.

This empirical equation is derived from ice-penetrating radar thickness measurements in different glaciers of western China (Reference Liu, Sun, Shen and LiLiu and others, 2003), and has been widely applied to estimate glacier volume in different regions of western China (e.g. Reference LiuLiu and others, 2006a,Reference Liub). The recent glacier changes of the YRSR for the period 1968–2002 were monitored using topographical maps based on aerial photographs acquired during the period 1968–71 and Landsat Enhanced Thematic Mapper Plus (ETM+) images obtained for the period 1999–2002, which indicated that the total glacier area in the YRSR had decreased by approximately 5.3% (Xu and others, in press). On the basis of YRSR glacier area change, the total decrease in YRSR ice volume is estimated from Equation (6) to be ∼454.4 × 10⁸ m³ w.e. (assuming an ice density of 900 kg m⁻³). That is, for the period 1968–2002, the water resource in the YRSR increased by approximately 454.4 × 10⁸ m³ because of the loss of ice mass. Note that the loss of ice volume includes glacier melt and evaporation. Previous studies have suggested that evaporation provides a relatively small contribution to the glacier surface mass loss and is of limited importance to the water balance of the YRSR glacierized regions (Reference ZhangZhang and others, 1997; Reference YaoYao, 2002). Furthermore, measurements in about 20 evaporation pans showed a gradual decrease in evaporation between 1956 and 2000 (Reference Yang, Ding and ChenYang and others, 2007). The evaporation contribution is not separately accounted for within the model. Simulated glacier runoff using the model is approximately 422.0 × 10⁸ m³ over the entire YRSR for the study period 1968–2002. This agrees closely with the increase in water resources due to ice loss (454.4 × 10⁸ m³ vs 422.0 × 10⁸ m³), with a relative error of 7.0%.

Our modelling approach seems to perform well with respect to glacier runoff in the YRSR using the above model parameters, and our results can be regarded as acceptable for the source region where there are no observational glacier runoff data.

Influence of glacier runoff during recent decades

Making use of precipitation and temperature, YRSR glacier runoff is estimated for 1961–2000. For this period, glacier runoff accounts for 11.0% of the total river runoff. As shown in Figure 4b, glacier runoff for the YRSR displayed a generally increasing trend during recent decades, becoming more pronounced in the 1990s. This is associated with an accelerated warming in the 1990s in the YRSR, where the rate of warming estimated from the Geladandong mountain ice core is 0.5°C (10 a)⁻¹ since the 1970s, increasing to 1.1°C (10 a)⁻¹ during the 1990s (Reference KangKang and others, 2007). Various studies have suggested that the response of glacier runoff to climate warming is immediate in the initial phase (Reference Braun, Weber and SchulzBraun and others, 2000; Reference Jansson, Hock and SchneiderJansson and others, 2003). With the rapid warming in the 1990s, additional water is released from glacier storage due to intensified melting, resulting in a significant increase in glacier runoff.

Fig. 4. Variation in annual river runoff (a) and glacier runoff (b) of the YRSR for the climatic normal period 1961–2000. Dashed lines are regression lines indicating trend.

Figure 4 shows that, while there was high interannual variability, YRSR river runoff decreased over the period 1961–2000, especially in the 1990s. During this period, river runoff decreased by 13.9%. Glacier runoff increased by 15.2%. As indicated by the decadal anomalies (Fig. 5), river runoff was more plentiful for the periods 1961–70 and 1981–90, while glacier runoff indicated negative anomalies. River runoff in the 1990s reached a minimum, whereas glacier runoff reached a maximum (Fig. 5).

Fig. 5. Decadal river runoff and glacier runoff anomalies for the YRSR during recent decades.

Meteorological data indicate that precipitation in the YRSR significantly increased in winter and spring and dramatically decreased in summer during recent decades (Reference Zhang, Shi, Li, Wang and DaiZhang and others, 2006a; Reference Yang, Ding and ChenYang and others, 2007). With rapid warming and decreasing precipitation in summer, the alpine ecosystem, marshland and permafrost are generally degraded and glaciers show rapid retreat (Reference Cheng and ZhaoCheng and Zhao, 2000; Reference Ye, Kang, Chen and WangYe and others, 2006), resulting in significant decrease in river runoff (Reference Ding, Ye, Han, Shen and LiuDing and others, 2007). Human activity has increased since the end of the 1990s, significantly increasing water demand while the supply decreases, and resulting in intensive degradation of the YRSR alpine ecosystem. Such changes affect not only the source region, but also the middle and lower reaches of the YRB. Glacier runoff in the YRSR, however, rapidly increased during the 1990s, to contribute ∼17% of the river runoff. The impact of glacier runoff on river runoff has become stronger than ever before in the YRSR, especially in the 1990s when river runoff exhibited a significant decrease. To a large extent, a significant increase in glacier runoff in the 1990s may have reduced the impact on the YRSR ecological environment of the decrease in river runoff.

YRSR climate scenario

On the basis of ECHAM5/MPI-OM projected temperature and precipitation, possible temperature and precipitation changes in the YRB from 2001 to 2050 were analyzed (Reference Zeng, Su, Jiang and WangZeng and others, 2007). Compared to the observed temperature and precipitation for the period 1961–2000, the projected temperature and precipitation of ECHAM5/MPI-OM are in good agreement over the entire YRB using the Mann–Kendall test, especially in the upper reaches of the river basin where the performance of ECHAM5/MPI-OM is slightly better than in the middle to lower reaches (Reference Zeng, Su, Jiang and WangZeng and others, 2007). According to the results of Reference Zeng, Su, Jiang and WangZeng and others (2007), ECHAM5/MPI-OM is capable of projecting air temperature and precipitation in the YRB, especially for the upper reaches. Nevertheless, meteorological outputs from climate models are not directly applied for impact studies in general, because climate models are unable to represent local subgrid-scale features and dynamics (Reference Giorgi and HoughtonGiorgi and others, 2001). Biases can develop in both temperature and precipitation, especially precipitation which highly depends on the local orographic conditions. Note that such differences will strongly affect simulated performance since the degree-day approach is particularly sensitive to the seasonal distribution of temperature. Downscaling techniques therefore need to be applied to the climate model output (Reference WilbyWilby and others, 1998; Reference Giorgi and HoughtonGiorgi and others, 2001).

In this study, we apply a simple statistical downscaling method, a non-linear regression method, which involves establishing a non-linear relationship between a large-scale average surface variable and the local-scale variable (e.g. Reference Kim, Chang, Baker, Wilks and GatesKim and others, 1984; Reference Wigley, Jones, Briffa and SmithWigley and others, 1990; Reference WilbyWilby and others, 1998). We can build a statistical relationship between the observed and projected climate variables based on the observations of the four meteorological stations in the YRSR. Consequently, future temperature and precipitation time series in A2 and B1 emission scenarios were calculated using

where y is a future temperature or precipitation time series, x is a projected temperature and precipitation time series (using ECHAM5/MPI-OM), and a ₁, a ₂, a ₃ and b are parameters given in Table 4. Correlations for precipitation are lower (as expected) than those for temperature (Table 4). This is because the spatial variability of precipitation in the catchment is high, and the reliability of the projected precipitation depends strongly on the local orographic conditions. On the other hand, climate models tend to underestimate large amounts of precipitation and over-estimate small amounts (Reference WilbyWilby and others, 1998; Reference Giorgi and HoughtonGiorgi and others, 2001). The histogram of the residuals indicates that the data are reasonably symmetric, there do not appear to be significant outliers in the tails, and it seems reasonable to assume that the data are from an approximately normal distribution (Fig. 6). It can be seen that the validity of the above downscaling method is acceptable for the present study. According to the above statistical relationship, future temperature and precipitation time series derived from downscaling from ECHAM5/MPI-OM agree closely with the observational data. Correlation coefficients are 0.94 and 0.52, respectively.

Fig. 6. Histogram of the residuals of (a, b) temperature and (c, d) precipitation in (a, c) A2 and (b, d) B1 emission scenarios.

As indicated in Table 5, there is a significant increase in temperature in the YRSR, as well as the YRB, in A2 and B1 emission scenarios for the period 2001–50 compared to 1961–90. Surface temperature in the YRSR by 2050 is likely to increase by approximately 1.44°C (A2) and 1.53°C (B1) compared to the climatic normal period (1961–90), especially for spring and winter (Fig. 7a). Precipitation in the source region, however, is likely to reduce by approximately 12.2% and 10.3% (Table 5) in A2 and B1 emission scenarios respectively, compared to the climatic normal period.

Fig. 7. Variations of temperature projected by ECHAM5/MPI-OM in A2 and B1 emission scenarios in (a, b) the YRSR and (c) the Dongkemadi glacier catchment (DGC) for 2001–50 relative to 1960–90.

Future trend of glacier runoff in the Dongkemadi glacier catchment

The Dongkemadi glacier catchment, situated in the south of the YRSR, is chosen as a representative glacierized region for this study (Fig. 1). Based on the meteorological observations in this catchment, we find an exponential relationship between daily glacier runoff, R _g, and mean daily air temperature, T _g, at 5600 ma.s.l. (Reference ZhangZhang and others, 1997):

An extensive discussion of this relationship can be found in Reference ZhangZhang and others (1997). As indicated in Figure 7c, future temperature in the Dongkemadi glacier catchment is likely to rise by approximately 1.57°C (A2) and 1.66°C (B1) during 2001–50, with respect to 1961–90, especially in spring and winter. Glacier runoff in the catchment, estimated using Equation (8), is very likely to significantly increase in both A2 and B1 emission scenarios (Fig. 8). Compared to the climatic normal period, glacier runoff by 2050 from this catchment is enhanced by 24.4% (A2) and 25.4% (B1) (Table 6). By 2050, a significant increase in glacier runoff occurs during April–June, with less increase during July–September compared to the climatic normal period (Fig. 8b). In the A2 emission scenario, the depth of glacier runoff varies from 0.06 to 7.8 mm, and from 0.07 to 8.1 mm in the B1 emission scenario. Contemporary comparison shows that the maximum glacier runoff in the B1 emission scenario occurs earlier than in the A2 emission scenario (Fig. 8a). On the other hand, climate warming in this catchment is expected to prolong the glacier melt season (Fig. 8a), thus reducing seasonal runoff concentration, i.e. the ratio of maximum monthly to mean monthly glacier runoff is expected to decrease from 3.5 during the period 1961–90 to 2.7 during the coming half-century in A2 and B1 emission scenarios.

Fig. 8. (a) Daily and (b) monthly variations of glacier runoff for the period 1961–90, in A2 and B1 emission scenarios in the Dongkemadi glacier catchment.

Future trend of glacier runoff in the YRSR

Based on the climate-change scenarios, glacier runoff in the YRSR was estimated for the period 2001–50 using Equations (1–3). The evolution of annual glacier runoff shows a uniform pattern in A2 and B1 emission scenarios for the period 2001–50 (Fig. 9a and b). Discharge gradually increases in the near future due to enhanced melting. Compared to the climatic normal period 1961–90, glacier runoff by 2050 is expected to increase by 29.2% (A2) and 29.8% (B1) (Table 6), i.e. almost twice the increase during the period 1991–2000 when the annual runoff was 15.2% higher than during the climatic normal period.

Fig. 9. (a) Evolution of annual glacier runoff in A2 and B1 emission scenarios (2001–50) and (b) annual cycle of glacier runoff (averaged over the period 1960–90), in A2 and B1 emission scenarios in the YRSR.

We have evaluated the annual glacier runoff cycle during the period 1961–90 and in A2 and B1 emission scenarios (Fig. 9c). By 2050, the results of the two emission scenarios show a significant decrease in the July–August runoff to much lower values than during 1961–90: −23.2% (A2) and −19.1% (B1). In contrast, there is a significant increase in the April–June runoff during an earlier onset of the melt season: +121.1% (A2) and +136.0% (B1). Due to enhanced temperature, there is a significant increase in the length of the melt season (Fig. 9c), reducing the seasonal runoff concentration, i.e. the ratio of maximum monthly to mean monthly glacier runoff decreases from 3.0 during the period 1961–90 to 2.2 during 2001–50.

The response of glacier runoff to climate warming is a matter of timescale. Although glacier runoff significantly increases through enhanced meltwater production and more efficient water transport through glaciers in recent decades, it will decrease due to some response variables changing sign at a later stage when enhanced melt rates have caused glacier volume to decrease significantly under a continuous warming in the coming decades (Reference Braun, Weber and SchulzBraun and others, 2000; Reference Jansson, Hock and SchneiderJansson and others, 2003). However, glacier runoff in the YRSR by 2050 still shows a significant increase as indicated in Figure 9a and b. It can be seen that the contribution of glacier runoff to river runoff of the YRSR is important during the first half of the 21st century.