Introduction

The Tibetan Plateau, the roof of the world, is home to more than 40,000 glaciers (Arendt et al., 2015). Glaciers on the Tibetan Plateau feed many important rivers (e.g., Brahmaputra, Ganges, Yellow, Yangtze, Indus, Salween, and Mekong), which are vital for agricultural irrigation and hence directly affect millions of people's lives. Under the impact of climatic warming, many glaciers in the Tibetan Plateau have experienced a rapid retreat in recent decades, although in some regions a minority of glaciers have advanced mainly due to the combined effects of precipitation and air temperature (e.g., Yao et al., 2012; Li et al., 2010;Yao et al., 2004; Bolch et al., 2012). Glacier loss has increased concern for water resources (Immerzeel et al., 2010; Lutz et al., 2014), but rates of mass loss vary and are greatest for the maritime glaciers and less in more continental regions (Rupper and Roe, 2008; Zhao et al., 2014b, 2016). A few studies have been done about the variations of equilibriumline altitude (ELA) and length for typical glaciers on Tibetan Plateau under climate change using a statistical model (e.g., Wang et al., 2010; Wang et al., 1996). Quantification of future Tibetan Plateau glacier response to climate change is also important given the lack of viable alternative water resources, but such quantification rarely has been attempted (e.g., Zhao et al., 2014a) using physically based simulations of typical glaciers.

To quantify changes in glacier volume, both the surface mass balance (SMB) determined by local accumulation and ablation rates on the glacier surface, and ice dynamics, which transfers mass from accumulation to ablation zones, need to be measured or estimated. However, in situ SMB data are only available for 15 glaciers (e.g., Yao et al., 2012; Li et al., 2010) on the Tibetan Plateau due to the inaccessibility and high elevation of most glaciers. Satellite surveys of glacier change have been used to map surface elevation changes across the region, but the data are limited to specific ground tracks and can only give information over a short period (e.g., Gardner et al., 2013). Repeated stereo imagery can also give estimates of elevation change over time with relatively low accuracy (e.g., Pieczonka et al., 2013). None of these methods is very useful in determining the SMB with respect to elevation over a glacier because they provide integrated estimates of regional mass budgets. SMB may be estimated using tailored surface energy-balance models (e.g., Jiang et al., 2010; Li et al., 2015). However, the use of a full surface energy-balance model is difficult because it requires a large number of input parameters that are neither observed in meteorological stations nor accurately reproduced by climate models. The degree-day model is a simple and popular approach (Braithwaite, 1984) to parameterize the SMB changes given that it only requires air temperature and precipitation and a linear coefficient (a so-called degree-day factor) between ablation and the sum of positive temperature.

Numerical ice-flow models of varying complexity have been applied on alpine glaciers (e.g., Gudmundsson, 1999; Le Meur and Vincent, 2003; Zwinger and Moore, 2009; Réveiller et al., 2015), but rarely on Tibetan glaciers. Simulations of a typical glacier (Gurenhekou) on the southeastern Tibetan Plateau suggested that ice dynamics was less important to future glacier evolution than the SMB (Zhao et al., 2014a) due to the slow ice-flow speeds. This implies that simplified parameterizations of dynamics based, for example, on volume/ area scaling may be used instead of ice dynamics models (e.g., Moore et al., 2013; Shi et al., 2016), but verification on specific glaciers across the region has not yet been done.

In this paper, we combine an SMB parameterization using the degree-day model with the full-Stokes ice-flow model to simulate the response of a representative glacier in the central Tibetan Plateau to different climate forcing from 2013 to 2050. The study region and glaciological measurements are described in the next section. The SMB parameterization, the model set-up, and climate forcing are then described, with the simulations and projections shown and discussed in the final sections.

Study Area and Observational Data

