1. Introduction

Glaciers are one of the major components of the cryosphere. They are ancient rivers of compressed snow that creep through the landscape, shaping the planet’s surface. They tend to be located in wet and cold places. Wet implies a region not too far from a moisture source; cold implies high latitude and/or altitude (Reference OerlemansOerlemans, 2001). Glaciers are considered excellent climate-change indicators (e.g. Reference Hall and FagreHall and Fagre, 2003) because their dimensions change in response to several climatic variables characterizing the local environment (Reference Kadota and AgetaKadota and Ageta, 1992). Currently, the Earth is at its warmest since record keeping began around AD 1850 (Reference Brohan, Kennedy, Harris, Tett and JonesBrohan and others, 2006). The same is true of Nepal, with the largest temperature increases found at higher altitudes (Reference Thomas and RaiThomas and Rai, 2005). Warming leads to the shrinkage of glacial systems, so it is not surprising that Himalayan glaciers have continuously receded over recent decades.

Glacier AX010, Shorong Himal, (27˚42′ N, 86˚34′ E; Fig. 1) is one of more than 3000 glaciers in the Nepal Himalaya. Sufficient data on this glacier exist in terms of areal extent, mass balance and ice flow (e.g. Reference Fujita, Kadota, Rana, Kayastha and AgetaFujita and others, 2001) to attempt a model-based study to investigate its response to climate change. Since the first detailed studies conducted in 1978/79 (e.g. Reference Ageta, Ohata, Tanaka, Ikegami and HiguchiAgeta and others, 1980; Reference Ikegami and AgetaIkegami and Ageta, 1991; Reference Kadota, Fujita, Seko, Kayastha and AgetaKadota and others, 1997), several state variables of the glacier have been recorded in 1989, 1991, 1995–99 (e.g. Reference Fujita, Kadota, Rana, Kayastha and AgetaFujita and others, 2001) and 2004 (Reference Shrestha and ShresthaShrestha and Shrestha, 2004). In 1978, intensive observations were conducted to determine its mass balance (Reference Ageta, Ohata, Tanaka, Ikegami and HiguchiAgeta and others, 1980), heat balance (Reference Ohata and HiguchiOhata and Higuchi, 1980), areal extent (WGMS, 1998) and surface velocity (Reference Ikegami and AgetaIkegami and Ageta, 1991). In 1989, only the terminus position was measured (Reference Fujita, Kadota, Rana, Kayastha and AgetaFujita and others, 2001). In 1991, a topographic mapping of almost the whole area of the glacier was carried out (WGMS, 1998). Annual monitoring of the glacier was initiated in 1995 and continued until 1999 by means of a ground survey, the stake method, in order to obtain the mass balance, surface velocity and areal-extent data (WGMS, 2005). The ice thickness was first recorded in June 1995 (Reference Kadota, Fujita, Seko, Kayastha and AgetaKadota and others, 1997). Radio-echo sounding was conducted at three different locations on the glacier. Its areal extent was also monitored in 2004 (Reference Shrestha and ShresthaShrestha and Shrestha, 2004).

Fig. 1. Location of Shorong Himal and topographical map of glacier AX010 as surveyed in 1979 (Reference Ikegami and AgetaIkegami and Ageta, 1991, modified). Solid circles on the glacier are positions of mass-balance stakes. The dotted line along the glacier is the assumed flowline. The solid lines across the glacier are 20m interval contour lines. The grey spot labelled ‘pond’ next to the glacier terminus is a glacier-fed lake.

Besides these observational records, there have been a few model-based studies of glacier AX010. The mass-balance model by Reference Ageta and KadotaAgeta and Kadota (1992) showed that its mass balance is highly sensitive to summer temperature, not only because of summer melting but also because the glacier is in the monsoon climate regime. It therefore receives almost all of its precipitation in summer, with temperature determining the fraction accumulating as snow. Furthermore, Reference Kayastha, Ohata and AgetaKayastha and others (1999), based on their energy-balance model, demonstrated the high sensitivity of the mass balance to global radiation and the surface albedo of the glacier. Reference Kadota, Fujita, Seko, Kayastha and AgetaKadota and others (1997) projected the surface profile for 2005 and 2015 for unchanged climatic conditions from 1995 onwards (summer temperature of 2.9˚C at the snout, and summer precipitation of 1.4 m). Reference Shrestha and ShresthaShrestha and Shrestha (2004) reported that the glacier surface measured in summer 2004 was remarkably close to Reference Kadota, Fujita, Seko, Kayastha and AgetaKadota and others’ (1997) prediction for 2005.

