Introduction

In high-altitude areas of the Hindu Kush-Himalayan (HKH) region, where meteorological observations are limited, glacier changes are valuable indicators of the effects of climate change (e.g. Reference OerlemansOerlemans, 2005). Furthermore, changing glaciers in this area influence the river discharge from one of the largest continental ice reservoirs outside the polar regions, which sustains the livelihood of several hundred million people (Reference Barnett, Adam and LettenmaierBarnett and others, 2005; Reference Dyurgerov and MeierDyurgerov and Meier, 2005; Reference Immerzeel, Van Beeke and BierkensImmerzeel and others, 2010). There has been recent debate on the speed of glacier changes in the HKH region (Reference CruzCruz and others, 2007; Reference RainaRaina, 2009; Reference Cogley, Kargel, Kaser and Van der Veen CJCogley and others, 2010), but most studies have documented widespread glacier retreat and thinning during the last century in the Himalaya and Tibet (e.g. Reference Mayewski and JeschkeMayewski and Jeschke, 1979; Reference Fujita, Nakawo, Fujii and PaudyalFujita and others, 1997; Reference Ageta, Iwata, Yabuki, Naito, Sakai and NaramaAgeta and others, 2000; Reference CogleyCogley, 2011). Some studies even suggest 21st-century acceleration of glacier thinning (Reference Berthier, Arnaud, Kumar, Ahmad, Wagnon and ChevallierBerthier and others, 2007; Reference Bolch, Buchroithner, Pieczonka and KunertBolch and others, 2008), but conclusive evidence is limited. In the Karakoram, however, several glaciers have been reported to thicken and advance (Reference HewittHewitt, 2005; Reference Scherler, Bookhagen and StreckerScherler and others, 2011a).

The reason for the different glacier behaviour in the Karakoram is not known, but there are hints that climatic changes associated with the mid-latitude westerlies could account for positive specific mass balances. For example, meteorological station data (Reference Archer and FowlerArcher and Fowler, 2004; Reference Fowler and ArcherFowler and Archer, 2006) and tree-ring records from the Karakoram (Reference TreydteTreydte and others, 2006) and the western Himalaya (Reference Yadav, Park, Singh and DubeyYadav and others, 2004) document decreasing summer temperatures and increasing precipitation during the second half of the 20th century. Concurrently, climate reanalysis data suggest strengthening and an altitudinal descent of the westerly jet stream over Central Asia (Reference Archer and CaldeiraArcher and Caldeira, 2008), which is the principal engine for transporting westerly moisture from the Mediterranean, Black and Caspian seas to the Karakoram. However, mass-balance measurements on glaciers are scarce and preclude a thorough assessment of trends in specific mass balances. In addition, the current political situation in northern Pakistan complicates extended field campaigns, and the installation and maintenance of dense monitoring networks on large and frequently debris- covered glaciers is extremely difficult.

Alternatively, remote-sensing based glacier surface velocities may be used as indicators of glacier conditions (e.g. Reference Scherler, Bookhagen and StreckerScherler and others, 2011a). For example, Reference Quincey, Copland, Mayer, Bishop, Luckman and BellaQuincey and others (2009) studied the surface velocity of Baltoro Glacier, central Karakoram, between 1993 and 2008 and linked velocity variations during this period to climatic changes. In particular, gradually increasing annual and winter velocities have been explained by a positive mass balance due to long-term increases in precipitation and thus accumulation in this region (Reference Quincey, Copland, Mayer, Bishop, Luckman and BellaQuincey and others, 2009). However, the velocity maximum of Baltoro Glacier in summer 2005 was attributed by the authors to high meltwater fluxes and increased basal sliding, following heavy winter snowfall. Such links between weather and short-term glacier sliding are well known from other temperate glaciers (e.g. Reference WillisWillis, 1995), but complicate the distinction between short-term anomalies and long-term trends. Another complicating fact is that surging glaciers are particularly common in the Karakoram (Reference HewittHewitt, 1998, Reference Hewitt2007; Reference Barrand and MurrayBarrand and Murray, 2006), which obey quasi-regular cyclic velocity variations, with a wide range of surge periods and magnitudes (Reference HewittHewitt, 1998;Reference Mayer, Fowler, Lambrecht and ScharrerMayer and others, 2011; Reference Quincey, Braun, Glasser, Bishop, Hewitt and LuckmanQuincey and others, 2011). It thus appears that differentiating velocity increases due to mass-balance effects from periods of increased sliding or surging phenomena are important, but may require long and detailed velocity records.