Qiangtang No. 1 (33.29°N, 88.70°E) is a small (2.4 km²) glacier located in the extremely continental Shuanghu, eastern Qiangtang Plateau of the central Tibetan Plateau (Fig. 1, part a). Qiangtang Plateau is surrounded by Tanggula, Nyainqentanglha, and Gandisê Ranges. It has a mean elevation of around 4800 m, and its climate is influenced by both southwest monsoon and westerly winds. Moisture feeding the glaciers in eastern Qiangtang Plateau comes mainly from the southwest monsoon (Chen et al., 2012), hence both accumulation and ablation occur overwhelmingly during the summer season. Data from the nearest meteorological station to the glacier, Bange (31.38°N, 90.0°E, 4700 m a.s.l., 245 km away from the glacier), shows that over the past three decades, the mean annual air temperature increased by 0.035 °C a^-1, while the annual precipitation increased slightly by 0.8 mm a ^-1 (Fig. 2).

Qiangtang No. 1 Glacier ranges between 5500 and 6050 m a.s.l. (Zhu et al., 2014) and is about 2 km in length. This glacier has been heavily studied on the Qiangtang Plateau due to its accessibility, and it is quite representative in terms of size, elevation and aspect of glaciers in the region. Differential GPS and ground-penetrating radar (GPR) surveys of the glacier were made in December 2011 and June 2013 (Zhu et al., 2014; Fig. 1, part b). The real time kinematic (RTK) GPS measurement was performed using a pair of E650 type dGPS units manufactured by Unistrong Company. The two dGPS units were connected by radio signal, enabling an RTK survey to be carried out with centimeter-level accuracy measurements. The GPR measurements were performed with a pulse EKKO 100™ system from Canadian Sensors & Software. We used 100 MHz frequency antennas, with a 4 m separation between the two antennas during the measurements. The GPS precision of single spot is 0.05 m and the ice thickness data has an accuracy of around 1–2 m. The maximum ice thickness observed was 132 m, 454 m down-glacier from the ice divide (Fig. 1, part c). An ice core was drilled to the bedrock (109 m) at the ice divide (33.297°N, 88.695°E, 5853 m a.s.l. Fig. 1, part b) on 4 May 2014, and an ice temperature profile was measured in the borehole on the same day (Fig. 3) by a KN2014-000 platinum resistance thermometer sensor with a precision of 0.01 °C. The borehole ice temperature was also measured several other times over the following year of 2015. From Figure 3, we can clearly see a seasonal air temperature impact in the upper 20 m. Our observed borehole ice temperature profile suggests a geothermal heat flux of 127 mW m^-2. This value is consistent with observations (100– 150 mW m^-2) across the central Tibetan Plateau (Hu et al., 2000), but higher than the value (79 mW m^-2) derived from Naimonanyi glacier (734 km from Qiangtang glacier) in the southwestern Tibetan Plateau (Zhao et al., 2014a).

FIGURE 1.

(a) The locations and image from Google Earth of Qiangtang No.1 Glacier, automatic weather station 1 (AWS1) (red dot), and AWS2. The glacier outline is marked by black curve. The inset map shows the Tibetan Plateau, the Google image location (red box, with sides enlarged by a factor of three for visibility), and Bange station (blue dot), (b) Glacier outline and surface elevation (m a.s.l.) interpolated from 2013 ground penetrating radar (GPR) measurement (blue curves), and global positioning system (GPS) profiles from 2011 and 2013 (black dots) with the glacier terminus mapped in 2013. The ice divide and ice drilling site is marked by a big green dot. (c) Ice thickness (m) interpolated from the radar and differential GPS data.

FIGURE 2.

(a) Annual mean 2 m air temperature and (b) annual precipitation and their trends at Bange station from 1980 to 2012.

FIGURE 3.

The measured borehole temperature (red) and that simulated by the dynamic flow model at the borehole (blue).

Daily air temperatures were measured by an automatic weather station (AWS1; Fig. 1, part a), which is about 10 m away from the ice coring site from November 2012 to November 2014. Daily air temperature and precipitation data were also recorded about 14 km away from the glacier (AWS2, 33.28°N, 88.84°E, 4974 m a.s.l; Fig. 1, part a) from October 2011 to September 2015.