In this study, the relation between glacier AX010 and local climate change is investigated under a series of reasonable assumptions. This provides an opportunity to understand how a small valley glacier in the Nepal Himalaya responds to climate change. In addition, we investigate which climatic parameter played the decisive role in driving the history of the glacier terminus fluctuations during the past two centuries.

2. The Numerical Ice-Flow Model

The numerical ice-flow model used in this study describes ice flow along a central flowline. The prognostic equation of this glacier model is the continuity equation, describing how the change in ice thickness is related to the flux divergence and the net surface mass balance. Along a flowline, conservation of mass can be approximated as (e.g. Reference OerlemansOerlemans, 2001):

where H is ice thickness, w is glacier width, M is mean annual net surface balance, t is time, x is distance along the flowline and U is vertically averaged velocity parallel to the bed. This velocity is comprised of two components, one originating from deformation U _d and the other from basal sliding U _s:

Reference OerlemansOerlemans (2001) reported that the flow parameter values f _d = 1.9×10^–24 Pa^–3 s^–1 and f _s = 5.7×10^–20 Pa^–3m² s^–1 provided good results when applied to several glaciers. However, it is not obvious a priori that these are also the best values for glacier AX010. These parameters depend on the specifics of a glacier, primarily its geometric setting, and may also vary temporally at a particular glacier, thereby serving a tuning purpose. We therefore adopt the above-mentioned values only as a reasonable starting point. A detailed description of how the flow parameters were tuned for our glacier is given in section 3. Equation (2) results from a generalized flow law and relates the total volume flux through a vertical section to the driving stress S _d:

where ρ is ice density (850 kgm ^–3), g is gravitational acceleration (9.81ms^–2) and h is surface elevation.

Although the model is basically one-dimensional, the two-dimensional glacier geometry is accounted for by a varying width factor for each gridpoint (see Reference OerlemansOerlemans, 1997). Assuming a trapezoidal-shaped valley, the width factor (w in Equation (1)) depends on the local ice thickness at the central flowline and the side slope of the valley cross-section:

Here w ₀ is the bottom width of the glacier. The second term on the right-hand side accounts for the width contribution of the valley’s side walls inclined at a slope of 45˚.

A finite-difference scheme was used to numerically integrate the continuity equation (Equation (1)) for ice thickness forward in time, under the additional condition that ice thickness cannot be negative (i.e. H _i≥0). The initial condition at t = 0 year is zero ice thickness everywhere (i.e. Hi(t=0) = 0m). The time integration is done with a forward explicit scheme. A time-step of 0.01 years was found to be sufficient to obtain a stable and accurate solution. This was checked by considering a steady-state glacier with a simple geometry (i.e. a constant bed slope and a unit glacier width with full lateral symmetry together with a linear mass-balance profile with elevation). Considering a linear stability criterion, the model time-step (Δt) should satisfy the condition that Δt≤(Δx)²/4D. Here Δx is the horizontal grid spacing and D is the average diffusivity (of the order of 2000 m² a^–1 on our glacier). A horizontal grid spacing of 10m was employed throughout this study, which satisfies the above criterion. Further steady-state experiments were performed to confirm whether the model conserves mass, and whether the ice flux, ice thickness and ice velocity along the central flowline were in agreement with theoretical considerations. As expected for our model set-up, all results exhibited appropriate symmetry with respect to the midway point. Furthermore, the ice flux anywhere in the glacier was equal to the integrated upstream mass balance, and the mean surface mass balance over the entire glacier was zero to within an acceptable accuracy.

3. Basic Input Data for Glacier AX010

In this section, we discuss how basic input data to run the model on glacier AX010 were collected. The input data related to the glacier geometry were extracted from a topographic map corresponding to the year 1996. Due to the fairly simple geometry of the glacier, a central flowline was easily constructed, ensuring that it remained perpendicular to all the contours. This line runs from the head of the glacier for 3 km downstream. A digital elevation model (DEM) was prepared based on the topographic map. From the DEM, we extracted the surface profile of the glacier along the assumed flowline. In addition, we accounted for the ice thicknesses measured in 1995 (Reference Kadota, Fujita, Seko, Kayastha and AgetaKadota and others, 1997). The expression that exhibits the relation between the ice thickness and the local driving stress (Equation (3)) was used to reconstruct the glacier bed profile in places where direct ice-thickness measurements were lacking (Reference OerlemansOerlemans, 2001). We also recorded the surface width of the glacier from the DEM by simply measuring the length of a contour (within the glacier domain) that passes through every gridpoint. The bottom width of the glacier at each gridpoint was then calculated using Equation (4).