In this study, we use a dense time series of Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) optical satellite images to conduct a detailed investigation of the recent (2000-09) velocity history of the 63 km long Biafo Glacier, central Karakoram (Fig. 1). This dataset includes information from 2005, when relatively high velocities were reported from Baltoro Glacier (Reference Quincey, Copland, Mayer, Bishop, Luckman and BellaQuincey and others, 2009). Because Biafo Glacier is relatively fast, accurate measurements with the feature-tracking software COSI-Corr (Reference Leprince, Barbot, Ayoub and AvouacLeprince and others, 2007) enable us to analyse the pattern and evolution of velocites at high spatial and temporal resolution, with the aim of potentially discriminating between different mechanisms accounting for any velocity variations. Such an assessment helps us to put the remotely sensed glacier velocities into context with ice dynamics and climatic changes.

Fig. 1. Biafo Glacier, central Karakoram. (a) Orthorectified ASTER satellite image (band 3N, 010101), with glacier outline. Longitudinal profile following the central flowline is shown in white, with points spaced 3 km apart. Black contours on glacier surface are at 100m intervals. ‘B’ marks Baintha profile, with velocity measurements in 1985-87 (Table 1). Inset shows location of study area. Grey areas >2000 m a.s.l.; black areas depict glaciers (based on Digital Chart of the World). Weather stations are 1. Gupis; 2. Gilgit; 3. Chilas; 4. Bunji; 5. Astore; 6. Skardu. (b) Surface elevation and slope along profile in (a). (c) Horizontal surface velocity along profile in (a). Grey circles depict velocity measurements with time-spans >0.5 years. Light grey area in background shows number of measurements per profile point, with y-axis to the right. (d) Time-spans over which velocities were measured. Triangles denote acquisition dates of ASTER images used in this study.

Study Area and Previous Studies