Surface Mass Balance Parameterization

There are no SMB data available for Qiangtang No. 1 Glacier. Thus we parameterize the glacier surface accumulation and ablation as functions of local air temperature and precipitation, which is then validated against the measured surface elevation changes between December 2011 and June 2013 (Fig. 4).

FIGURE 4.

Measured surface elevation change (dots and crosses) versus modeled surface mass balance (SMB) (curve) from December 2011 to June 2013. Red crosses show data below 5760 m, close to the equilibrium line altitude (ELA).

We use the precipitation data from AWS2 between October 2011 and September 2015 to estimate the annual precipitation on the glacier. A precipitation lapse rate of around 8.6 mm (100 m)^-1 is estimated using the precipitation records at AWS2 and Bange Station.

We employ the positive degree-day (PDD) model (Braithwaite, 1984) to estimate the ablations at different elevations. The annual ablation rate α is calculated by multiplying the sum of positive daily mean air temperatures in the entire summer (June, July, and August, or JJA) by a suitable degree-day factor (DDF). From the summer (JJA) daily air temperatures in the years 2013 and 2014, we find a lapse rate of 0.81 °C (100 m)^-1 between AWS1 and AWS2, close to the lapse rate, 0.84 °C (100 m)^-1, between AWS2 and the Bange station.

There is a large regional variability in DDF across the Tibetan Plateau, with reported values ranging from 3 to 15 mm day^-1 °C^-1 (e.g., Kayastha et al., 2003; Zhang et al., 2006). Zhang et al. (2006) analyzed the spatial variability of DDFs obtained from 15 monitored glaciers in different regions of western China, and suggested a value of 7 mm day^-1 °C^-1 for the region around Qiangtang No.1 Glacier. Therefore, we use a DDF of 7 mm day^-1 °C^-1 in our PDD model.

We compared the observed surface elevation changes between December 2011 and June 2013 with the calculated SMB values (Fig. 4). In general, they are in good agreement below the equilibrium-line altitude (ELA), but above the ELA the points of surface elevation change show a scattering pattern and appear to be less than the SMB values we estimate.

Summer daily air temperature data at AWS1 and summer daily air temperature and precipitation data at AWS2 were used to model annual SMB for the period 2012–2015 (Fig. 5, part a). To extend the record to 1980, we use daily temperature and precipitation data from the Bange station (Fig. 5, part b). The simulated SMB is repeatable: SMB increases linearly with a gradient of (3.6 ± 0.1) mm m^-1 below the ELA, and negative exponentially above the ELA to 6000 m with an asymptotic SMB of 0.450+0.010 m yr^-1. Therefore, following Zhao et al. (2014a), the SMB can be parameterized as a function of elevation and ELA (where SMB is zero) :