Besides the extraction of the geometric input, we plotted the mass-balance data recorded during the years 1995–99 against the difference between their elevation and the corresponding equilibrium-line altitude (ELA). Before plotting, the 1998 dataset was smoothed by taking a two-point average. A few outliers were removed to ensure a good fit around the ELA. The measurements were well fitted linearly, giving a value for the balance gradient of 0.1 mw.e. a^–1m^–1, with a correlation coefficient of 0.91 (Fig. 2). This value for the balance gradient is corroborated by the value of 10mmm^–1 found by Reference Harper and HumphreyHarper and Humphrey (2003). At any time t, the mass-balance profile at the glacier surface, M(h), can now be described as:

Fig. 2. Extraction of the balance gradient: mass-balance data recorded during 1995–99 (WGMS, 2005) are plotted against the elevation with respect to the ELA. This gives a balance gradient of 0.01 mw.e. a^–1m^–1 (R ² = 0.91). Triangles represent data outliers that were not used for the analysis.

Here the first term on the right-hand side is the reference balance profile, independent of time, and the second term represents the balance perturbation, assumed independent of altitude h.

With these geometric data and the associated mass-balance profile, the model was initially run with the flow parameters reported in section 2. ΔM was set to 0 mw.e. a^–1. By slightly adjusting the ELA, a steady-state glacier was produced, whose terminus position was at the same location as on the 1996 topographic map. The corresponding value of the ELA was 5174ma.s.l. Comparing the modelled with the observed ice surface gave a poor fit for several parts of the glacier. Although some of these mismatches may merely reflect the non-steady-state conditions of the real glacier, we chose to fine-tune the model by adjusting the flow parameters until a good fit was obtained. We could have obtained a similarly good match between observed and simulated surface elevation by adjusting the bed topography (Reference OerlemansOerlemans, 1997) but opted to respect the known ice-thickness measurements. A best fit between the modelled and observed ice surface was obtained by dividing both flow parameters by a factor of 30. We conserved the ratio between the sliding and deformation parameters f _s and f _d, mainly because we had no field evidence to support a different partitioning between these two flow mechanisms.

Finally a few small adjustments were made to the bed profile in zones without thickness measurements. The final adjusted bed is shown in Figure 3. Here the ice surface results from a time-dependent simulation with historic mass-balance forcing as discussed further below (section 5). Down from x = 1300 m, the bottom width was taken as a constant value of 80 m. These adjusted parameters were used for all the subsequent experiments.

Fig. 3. Longitudinal profile of the glacier along the central flowline for AD 1996. The black bold curve is the bed elevation reconstructed from the observed surface elevation (thin curve) using Equation (3). The dotted curve is the simulated surface elevation. Open circles along the bed profile represent the locations where ice thickness was measured.

4. Basic Sensitivity Experiments

In the model, the climate forcing is introduced by variations in the local net mass balance. Starting from the net balance curve, mass-balance variations can be imposed by assuming either a balance shift or an ELA shift. Since the mass balance versus elevation with respect to the ELA exhibits a linear relationship (Fig. 2), the two approaches are equivalent. Here we express mass-balance changes as a balance shift as was done by Reference Huybrechts, de Nooze, Decleir and OerlemansHuybrechts and others (1989).

Basic model sensitivity versus stepwise perturbations in the net mass balance is shown in Figure 4a. It is clear that the model exhibits a similar behaviour for both positive and negative shifts in the mass balance except for a higher perturbation in negative mass balance (<–0.20mw.e. a^–1). An abrupt increase in sensitivity can be seen for a mass-balance perturbation between –0.20 and –0.25ma^–1. The total length of the glacier suddenly reduces from 1050 to 650 m for the respective perturbations. Closer inspection of the ice-thickness profile reveals a geometric origin. There exists a thin ice layer downstream from x∽900m for a balance perturbation of –0.20ma^–1 that is associated with a bedrock sill (Fig. 3). Increasing the negative mass-balance perturbation from –0.20ma^–1 to –0.25ma^–1 is apparently enough to melt the entire tongue over a distance of almost 500 m.