The Karakoram, northern Pakistan, is the most heavily ice- covered mountain range in Asia and hosts several of the world’s longest valley glaciers (Reference von Wissmannvon Wissmann, 1959). Snow and glacier melt collectively contribute ˜60-70% of water discharge in the Indus River (Reference Karim and VeizerKarim and Veizer, 2002; Reference Bookhagen and BurbankBookhagen and Burbank, 2010), which supplies fresh water for drinking, irrigation and power supply to millions of people living in the Indus Plains (Reference Immerzeel, Van Beeke and BierkensImmerzeel and others, 2010). Situated in the western part of the HKH region, the Karakoram receives most snowfall from the mid-latitude westerlies (e.g. Reference WakeWake, 1989; Reference Barry and ChorleyBarry and Chorley, 2003). The high topography of the mountain range forms an orographic barrier that intercepts atmospheric moisture at high elevations, whereas the valley bottoms remain relatively arid (Reference Hewitt, Wake, Young and DavidHewitt and others, 1989; Reference HewittHewitt, 2005). For example, meteorological stations in Skardu (2209ma.s.l.) and Gilgit (1459m a.s.l.), located immediately southwest of the Karakoram (Fig. 1, inset), record long-term mean annual precipitation of <0.2 m (Pakistan Meteorological Department, http://www.pakmet.com.pk/), whereas annual snowfall of up to ˜2mw.e. has been reported from elevations of ˜5000m (e.g. Reference WakeWake, 1989; Reference Winiger, Gumpert and YamoutWiniger and others, 2005).

Biafo Glacier is situated in the central Karakoram (Fig. 1). It has an extensive accumulation area at >4650 m a.s.l. (as of 1985/86; Reference Hewitt, Wake, Young and DavidHewitt and others, 1989) that comprises two major glacier branches. At its western edge, Biafo Glacier shares an ice divide with the westward-flowing Hispar Glacier. Biafo Glacier flows southeast into a long (>40km), deeply incised valley that approximately follows the orientation of the orogen, and terminates at 3030 m a.s.l., close to the Braldu River. Reference WakeWake (1989) measured annual accumulation rates of 0.9-1.9mw.e. in six snow pits dug in the accumulation area of Biafo Glacier between 4650 and 5520 m a.s.l. in 1985 and 1986. Coeval stake measurements in the ablation area (Reference Hewitt, Wake, Young and DavidHewitt and others, 1989) indicate high melt rates of 5-7mw.e. a^-1 at ˜4000m a.s.l. In summary, like other Karakoram glaciers, Biafo Glacier is characterized by heavy winter snowfall followed by intense summer melting (Reference HewittHewitt, 2005). These conditions promote the formation of significant amounts of meltwater, which influence basal sliding (Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2008).

Reference Hewitt, Wake, Young and DavidHewitt and others (1989) reported glacier surface velocities from ground measurements along a transverse profile (near km40; labeled ‘B’ in Fig. 1a), collected during midsummer 1985, early summer 1986 and winter 1985/86 (Table 1). Reference CoplandCopland and others (2009) compared Reference Hewitt, Wake, Young and DavidHewitt and others’ (1989) measurements with their own measurements obtained from correlation of ASTER images from summer 2006 and 2007, i.e. about 20 years later, and concluded that velocities from both studies overlap within uncertainties, suggesting no or only minor changes.

Materials and Methods

We used the 15 m resolution near-infrared (NIR) near-nadir 3N band from ASTER L1A satellite images and the orthorectification and correlation tool COSI-Corr (Reference Leprince, Barbot, Ayoub and AvouacLeprince and others, 2007, Reference Leprince, Berthier, Ayoub, Delacourt and Avouac2008) to obtain glacier surface velocities. The technique has been tested and used extensively in high- altitude terrain where it was found to provide surface displacements with an accuracy of one-fifth to one-tenth of the pixel resolution (Reference Leprince, Berthier, Ayoub, Delacourt and AvouacLeprince and others, 2008; Reference Scherler, Leprince and StreckerScherler and others, 2008, Reference Scherler, Bookhagen and Strecker2011b; Reference Herman, Anderson and LeprinceHerman and others, 2010). Prior to sub-pixel correlation, the images were co-registered and orthorectified using a Shuttle Radar Topography Mission (SRTM)-based gapless digital elevation model (DEM;J. de Ferranti, 3’’ resolution digital elevation data for High Asia, 2007, http://www.viewfinderpanoramas.org/dem3.html).

Table 1. Surface velocities (m d^–1) at Baintha profile, Biafo Glacier,

from Reference Hewitt, Wake, Young and DavidHewitt and others (1989)

The available satellite images (n = 46; Table 2) were acquired on 34 different days between September 2000 and March 2009 (Fig. 1d). Not all images cover the entire area of Biafo Glacier and some are partly obscured by cloud cover. In addition, extensive snowfall or snowmelt between two image acquisitions can make their correlation difficult. In total, we selected 166 image pairs for correlation, with time- spans ranging between 0.02 years (i.e. 7 days) and 5 years (mean 1.01 years). We found that time-spans between 3 months and 3 years yielded mostly good results with relatively low uncertainties, whereas for longer time-spans the correlation often fails and shorter time-spans often result in very noisy data. Our detailed analysis of the velocity history of Biafo Glacier is mostly based on 13 and 24 image pairs covering time-spans of 1-3 years and 1-7 months for interannual and seasonal velocity variations, respectively. To study glacier velocities in the vicinity of Biafo Glacier, we used 54 different image pairs with time-spans of 0.9-3 years, because of varying spatial and temporal coverage.

Post-processing steps involved the removal of residual attitude effects in the displacement maps and excluding miscorrelations with a filter tool based on a priori knowledge of flow directions from surface features, and the distribution and lateral continuity of flow speeds (Reference Scherler, Leprince and StreckerScherler and others, 2008, Reference Scherler, Bookhagen and Strecker2011b). The residual uncertainties in the displacement measurements obtained with COSI-Corr are Gaussian-distributed, with zero mean and a standard deviation that depends on the image resolution (Reference Leprince, Barbot, Ayoub and AvouacLeprince and others, 2007). These uncertainties may be due to slight physical changes of the surface between the correlated scenes, shadowing, lack of topographic resolution, uncorrected satellite-induced artefacts, aliasing, quantization, noise of the sensors, and intrinsic correlation accuracy (Reference Leprince, Barbot, Ayoub and AvouacLeprince and others, 2007). To obtain an estimate of the residual uncertainty in the displacement measurements, we took the mean and standard deviation in the east-west and north-south displacement maps over all pixels with <10 m absolute displacement (Reference Scherler, Leprince and StreckerScherler and others, 2008, Reference Scherler, Bookhagen and Strecker2011b).

We acknowledge that better estimates could be obtained by excluding all ice-covered areas and calculating residual displacement statistics with the remaining pixels. However, this step is very labour-intensive as it requires a manually generated mask hiding all glaciers, moving shadows or other moving surface features for each displacement map individually. With a large number of different correlations, this step is unfeasible. The 10m absolute displacement threshold serves as a first-order mask for miscorrelated pixels and ice-covered areas, but unintentionally leads to inclusion of slow-moving ice or shifting shadows, and results in a nonzero mean and a higher standard deviation while potentially excluding stable ground with residual offsets exceeding 10 m a^-1. However, comparing the extents of our masks with the satellite images shows that ice-free areas are much larger than ice-covered areas and only a few, isolated pixels on stable ground were omitted. Therefore, we consider our estimates to be reasonable approximations of the displacement uncertainties. We computed the uncertainties in the annual velocities by adding the mean and standard deviation of the total displacement uncertainties and dividing by the temporal separation of the correlated images. As expected, the residual uncertainties in the annual velocity measurements scale with the length of the measurement period (Fig. 2). For time-spans of 1-7 months (i.e. ˜0.1- 0.6 years), the uncertainties are ˜10-50 m a^-1, and for time- spans of 1-3 years the uncertainties are ˜1-5 ma^-1.

Table 2. List of ASTER images used in this study

Fig. 2. Residual uncertainties of the remotely sensed glacier surface velocity measurements plotted against the measurement time-span.

Results

In the following we report all velocity data collected along a profile that follows the central flowline of the glacier. In the upper part of Biafo Glacier (0-20 km from the glacier head), frequent elliptical depressions with up to 100 m and 250 m diameters, which appear to reflect strained crevasses, provide abundant and easily recognizable surface features suitable for tracking. Despite poorer image coverage, the analysis of these features results in ˜20-30 measurements per point along the profile (with measurement periods >0.5 years) and indicates relatively constant velocities during the period 2005-08 which is covered by the image pairs in this part of the glacier (Fig. 1c). The same observations, which are not shown here, apply to the eastern branch (Sim Gang) of Biafo Glacier. We have almost no data coverage at 2033 km distance from the glacier head, because of rarely overlapping images in this part of the glacier and the lack of surface features for tracking. In much of the ablation zone, down-glacier of ˜33 km from the glacier head, we have good image coverage and obtained 50-100 measurements (with 0.5 year measurement periods) per point along the profile. Velocities in this part of the glacier show high variability, with variations of up to 150 m a^-1, or ˜100% in some places (Fig. 1c). In the following, we investigate velocity variations in the lower part (>33 km from the glacier head) of Biafo Glacier in more detail.

To depict the velocity variations more accurately, we first present velocity measurements covering periods close to 1 year. However, because we lack enough images and thus velocity data during certain periods (e.g. 2002-04), we also included some measurements made over time-spans of up to 3 years, and acknowledge that these represent average velocities during the respective time intervals. The earliest measurements with time-spans of 1-3 years cover the period 2001-03 and indicate velocities of ˜150-200 m a^-1 between 43 and 58km from the glacier head (Fig. 3a). Between 2003 and 2004, Biafo Glacier accelerated in the upper part of its ablation zone by >50 m a^-1, whereas velocities in the lowermost 12 km (km 55-67) were identical to those in the first period. In the following year (2004/05), increased velocities appear further down-glacier and Biafo Glacier attained the highest velocities during the entire observation period. Starting in 2005, velocities declined again throughout the lower part, and during 2006 reached rates comparable to the initial phase (2001-03). Between 40 and 45 km from the glacier head, minimum velocities were reached by the end of 2006, while further down-glacier velocities continued to decrease until 2008, when they were up to 50 m a^-1 lower than in 2006.

Fig. 3. Horizontal surface velocities of Biafo Glacier along its lower ˜36 km. (a) Absolute annual velocities between May 2001 and March 2009, measured from 13 image pairs with time-spans of 1–3 years. Colour coding of data points refers to measurement period. (b) Time and duration of measurement periods shown in (a). Black triangles indicate acquisition dates of the satellite images. (c) Absolute annual velocities between May 2006 and June 2007, measured from 20 image pairs with time-spans of 1–7 months. (d) Time and duration of measurement periods shown in (c).

Most of the velocities shown in Figure 3a are measured over time-spans close to a year or an integer multiple of a year, so potential seasonal velocity variations should not greatly affect the remote-sensing-based velocity history. However, we acknowledge that interannual variations in the amount of sliding, for example, can potentially influence the mean annual velocities that we report. To test if Biafo Glacier exhibits seasonal velocity variations, we also correlated images only a few months apart. Velocities measured over periods of 1-7 months are noisier than for longer periods, and higher uncertainties make detailed assessment of velocity variations more difficult (Fig. 3c). Nevertheless, significant sub-annual velocity changes can be identified. Between May/June and December 2006, velocities steadily decreased in the lower ˜33 km of Biafo Glacier, with greater changes further up-glacier from the terminus. Velocities remained low during winter 2007 and started to increase over the lower 10 km of the glacier during spring and early summer 2007. Velocity differences between early summer 2006 and winter 2007 span 150-200 m a^-1 at the upper end of this zone. We note that the sub-annually resolved velocity variations are possibly superimposed by the previously identified interannual velocity changes, and therefore may contain a ‘background’ deceleration during the observation time.

Figure 4 shows the velocity history of Biafo Glacier at different positions along its lower part, measured at different time intervals and, to take the scatter of the data into account, averaged over a longitudinal distance of 2 km. The highest annual velocities during the period 2002-08 occurred in early summer 2005 at 40-42 km from the glacier head (Fig. 4a, position A), but slightly later, in autumn 2005, further down-glacier at km 52-54 (C) and km 58-60 (D). Based on the available measurements, the apparent down- glacier migration of the acceleration front is ˜70md^-1 (˜24 km a^-1). The mean deceleration during 2006 and 2007 at position A-D is about 36 m a^-2 (Fig. 4a). Velocities at the terminus (E) remained more-or-less constant throughout the observation period.

Fig. 4. Velocity variations at different positions along Biafo Glacier. Velocity measurements are based on (a) 13 image pairs that are 1-3 years apart, obtained between 2002 and 2008, and (b) 24 image pairs that are 1-7 months apart, obtained during 2006 and 2007. (c) Location of measurements shown in (a) and (b). Measurements are averaged over 2 km long down-glacier reaches, and error bars in (a) and (b) reflect corresponding standard deviation within the reaches.

We used the estimated deceleration of 36m a^-2 to detrend the sub-annually resolved velocity history between winter/spring 2006 and spring/summer 2007 (Fig. 4b;for comparison, the originally measured velocities are shown in light grey). In contrast to the data shown in Figure 3c, we included four more image pairs to extend the observation period. Due to removal of the deceleration trend, the maximum velocity differences slightly decrease, but remain high, up to ˜150ma^-1. By extending the observation period, an increase in velocities during spring and early summer 2007 also becomes apparent in the detrended data, at greater distances from the terminus, i.e. at positions A-C (Fig. 4b). Compared to the multiple-year velocity history, down-glacier differences in the timing of peak velocities are less obvious. The velocity variations appear to occur simultaneously, with maximum speeds during May/June 2006, minimum speeds during January 2007 and decreasing amplitudes down-glacier.

Selected transverse surface velocity profiles from the period 2001-08, at 47 and 56 km from the glacier head, are shown in Figure 5. Although the spatial coverage is not complete, both interannual and seasonal velocity changes are well depicted. At time-spans of 1 year and more (Fig. 5a and b), the correlation often failed near the edges of the profile, presumably due to high strain rates. Nevertheless, the flat shape of the central ˜50% of the profiles does not appear to change significantly between years. The seasonal velocity changes contain more scatter but better depict the steep decline towards the edges (Fig. 5c and d). Like the interannual changes, the maximum seasonal velocity changes have similar magnitude over most of the profile’s widths.

Fig. 5. Transverse horizontal surface velocity profiles at (a, c) 47 km and (b, d) 56 km from the glacier head, showing (a, b) interannual velocity variations and (c, d) seasonal velocity variations. Date format is mm/yy.

Motivated by our findings at Biafo Glacier, we also investigated the velocity histories of other glaciers in the vicinity, to test if similar interannual velocity variations could be detected. Of the fourteen glaciers for which we obtained usable velocity data, seven show marked interannual velocity variations during the study period (Fig. 6). Three more show some velocity variations, but of much smaller magnitude, and another one, whose velocity history is influenced by two surging tributary glaciers (Reference HewittHewitt, 2007), shows velocity variations, upstream of where the surging tributaries join the glacier, that appear similar in timing and magnitude to those of the other glaciers. The other three glaciers do not show any significant interannual velocity variations. Including Biafo Glacier, three of the studied glaciers slowed down between 2005 and 2008 (Fig. 6, A-C). Unfortunately, the data coverage does not allow reconstruction of the velocities throughout the study period for all glaciers, so we cannot determine whether this deceleration followed a period of acceleration as observed for Biafo Glacier. The other four glaciers shown in Figure 6 (D-G) accelerated, sometimes throughout the study period (2001-08). The three glaciers with minor variations and the glacier influenced by surging tributaries also accelerated.

Fig. 6. Regional overview of glacier surface velocity variations between 2001 and 2008 in the central Karakoram. Map in upper left corner gives location of profiles shown below and referenced by capital letters. Colour bar denotes the central date of the 54 measurement periods with lengths of 0.9–3 years. Note that the data coverage is not the same for all glaciers, but that the general velocity trends are visible: the velocities of glaciers A–C decrease from 2005 to 2008, whereas those of glaciers D–G increase

Discussion

Seasonal and interannual velocity variations

Our remotely sensed glacier velocity data highlight two different kinds of velocity variations, Δu. The first are monthly to seasonal (Δu ≤ 150ma^-1), with peak velocities during late spring and summer and a minimum during winter. Similar velocity variations were observed by Reference Hewitt, Wake, Young and DavidHewitt and others (1989) in 1985/86. Such seasonal velocity variations are common for temperate glaciers, where they have been linked to the development of the subglacial hydrological system (e.g. Reference Iken and BindschadlerIken and Bindschadler, 1986; Reference WillisWillis, 1995; Reference AndersonAnderson and others, 2004). According to these models, high melt rates during spring increase the flux of meltwater that penetrates to the glacier bed, increasing hydrostatic pressures and therefore basal sliding rates. This appears to be a reasonable explanation for the observed seasonal variations, given their timing and that they appear only in the ablation area of Biafo Glacier. In addition, the shape of the transverse velocity profiles, with their almost constant velocities in the centre and steep velocity gradients near the edges (Fig. 5), is indicative of so-called plug flow (Reference Fatland, Lingle and TrufferFatland and others, 2003), which has been interpreted as resulting from mostly basal sliding (Reference CoplandCopland and others, 2009;Reference Quincey, Copland, Mayer, Bishop, Luckman and BellaQuincey and others, 2009).

The second kind of velocity variation has slightly smaller amplitude (Δu ≤ 100 m a^-1) and occurred over several years, with peak velocities during 2005. Our measurements, mainly covering full years, indicate that these velocity changes are not simply related to seasonal variations. However, year-to- year differences in the total amount of basal sliding, following variations in annual meltwater fluxes, for example, are possible. Unfortunately, we do not have sub-annually resolved velocity measurements for the entire period and therefore cannot recover the full history of summer and winter velocities. Because the observed acceleration and deceleration are pronounced, occurred over relatively short periods and are absent in the accumulation area of Biafo Glacier, we deem it most unlikely that they are linked to mass-balance-related ice thickening and thinning, also considering the reported ablation rates of 5-7 m a^-1 at 4000ma.s.l. (Reference Hewitt, Wake, Young and DavidHewitt and others, 1989). Centre-line velocities measured by Reference Hewitt, Wake, Young and DavidHewitt and others (1989) ˜40km from the glacier head range from 220 m a^-1 during summer to 110m a^-1 during winter, with mean annual velocities of 165 m a^-1 (Table 2). This is close to the velocities observed at the beginning and end of our study period at this position, which suggests that the higher speeds of Biafo Glacier around 2005 mark an anomalous behaviour.

Although we cannot distinguish between ice movement due to deformation or sliding, based on our measurements of surface velocities, the spatial and temporal patterns of the speed-up episode around 2005 provide some clues to the processes that could potentially account for such large annual velocity variations. During 2005, and over a period of ˜6 months, we observe a down-glacier migration of the acceleration front and a decline of velocities after a short period of peak velocities (Fig. 4a). The subsequent deceleration was more rapid and stopped earlier at higher parts of the glacier than at lower parts, thus resembling the passing of a wave of high velocities, with some diffusion down-glacier. This pattern indicates that the forces behind the acceleration originated up-glacier and migrated down-glacier. Such a high-velocity wave could be linked to wave-like zones identified from ground observations made by Reference HewittHewitt (2005) at Biafo Glacier between 1985 and 2001. He described wave-like zones of higher ice moving down-glacier which originated from the upper ablation zone, where he observed 10-25 m of ice thickening. Could such waves of thicker ice potentially account for the observed velocity changes of up to 50% during both the acceleration and the deceleration phase (Fig. 3a)?

The glacier surface velocity, u _s, is the sum of ice flow due to shear deformation, u _d, and basal sliding, u _b. Standard ice- flow laws (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 310) relate u _d at the glacier surface to the basal shear stress, τ_b, according to

(1)

where A(T) is the temperature-dependent flow-rate factor, n is a constant usually taken to be 3 and h is the ice thickness. A typical relation for basal sliding has the form

(2)

where N is the effective basal water pressure, k is a sliding coefficient, and p and q are constants typically taken to be 3 and 1 (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 240). For wide and shallow sloping glaciers, τ_b can be approximated by

(3)

where ρ is the ice density, g is gravitational acceleration and α is the slope of the ice surface. Combining Eqns (1) and (2) with Eqn (3), differentiating with respect to the ice thickness and rearranging gives

(4)

(5)

At 41 km from the glacier head, Reference Hewitt, Wake, Young and DavidHewitt and others (1989) measured ice thicknesses of ˜700-500m. There, an ice thickness increase of ˜25 m will result in around 15-20% higher u _d and 10-15% higher u _b. Assuming that u _d and u _b are of similar magnitude (Fig. 3b) according to Eqns (3) and (4), the surface velocity will thus be 25-35% higher, less than the observed 50%. Under the same assumptions, however, an ice thickness increase of ˜40 m will result in 40-56% higher surface velocity, similar to what we observe. Such an amount of ice thickening is ˜60% greater than what Reference HewittHewitt (2005) reported between 1985 and 2001. We conclude that the observed interannual velocity variations could indeed be related to waves of thicker ice moving down-glacier, but this would require at least 40 m wave height and also variations in basal sliding. Such an elevation increase may even be detectable from remote-sensing-derived DEMs in the future.

The observed interannual variations occur only in the lower part of Biafo Glacier, similar to the seasonal velocity variations. Because we have a large data gap 20-33 km from the glacier head, we cannot depict the exact transition from variable to stable velocities. According to Reference Hewitt, Wake, Young and DavidHewitt and others (1989), the equilibrium-line altitude (ELA) of Biafo Glacier was located at 4650 m a.s.l. during 1985 and 1986, which is at the upper end of our data gap. If there was no ELA lowering in excess of ˜50-100 m - and there is no evidence for this - the velocity variations would only occur in the ablation zone of Biafo Glacier, which suggests that they are indeed related to processes involving meltwater (e.g. basal sliding).

The apparent wave speed of ˜24 km a^-1 seems to be too high for a kinematic wave, which theoretically should travel with a speed five times the surface velocity, or up to eight times the surface velocity in numerical experiments (Reference Van de Wal RSWVan de Wal and Oerlemans, 1995), i.e. with a maximum speed of ˜2 km a^-1 at Biafo Glacier. Reference HewittHewitt (2005) points out that the migrating waves he observed appeared in the last few decades, which he associates with recent climatic changes. An important question is whether these waves are regularly reoccurring features, much like glacier surges. At present, this is uncertain due to lack of long-term observations, but we suggest that the velocities near the beginning and end of our study period are probably closer to the long-term average velocity than the high velocities during 2005. Longer observation periods are needed to better constrain how regular and frequent such speed-up episodes are and what other durations and magnitudes may exist.

Meteorological conditions

The Pakistan Meteorological Department (PMD) operates six permanent weather stations near the Karakoram, at distances of 50-200 km from Biafo Glacier (Table 3). These stations are all located below 2500 m a.s.l. and therefore do not record conditions at the high elevations where the glaciers are located. However, we assume that trends in these data can be used as indicators for trends at higher elevations, because winter precipitation, the major source for accumulation (Reference WakeWake, 1989), is (1) correlated between different weather stations in the upper Indus River basin (Reference Archer and FowlerArcher and Fowler, 2004), suggesting regional signals to be reflected, and (2) correlated with summer runoff at several valley stations in the south adjoining the Jhelum River basin, which is mostly fed from snow cover at elevations higher than the weather stations (Reference Archer and CaldeiraArcher and Fowler, 2008). Summer precipitation is more convective and local, and hence not well correlated between different weather stations (Reference Archer and FowlerArcher and Fowler, 2004), so it is not clear if the weather station data are representative of summer precipitation at higher elevations. Because lateral variations in temperature between the weather stations are low (Fig. 7), we assume that the temperature records, although located at lower elevations, are representative of temporal variations over a larger region.

Fig. 7. Meteorological records from northern Pakistan spanning the years 2000-08. (a) Mean monthly precipitation from six meteorological stations. For station locations, see inset in Figure 1a. (b-d) Mean monthly anomalies of precipitation (b), maximum temperature (c) and minimum temperature (d), for the meteorological stations in Gilgit, Astore and Skardu. Black lines represent 1 year running means. Grey shaded boxes indicate anomalous periods present in almost all records. See Table 3 for a list of the stations.

Table 3. Locations of the PMD meteorological stations used in this study

Monthly precipitation between 2000 and 2008 at the six PMD stations was usually highest during winter and spring (Fig. 7a). There is also significant summer precipitation at some of these stations, which is related to the Indian monsoon. For three of the stations (Astore, Gilgit and Skardu), long-term (>30 years) average values are available, allowing us to focus on anomalies (Fig. 7b-d). Between roughly 2001 and 2005 there was an alternation of relatively cold summers and warm winters. Whereas cold summers reduce melting, warm winters can be expected to promote higher snowfall. Although mean monthly precipitation was often higher than the average value during winter and spring between 2003 and 2005 at Gilgit and Astore stations, this cannot be confirmed for Skardu station, which is closest to Biafo Glacier. After 2005, positive precipitation anomalies appear to have been less pronounced and less frequent and temperatures were higher in both winter and summer.

In summary, compared to the long-term average, meteorological conditions during the years 2001-05 may have resulted in higher amounts of snowfall, as indicated by local reports and reanalysis data (Reference Quincey, Copland, Mayer, Bishop, Luckman and BellaQuincey and others, 2009) and partly supported by the weather station data presented here. These can generate large amounts of meltwater to assist sliding processes. One could also expect that relatively colder summers may lead to a slower than normal increase in melt rates, and thus slowly increasing subglacial water discharge. This process may have sustained high basal water pressures over a relatively longer period, inducing greater amounts of sliding (Reference SchoofSchoof, 2010). While such meteorological conditions theoretically could have promoted long- lasting sliding events during spring, we find that this scenario does not explain well the wave-like passage of higher velocities, with regular acceleration and deceleration phases, as observed at Biafo Glacier. Nor does it explain why the acceleration continued at least until 2009 for some of the neighbouring glaciers. Given the scarcity of detailed ground observations, we find it difficult at present to evaluate the true nature of the velocity variations, and whether these are directly related to annual atmospheric conditions (e.g. Reference Span and KuhnSpan and Kuhn, 2003), or mostly governed by glacier dynamics, in which weather and climate play secondary roles.

Conclusions

We have presented remote-sensing-based glacier surface velocity data from Biafo Glacier, which provide evidence for high seasonal as well as interannual velocity variations occurring over a distance of 33 km up-glacier from its terminus. The interannual velocity variations occurred over a period of at least 4 years with an acceleration phase of at least 1-2 years, a short period of peak velocities, and a deceleration phase of at least 1.5 years, depending on the longitudinal position along the glacier. The wave-like down- glacier migration of the speed-up episode suggests linkages with waves of thicker ice reported previously from Biafo Glacier (Reference HewittHewitt, 2005). The causes for these interannual velocity variations are not known, but the restriction to the ablation zone, and the magnitude and pattern suggest that they are at least partly related to variations in basal sliding, although a combination with thickness and thus ice-flow changes is possible.

Apparently similar, but mostly longer, multiple-year acceleration phases are observed at other glaciers in the vicinity of Biafo Glacier, suggesting that this speed-up episode is a regional rather than a local feature. In the early 2000s up to 2005, several glaciers in this region attained very high velocities and several glaciers started to surge. Analysis of mean monthly meteorological data from nearby weather stations indicates that during the early 2000s, and compared to the long-term average, winters have been warmer, summers cooler and winter precipitation may have been higher. While these meteorological conditions could promote long-lasting sliding events during spring, they do not explain the wave-like and regular nature of the speed-up episode and why the acceleration continued until 2009 for some neighbouring glaciers. The causes and range of such significant interannual velocity variations, as presented in this paper, need to be better understood before inferences about the long-term effects of mass-balance changes can be made in data-scarce regions like the Karakoram.