where β = 3.6 ± 0.1 mm m^-1 and (φ = 0.450 ± 0.010 m yr^-1 represent SMB gradient and the upper limit of SMB, respectively. The ELA time series (Fig. 6) derived from the modeled SMBs over the period 1980–2012, together with time series of summer mean air temperature from Bange station over the same period, yielded a modeled ELA sensitivity to summer temperature of 130 ± 15m°C^-1.

FIGURE 5.

(a) The modeled SMB in the year 2012 (blue), 2013 (red), 2014 (black), and 2015 (green), (b) The modeled SMB from 1980 to 2011. The 1980s, 1990s, 2000s, and 2010s are plotted in blue, red, green, and yellow, respectively.

Rupper and Roe (2008) assessed the sensitivity of ELAs to changes in regional climate for Central Asia and found that ELAs in melt-dominated regions (as is Qiangtang No. 1 Glacier) were most sensitive to interannual variability in air temperature, while ELAs in sublimation-dominated regions were most sensitive to interannual variability in precipitation. Zhao et al. (2016) considered ELA parameterizations with and without ELA sensitivity to precipitation for glaciers in high mountain Asia, and found a 12% difference in the glacier mass loss. Therefore, if we assume that changes of ELA for Qiangtang No. 1 are mainly controlled by summer temperature variation, then the ELA in the k _th year after a starting year (2013) can be estimated by

where Δ T _k is the net change of mean summer temperature between 2013 and the kth year and α is the sensitivity of ELA to temperature changes, and α= 130 ± 15 m °C^-1.

FIGURE 6.

The ELA time series during 1980–2012 resulting from modeled SMB.

Modeling Approach

Qiangtang No. 1 Glacier flow dynamics was modeled using the thermomechanically coupled full-Stokes model, Elmer/Ice, an open source finite element code that has been successfully applied to several mountain glacier flow modeling studies (e.g., Gagliardini et al., 2007; Zwinger et al., 2007; Zwinger and Moore, 2009; Jay-Allemand et al., 2011; Adhikari and Marshall, 2013; Gagliardini et al., 2013; Zhao et al., 2014a; Réveillet et al., 2015). We use it here because previous experience on a similar glacier on the Tibetan Plateau showed that the full stress components were needed to capture the glacier dynamics (Zhao et al., 2014a).

Field Equations

The thermomechanically coupled full-Stokes model consists of the following equations:

Equation 3 expresses the conservation of mass in an incompressible fluid and is called the incompressibility condition. The ice density ρ is assumed to be a constant 910 kg m^-3; v = (u, v, w) is the vector of ice flow velocity. Equation 4 expresses the conservation of momentum, where τ is the deviatoric stress tensor, p is the isotropic pressure, and g is the vector of the acceleration due to gravity and equal to (0, 0, -9.81) m s^-2 The deviatoric stress tensor, τ, is related to the deviatoric part of the strain rate tensor, D, via Glen's flow law

and the ice viscosity, η, is given by

where the flow rate factor, A(T), is derived from the Arrhenius law (Paterson, 1994); d _e is the second invariant of the strain-rate tensor, d _e = (0.5tr D ²)^1/2; and the Glen exponent n is taken to be 3. Equation 5 is the heat transfer equation where T is the ice temperature. The heat capacity c and heat conductivity κ are functions of the temperature (Paterson, 1994). The temperature in Equation 5 is constrained by the upper limit imposed by the pressure melting point.

FIGURE 7.

(a) Two June–July–August mean air temperature warming scenarios (with respect to the temperature in 2013) used for prognostic simulations during 2013–2050. The extrapolated historical warming trend (red curve) with averaged rate of 0.035 °C a^-1 and sequential annual anomalies observed at Bange station over the period 1975–2012, and projections from RegCM3 under the A1B scenario for Qiangtang No. 1 Glacier (blue curve), (b) The mesh on the 2-D footprint with resolution of 50 m.

Boundary Conditions

The borehole temperature profile shows a frozen ice-bedrock interface (Fig. 3). As an initial assumption, we assume that the whole glacier is “cold” and hence no sliding occurs at the ice-bedrock interface. We shall show that this assumption leads to consistent results in the Simulation Results section. We set a Neumann boundary condition at the ice/bedrock interface for the temperature simulation with a geothermal heat flux of 127 mW m^-2.

At the free surface, a vanishing deviatoric stress vector, τ, is imposed, and the atmospheric pressure and its change over the glacier surface are taken to be zero. We estimate the mean annual surface temperature at the borehole (5856 m) from the measured temperature (-9 °C; Fig. 3) at 10 m depth, and glacier surface temperatures using a lapse rate of -0.008 °C m^-1, which is the annual average daily temperature lapse rate found between AWS1 and AWS2,

where T is the ice surface temperature (°C) and z is the ice surface elevation (m a.s.l.).

The evolution of the free surface during the prognostic simulation 2013–2050 was described using a kinematic boundary condition:

where z is surface elevation; u, v, and w denote ice flow velocity components in the x, y, and z directions, respectively; and SMB is calculated as described in the Surface Mass Balance Parameterization section.

Simulation Results

Glacier Evolution to 2050

Using a fixed glacier geometry measured in 2013, we first performed a diagnostic simulation to build the steady-state temperature and velocity fields, which were then used as an initial condition for prognostic simulations during the period of 2013–2050 to predict glacier evolution driven by the parameterized SMB profiles (see details in Surface Mass Balance Parameterization section) under the two climate scenarios described in Boundary Conditions above. Here we use the means of parameters in SMB profiles (see Equations 1 and 2).

The simulated and measured ice temperature distributions in the borehole are shown in Figure 3 for the steady-state simulation. The modeled temperature fits well with the measurement at depths below 20 m where seasonal temperature fluctuations are damped. The modeled surface and vertical distribution of velocity are shown in Figure 8. The maximum speed (8.5 m a^-1) appears at around the middle of the glacier. The modeled basal temperature distribution of the glacier (Fig. 9, part a) shows that the whole glacier is cold (below pressure-melting-point), consistent with our initial no-sliding assumption.

Before running the prognostic simulations, we performed a “surface relaxation” experiment (Zwinger and Moore, 2009) over five years without SMB forcing to see if parts of the glacier exhibit unrealistic surface velocities that were likely caused by errors in the surface or bedrock DEMs. The maximum deviation of surface velocities is only 0.1 m a^-1. Therefore, this procedure did not detect systematic errors in the DEMs.

The steady-state solution of the diagnostic run is the starting point for the prognostic simulations from 2013 to 2050. The evolution of the glacier was forced by the time-changing SMB distribution over the whole glacier, and driven by both the RegCM3 projected scenario and the historical warming scenario (Fig. 7, part a).

The simulated changes in surface elevation over the period 2013–2050 and glacier margin outlines in 2050 under two scenarios are shown in Figure 10. The measured mean retreat speed of the glacier front between 2000 and 2011 was around 5 m a^-1. The modeled mean retreat speed is around 3.8–4.1 m a^-1 between 2013 and 2027. After 2027, the modeled terminus retreat speeding accelerates to around 6.6–7.3 m a^-1 (Fig. 11).

Changes of area and volume under both scenarios are shown in Figure 12. Because the warming in the RegCM3 scenario is stronger than the historical warming scenario, area and volume decrease more rapidly under the RegCM3 forcing than the historic warming trend. Under the RegCM3 projection, the modeled glacier area and volume in 2050 drop to 1.4 km and 0.069 km³—that is, around 82% and 70% of that in 2013, respectively. Under the extrapolated historical warming scenario, the modeled glacier area and volume in 2050 reduces to 1.5 km² and 0.08 km³—that is, around 89% and 81% of that in 2013, respectively. The mean volume loss rate increases from (0.035 ± 2.5) × 10^-4 km³ a^-1 for 2013–2015 to (∼7.5–10.7) × 10^-4 km³ a^-1 for 2040–2050. This increase can be compared with the calculated summer melt rate (6.9 ± 1.7) × 10^-4 km³ a^-1 below the ELA for 2013–2015.

FIGURE 8.

The modeled ice flow velocity for the diagnostic run (a) on the surface and (b) on the transect along the central line (the vertical coordinate is stretched by 3).

FIGURE 9.

Distribution of temperature relative to pressure-melting point (a) at the bedrock and (b) on the transect along the central line in the diagnostic run (the vertical coordinate is stretched by a factor of 3).

FIGURE 10.

Modeled ice surface elevation change (positive value for thinning) in 2050 for (a) RegCM3 scenario and (b) historical warming scenario. The black lines indicate glacier outlines in 2050.

Uncertainty in Glacier Evolution

Due to difficulties with in situ measurements, the GPR data on this glacier is relatively sparse and may cause some errors in the ice thickness interpolations and geometry generation. In particular, the lower part of glacier lacks measurements because it is steep and slippery, which made it dangerous to measure (Fig. 1, part b). To investigate at what size the uncertainty of geometry may influence the simulation results, we derive an alternative ice thickness distribution and generate a corresponding glacier geometry from surface topography and surface mass balance using the method described by Huss and Farinotti (2012). There are differences up to tens of meters between this alternative ice thickness and the one we derive from measurements. Then we perform glacier evolution simulations based on this alternative “Huss and Farinotti” geometry. The maximum speeds in steady-state simulation are reduced from about 8.5 m a^-1 to 2.0 m a^-1 with the Huss and Farinotti derived geometry. Simulated temperatures at the borehole in steady-state simulation are almost the same as with the measured geometry over the common depth interval, though the Huss and Farinotti ice thickness at the borehole is only about 100 m. Despite these dynamic differences, in the prognostic run, the volume changes based on two geometries are almost the same and the area changes are also quite similar. This implies that SMB is the dominant factor in glacier change rather than ice dynamics. However since the ice thickness in the two geometries differs appreciably, the glacier terminus retreat speeds also differ. The modeled terminus retreat speed is around 2.9– 3.4 m a^-1 based on “Huss and Farinotti” geometry, while 3.8–4.1 m a^-1 using the measured geometry between 2013 and 2027. After 2027, the modeled terminus retreat speed slows down to 1.9–2.2 m a^-1 with “Huss and Farinotti” geometry, but accelerates to around 6.6–7.3 m a^-1 using the measured geometry.

FIGURE 11.

Cumulative retreat distance of Qiangtang No.1 Glacier for the RegCM3 scenario (blue curve) and the historical warming scenario (red curve).

FIGURE 12.

(a) Area and (b) volume evolution over time of Qiangtang No.1 Glacier and its estimated uncertainty under the RegCM3 scenario (blue band), and under the historical warming scenario (red band); (c) volume loss rate over time of Qiangtang No.1 Glacier under RegCM3 scenario (blue curve) and the historical warming scenario (red curve).

SMB is the most important factor in control of this glacier's evolutions. There are three key parameters in our SMB parameterizations (see Equations 1 and 2) : the SMB altitude gradient (3.6 ± 0.1 mm m^-1), the maximum upper limit of SMB (0.450 ± 0.010 m yr^-1), and ELA sensitivity to summer mean temperature changes (130 ± 15 m °C^-1). We study the sensitivity of glacier area and volume change to these parameters. We find (1) the upper limit of SMB has negligible effect on the uncertainty of glacier change; (2) a change of SMB gradient by 1 standard error results in ±0.66% and ±1.22% change in volume projection under historical warming scenario and RegCM3 scenario, respectively; and (3) a change of ELA sensitivity to summer mean temperature by 1 standard error results in ±0.40% and [-3.16%, 2.50%] change in the volume projection under historical warming scenario and RegCM3 scenario, respectively. The sensitivity of area change to each parameter is similar to that of volume projection.

Monte Carlo sampling approach is a robust way to investigate the uncertainty of glacier volume or area projection caused by joint uncertainties of all three parameters in Equations 1 and 2. This would, however, be computationally prohibitively expensive. As a guide to the joint sensitivity, we make combinations of the best estimate ± 1 standard error, of pairs of parameters (producing 12 extra sets of simulations per climate scenario); these combinations would be relatively rare (p = 0.05) if the parameters were normal and independent variables, but we cannot assume that and so make no claim for a statistical confidence interval. We simply take the maximal interval of the resulting uncertainty ranges as defining area and volume uncertainty defined by SMB. The resulting uncertainty range is [-0.96%, 0.90%] and [-4.45%, 3.64%] for volume projection, and [-0.65%, 0.74%] and [-1.87%, 2.30%] for area projection (Fig. 12), under the historical warming and RegCM3 scenarios, respectively.

Discussions and Implications

We have applied a 3-D thermomechanically coupled full-Stokes model to a mountain glacier, Qiangtang No. 1 Glacier, on the Qiangtang Plateau of the central Tibetan Plateau. As with almost all glaciers on the Tibetan Plateau, there is no published SMB data for this glacier. To estimate the SMB for the past 30 years, we used a degree-day model based on air temperature and precipitation data at two AWSs and one meteorological station near the glacier. Then we use the modeled SMB characteristics (the SMB-altitude gradients and ELA sensitivity to air temperature) to parameterize the future changes of SMB. The modeled SMB profile (a function of elevation) between December 2011 and June 2013 fits reasonably well with the observed surface elevation lowering below the ELA, while above the ELA the elevation changes are in general spatially scattering, possibly due to firn compaction and snow redistribution by winds.

The lack of SMB observation is the main limitation for this study. More field surveys of SMB and ice flow velocity in the future will help to calibrate and validate the model. The uncertainty of the modeled SMB mainly comes from the accuracy of (a) the degree-day factor DDF and (b) the extrapolation of temperature and precipitation from the meteorological station or AWSs near the glacier. Though our model approach of the ELA change is parameterized only by the variation of air temperature, it also incorporates the covarying change of precipitation. In accordance with Zhao et al. (2016), we assume the glacier changes in the central Tibetan Plateau are mainly controlled by the summer air temperatures. For the Qiangtang No. 1 Glacier, we find an ELA sensitivity of 130 ± 15 m °C^-1, smaller than the one (172 m °C^-1) for Qiyi Glacier (39.23°N, 97.75°E) over the period 1958–2008 (Wang et al., 2010) and the one (321 m °C^-1) calculated from energy balance methods for central Asia (Rupper and Roe, 2008). Note that Zhao et al. (2014a) obtained similar ELA sensitivities (79 and 135 m °C^-1) for Gurenhekou Glacier in the southern Tibetan Plateau.

An ice core to bedrock was drilled near the top of the glacier in the year 2014. The measured ice temperatures point to a relatively high geothermal heat flux of 127 mW m^-2. The modeled steady state temperature profile at the borehole fits well with the measured one. The agreement between the measured and the modeled borehole temperature implies that the heat sources which are not considered (such as latent heat released by refreezing of meltwater in the firn) in the model can be ignored in this location and under present climate; of course, that may not be the case over the whole glacier over the coming century. Now, the whole glacier seems cold-based despite the warmest basal temperature being only -0.2 °C below pressure melting point, consistent with the absence of any subglacial outlet stream. The geothermal heat flux from the measured temperature profile is higher than the southeast of Tibetan Plateau (Hu et al., 2000), which leads to a relatively warm bed. We assumed a no-slip boundary condition for all the simulations, consistent with no basal temperate ice during the prognostic experiments, and the absence of a subglacial water outflow. But as the basal temperature is locally very close to the pressure melting point, sliding of the ice on its bed is possible in some locations. Surface velocities would be required to determine the basal friction distribution and hence deduce sliding. Trim lines visible for the lower part of the glacier show the glacier was thicker in the past, and parts of the bed may have been at the pressure melting point. Surface warming is unlikely to cause basal melting in the future, however, because the thermal diffusion time is likely longer than the survival of the glacier under the warming.

We simulated glacier evolution from 2013 to 2050 under two future warming projections—temperatures from the high-resolution RegCM3 under the A1B scenario, and the historical warming trend over the past three decades. Almost the whole glacier is subject to ablation in 2050, with losses of 11%–18% in area and 19%–30% in volume under the RegCM3 and historic trend forcings. The modeled retreat speed of glacier terminus before 2027 is 3.8–4.1 m a^-1, which is slightly less than the observed rates from 2000 to 2013. We suggest retreat will accelerate to 6.6–7.3 m a^-1 after the year 2027. Runoff by 2040–2050 increases under both warming scenarios by 0.00043–0.00075 km³ a^-1, which is comparable to the present rate of summer melt on the glacier. It is thus likely that downstream infrastructure would have to cope with volumes of water flow that may be twice as large as present over the next 50 years.