Fig. 4. (a) Steady-state glacier length versus mass-balance perturbations; and (b) reaction of the glacier front position to a stepwise change in surface mass balance of given magnitude (mw.e. a^–1). The corresponding e-folding length response time (τ) is given in years.

Next we explore the length response time of glacier AX010. This timescale is defined as the e-folding time period required to reach (1 – e^–1), i.e. 63% of the final response after a stepwise perturbation in climate conditions. This intrinsic timescale yields information on how quickly a glacier adjusts to a sudden change in mass balance. This knowledge is very useful to determine how far back in time the model needs to be spun up (e.g. to simulate the observed historical front variations, as discussed in section 5).

To estimate the length response time, stepwise changes in the mass balance in the range –0.5 to +0.5mw.e. a^–1 were applied to a steady-state reference glacier corresponding to the year 1996. Changes in glacier length and their corresponding response times are shown in Figure 4b. The presence of the icefall between x∽900 and x∽1200m along the longitudinal glacier profile causes a rapid retreat of the glacier starting from ∽1200m for the larger negative mass-balance changes. This can be clearly seen in the figure for a mass-balance perturbation of –0.50ma^–1. The length response time of glacier AX010 was found to be in the range 50–85 years and similar for both positive and negative perturbations. Although this estimate matches well with the length response time of other glaciers as estimated using the same method for Glacier d’Argentière, France, (27–45 years; Reference Huybrechts, de Nooze, Decleir and OerlemansHuybrechts and others, 1989) Nigardsbreen, Norway, (63– 73 years; Reference OerlemansOerlemans, 1997) and Sofiyskiy glacier, Russian Altai, (73–84 years; Reference De Smedt and PattynDe Smedt and Pattyn, 2003), we believe that our result is quite high for such a small subtropical valley glacier. Such a high value results mainly from the adjustments made in the flow parameters (reduction by a factor of 30): the smaller these parameters, the stiffer the ice becomes, thereby making the glacier respond more slowly to a climate forcing. In fact, reducing the flow parameters is equivalent to making the ice colder. This might, however, be true for high-altitude Himalayan valley glaciers such as glacier AX010.

5. Simulation of Historical Front Variations

The fairly short record of the terminus positions clearly reflects a continuous receding trend of glacier AX010 since 1978 (Reference Fujita, Kadota, Rana, Kayastha and AgetaFujita and other, 2001; Reference Shrestha and ShresthaShrestha and Shrestha, 2004). In addition, it is believed that the glacial systems in the Nepal Himalaya were much longer during the Little Ice Age (LIA) than today (Reference Thomas and RaiThomas and Rai, 2005). With this information, an attempt was made to simulate the historic record of the glacier’s front positions. A comparison between modelled and observed front positions is expected to explain more about the type of climatic variables decisively affecting the glacier’s mass balance (Reference Huybrechts, de Nooze, Decleir and OerlemansHuybrechts and others, 1989). Since the mass-balance record of the glacier is limited to a few years, a proper forcing function has to be established. Here we implemented mass-balance perturbations as a function of temperature and precipitation separately. As glacier AX010 is a summer-accumulation-type glacier, summer balance is considered representative for the annual balance. The winter precipitation, which accounts for <20% of the annual value, does not contribute much to accumulation, as it is blown away by high-speed gusts that can exceed 140 kmh^–1, especially in January (Reference Tartari, Verza and BertolamiTartari and others, 1998), and ablation in winter is negligible.

A first experiment was conducted under the assumption that the mass-balance perturbation is a sole function of summer temperature T _s. We combined T _s data from Darjeeling and Kathmandu stations (Fig. 5a), which were selected as the nearest stations to the glacier (∽150km away, to the east and west respectively) that have a long series of data. As the corresponding climatic anomalies of these stations match well for overlapping periods, one can assume that these datasets are likely to represent climatic variations around the glacier well. These data only extend back to 1882, so proxy data have to be used for the preceding period. We have used the general trend of the temperature anomaly derived from the Dasuopu (Qinghai– Tibetan Plateau) ice core for the period 1600–1881 (Reference Feng and HuFeng and Hu, 2005).

Fig. 5. (a)Darjeeling and Kathmandu summer mean temperature anomalies with respect to the 1951–80 mean. (b) Mass-balance perturbation solely based upon the temperature anomaly, constructed by combining a general trend of T _s anomaly derived from the Dasuopu ice core (Reference Feng and HuFeng and Hu, 2005) for the period before 1882, the Darjeeling anomaly for 1882–1980, and the Kathmandu anomaly thereafter. (c) Simulated and observed front positions.

Assuming that the same linear balance model also applies to the past, the mass balance is perturbed, ΔM (mw.e. a^–1), equally for all altitudes in the experiments as follows:

Here ΔT (˚C) is the temperature anomaly, C ₁ (mw.e. a^–1 ˚C^–1) controls the terminus response amplitude, and C ₂ (˚C) is a parameter to be optimized in order to obtain the right glacier length (Reference Huybrechts, de Nooze, Decleir and OerlemansHuybrechts and others, 1989). Parameter C ₁ gives the sensitivity of the mass balance to the summer temperature. This glacier property depends on the geometric and climatic setting, and is unique for each glacier. Parameter C ₂, on the other hand, serves a tuning purpose. It controls the absolute length of the glacier but hardly affects the amplitude of length variations. Both parameters are considered invariant in time.

The value of C ₁ for glacier AX010 has been derived from the mass-balance data recorded during the years 1995–99 assuming an atmospheric lapse rate of –6.5˚Ckm^–1. Excluding the data of 1998 when the ELA was higher than the highest altitude of the glacier, the maximum difference in ELA is 23 m. The corresponding difference in mass balance is 0.111 m (WGMS, 2005). Consequently, a 1˚C change in temperature causes the mass balance to change by 0.7 mw.e. a^–1. Our value of C ₁ = 0.7 mw.e. a^–1 ˚C^–1 agrees well with the range of mass-balance sensitivities to summer temperature (0.7–1.0mw.e. a^–1 ˚C^–1) as estimated by Reference Kayastha, Ohata and AgetaKayastha and others (1999).

Our simulation of historic front variations started at AD 1200 (t = 0 year) with zero ice thickness. For the mass balance, the ELA was set at 5210ma.s.l. (the 1995–99 mean) with a value of ΔM = +0.32mw.e. a^–1 until AD 1600. This period to spin the model up is long enough to have negligible effect on the reconstructed history of the terminus position since the LIA. Normally, a glacier carries in its memory the mass-balance history over a period equivalent to one or two times the response time (Reference OerlemansOerlemans, 1997). After AD 1600, the model was forced with ΔM values (Equation (6)) derived from climate reconstructions. The best fit between the modelled and observed front positions (Fig. 5c) was achieved for C ₂ = +0.05˚C. The corresponding evolution of mass-balance forcing is shown in Figure 5b.

As shown in Figure 5c, the model projects that the glacier was slightly advancing up to a maximum length of around 2500 m until ∽AD 1850, with a general retreat thereafter. The simulated and the observed history of retreat match remarkably well (root-mean-square error (rmse) ∽27 m; linear correlation coefficient 0.95). This seems to suggest that the glacier length history can be defined well by the temperature anomaly alone. In Figure 3, a comparison between the modelled and observed AD 1996 glacier stand can be seen in terms of ice thickness (rmse ∽11 m). The model underestimates the ice thickness in the accumulation zone, but this disagreement will probably have little effect on the results of the historical front simulations, as variations in front position are mostly governed by changes in ablation area (Reference De Smedt and PattynDe Smedt and Pattyn, 2003).

Another, similar experiment was conducted under the assumption that the balance perturbation is a sole function of summer precipitation P _s. For this purpose, we used the precipitation anomalies from Kathmandu and Darjeeling stations (Fig. 6a), as well as proxy data derived from the Dasuopu ice core that have a correlation coefficient of 0.57 with Nepal rainfall for the period 1901–95 (Reference Duan, Yao and ThompsonDuan and others, 2004). Here also a linear mass-balance model is assumed to be applicable to the past:

Fig. 6. (a)Darjeeling and Kathmandu summer precipitation anomalies (by ratio) with respect to the 1951–80 mean. (b) Mass-balance perturbation solely based upon the precipitation anomaly, constructed by combining a general trend of P _s anomaly derived from the Dasuopu ice core (Reference Duan, Yao and ThompsonDuan and others, 2004) for the period before 1870, the Darjeeling anomaly for 1870–1972 and the Kathmandu anomaly thereafter. (c) Simulated and observed terminus positions.

where ΔP is the dimensionless precipitation anomaly expressed as a ratio. The value of C ₃ (mw.e. a^–1) describes the amplitude of the terminus response. Likewise, C ₄ (mw.e. a^–1) is a tuning parameter controlling mean glacier length.

Parameter C ₃ was estimated from the empirical equations defined by Reference Kadota and AgetaKadota and Ageta (1992) based on records of mass balance and climatic parameters in 1978/79. During that season, T _s and P _s at the glacier terminus (4958ma.s.l.) were recorded as 2.4˚C and 1.45 m, respectively. Using the same atmospheric lapse rate of –6.5˚Ckm^–1, we found that the mean T _s at the glacier head (5360ma.s.l.) would be about – 0.2˚C. These equations were used to plot the mass balance against temperature on glacier AX010 for a varying precipitation. The plot revealed that a 10% change (both positive and negative) in precipitation, i.e. 0.145 m, causes a variation in mass balance between 0.145 m (for T _s≤–3.2˚C) and 0 m (for T _s≥3.4˚C). For the same change in P _s, the average value of mass-balance change for a temperature range –0.2˚C (at the glacier head: 0.123 mw.e. a^–1) to +2.4˚C (at the terminus: 0.036 mw.e. a^–1) was found to be 0.08 mw.e. a^–1. This means that, assuming the linear balance model, a 1 m change in precipitation results in a 0.55 mw.e. a^–1 change in mass balance, which is used as the value of C ₃. The physical meaning is that on average 55% of the summer precipitation falls as snow and nourishes the glacier.

Further experimentation revealed a suitable value of C ₄ = 0.40mw.e. a^–1. The resulting mass-balance perturbations and front positions are shown in Figure 6b and c. The figure reveals that the length history of glacier AX010 as solely defined by the precipitation forcing does not match the recorded length history appreciably (rmse ∽110 m). This suggests that precipitation by itself is not a good predictor to explain the glacier history. This result could have been anticipated since there is a little correspondence between the longer-term precipitation trend and glacier length. To investigate the future evolution of glacier AX010, it may thus suffice to assume that temperature can solely define the glacier’s response.

6. Future Evolution of the Glacier

Reference Shrestha and ShresthaShrestha and Shrestha (2004) mentioned that the observed trends in temperature and precipitation in Nepal are in good agreement with most climate-model predictions. We therefore opted to choose projected scenarios rather than prescribe fully arbitrary ones to force the model in the future. All the scenarios considered in this study are based upon the Intergovernmental Panel on Climate Change (IPCC) projections for the high-altitude Himalaya and the Tibetan Plateau (Reference Christensen and SolomonChristensen and others, 2007). We selected three different scenarios of temperature change by the end of the 21st century: +3˚C, +4.5˚C and +6˚C with respect to the 1980–99 mean. Likewise, five different scenarios of precipitation change were chosen in the range –30% to +30% with respect to the 1980–99 mean. The future evolution simulations follow on from the experiment of the historical front variations as solely defined by the temperature forcing.

In the first series of experiments we keep the precipitation constant at the 1980–99 mean. Only the temperature anomaly is varied, linearly by +3˚C, +4.5˚C and +6˚C by AD 2100, respectively. Smaller temperature rises are excluded because +3˚C and +6˚C are the lower and upper limits on IPCC projections for the region considered. Figure 7a displays the applied mass-balance perturbation throughout this century. The resulting evolution of the glacier length is shown in Figure 7b. As shown by the solid curve in this figure, the glacier length is reduced to 940 m by AD 2100 if the temperature remains constant at the 1980–99 mean. However, the glacier ceases to exist at AD 2083, 2056 or 2049 if the temperature increases linearly by the end of the 21st century by +3˚C, +4.5˚C or +6˚C, respectively. During its decay, the glacier retains ice mass in a number of segments depending upon the basal topography. The possible number of segments for different scenarios of the temperature rise is shown in Figure 7c. It is seen that if the temperature remains constant, the glacier will not disintegrate into segments within this century. For all of the rising-temperature scenarios considered here, however, the glacier breaks into a number of pieces before complete drawdown. Fragmentation starts at around AD 2020, with the maximum possible number of segments equal to three.

Fig. 7. (a)Applied mass-balance perturbations, dM; (b) the corresponding evolution of glacier length; and (c) the number of glacier segments for the first series of scenarios with constant precipitation (at the 1980–99 mean) and temperature varied in the range +3 to +6˚C by AD 2100.

Before considering another series of experiments with a varying precipitation rate, we explore the threshold temperature for the complete drawdown of the glacier within this century, by imposing temperature rises of 0–6˚C in 0.1 ˚C increments. We determine +2.5˚C as the threshold value above which the glacier disappears completely within this century. The year of glacier drawdown for each 0.1˚C increment in temperature (above the threshold value) is plotted in Figure 8. The general nature of the curve is exponential. If the temperature increases by less than 2.5˚C with respect to the 1980–99 mean, glacier AX010 is expected to survive until at least the end of this century. On the other hand, if the temperature rises linearly by more than 6˚C by AD 2100, the glacier will only last until the 2040s.

Fig. 8. Time (year AD) of complete glacier drawdown as a function of a linear temperature rise until AD 2100. For complete melting of glacier AX010 within this century, the temperature rise by AD 2100 must be higher than a threshold of 2.5˚C with respect to the 1980–99 mean.

In the next series of experiments, we maintain a constant rise in temperature but additionally vary the precipitation in the range –30 to +30%. For a constant rise in temperature of 3˚C by AD 2100, the model projects that the glacier will vanish during the 18 year period AD 2075–93, the timing being determined by the percentage of precipitation change: AD 2075, 2079, 2083, 2088 or 2093 if the precipitation linearly changes by –30%, –15%, 0%, +15% or +30% respectively by AD 2100. In all these scenarios, the glacier is likely to break into two pieces between ∽AD 2030 and ∽AD 2065, thereafter retaining ice only in a single segment.

If the temperature rises linearly by 4.5˚C by AD 2100, for all scenarios of precipitation change the glacier will vanish within the 6 year period AD 2054–60. The model projects fragmentation into a maximum of three pieces between AD 2026 and 2029 for all scenarios of precipitation change, and then ice retention in two segments until the year of glacier drawdown.

If the temperature rises linearly by 6˚C by AD 2100, for all scenarios of precipitation change the glacier will vanish within the 2 year period AD 2048–50. It is likely to fragment into a maximum of three segments between AD 2025 and 2027, then retain the ice mass in two patches until AD 2038–42 depending on the percentage of precipitation change. The remaining ice is then retained in a single segment until the year of glacier drawdown.

These results demonstrate that the given range of precipitation change (–30 to +30%) produces a narrower range of year of glacier drawdown (2 years) for a higher temperature rise (+6˚C) than for a smaller temperature rise (18 years at +3˚C). This implies that for a higher temperature change, precipitation will have less effect on mass balance, and thus affect the glacier’s behaviour minimally. This could be mainly due to the overriding effects of temperature over precipitation not only by enhancing the rate of ablation but also by reducing the fraction of snow accumulation.

7. Conclusions and Recommendations

The simulation of historical front variations of glacier AX010 shows that temperature, but not precipitation, can solely define the history of glacier fluctuations. This does not imply that precipitation plays no role in defining the evolution of Himalayan glaciers such as glacier AX010 that are mainly fed by the Asian summer monsoon; but it indicates the overriding effect of summer temperature over summer precipitation in defining ice melting and the snow accumulation fraction. Had we had better quantitative information regarding 19th-century glacier front positions, we could have refined the ratio of the contribution of temperature to that of precipitation in defining the historical fluctuations of the glacier. Investigation of the future behaviour of the glacier confirms that it will shrink continuously throughout this century, even if the climatic situation remains constant at the 1980–99 mean. The glacier is projected to vanish completely within this century for all trends with a temperature rise of ≥2.5˚C by AD 2100 with respect to the 1980–99 mean and precipitation changes tested in this study. The fidelity of these projections is restricted by the assumptions made in this study. An important constraint of the model is probably the choice of rather low flow parameters. A more extensive set of ice-thickness measurements on more locations of the glacier would also have helped to minimize possible errors made in defining the bed profile. In addition, the collection of any information left by end moraines would have been useful to better validate the model during the LIA glacier maximum.