Introduction
Surging glaciers occur in most mountain ranges of the world, but in the Karakoram a concentration of periodically fast-moving glaciers can be observed (Reference HewittHewitt, 2007). About 37% of this extensive mountain range between the Hindu Kush in the west and the Himalaya in the south and east is covered with glaciers (Reference Shroder, Bishop, WIlliams and FerrignoShroder and Bishop, 2010) and provides many of the largest glaciers outside the polar regions. In the centre, major glaciers, like Baltoro Glacier (∼64 km) and Siachen Glacier (∼72 km), drain from a small mountain region consisting of some of the highest mountains in the world, with four peaks higher than 8000 m a.s.l. Compared with the high elevations in the Himalaya, this region is dry, with almost desert-like valley floors. However, a strong precipitation gradient exists between the low parts of the valleys and the high mountain slopes, where the maximum precipitation is expected to occur between 5000 and 6000 m a.s.l. (Batura Glacier Investigation Group, 1979; Reference WakeWake, 1989; Reference Miehe, Winiger, Böhner and ZhangMiehe and others, 2001) and reaches 2000 mm a−1 . Most of the precipitation occurs during winter and spring, when westerly winds dominate the large-scale weather system. However, variations are large, and, depending on the strength of the monsoon, considerable precipitation amounts can also be accumulated at high elevations during summer (Reference WakeWake, 1989). Since the 1960s a positive change has been observed for annual precipitation, while winter precipitation seems to be declining, at least at some meteorological stations in the region (Reference Archer and FowlerArcher and Fowler, 2004). The topography is characterized by extensive valley systems and large elevation gradients in the headwaters. Often the accumulation basins, which generally occupy the elevation band of maximum precipitation, are still towered over by mountain slopes 2000–3000 m high. Consequently, many of the glaciers are fed, to a considerable extent, by avalanche snow (Reference HewittHewitt, 1998; Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006).
From 1978 onwards (with the construction of the Karakoram Highway) glacier research activities increased in the Karakoram, starting with the Chinese Batura Glacier Investigation Group (1979) and other projects like the Canadian Snow and Ice Hydrology project (e.g. Reference Hewitt and YoungHewitt and Young, 1990). In more recent times these activities have also included remote-sensing studies, focusing on the glaciers in this region (e.g. Reference Bishop, Kargel, Kieffer, MacKinnon, Raup and ShroderBishop and others, 2000; Reference Osipova and TsvetkovOsipova and Tsvetkov, 2002; Reference Quincey, Copland, Mayer, Bishop, Luckman and BelóQuincey and others, 2009). Today the overall glacier mass balance is thought to be negative (Reference Shroder, Bishop, WIlliams and FerrignoShroder and Bishop, 2010), but a considerable number of advancing glaciers were observed during the 1990s (Reference HewittHewitt, 2005). At the northeastern side of the Karakoram the conditions changed from warm/dry to warm/humid with a clear trend towards more negative mass balances in the early 1990s, estimated from river discharge measurements (Reference Gao, Zhang, Ye and QiaoGao and others, 2010).
The glaciers in the Karakoram and in Alaska, USA, represent ∼90% of the known surging glaciers (Reference HewittHewitt, 1969); in the Pamirs a large number of glaciers are potential surge-type glaciers, with >60 observed surge events (Reference Kotlyakov, Osipova and TsvetkovKotlyakov and others, 1997). In the Karakoram >30 surges have been reported since the 1860s (Reference HewittHewitt, 1998), and many more probably went unnoticed. During recent times it seems that surge activity has been increasing (Reference HewittHewitt, 2007). Many of the reported surges in the Karakoram occur in tributaries of large glaciers, while the main glaciers show no signs of strong periodic advances. These tributaries are often steep, with an altitude range of >2000–3000 m within only 10–15 km (Reference HewittHewitt, 2007). Most of these surge events last only a few months. Even for normal dynamic conditions, the surging glaciers in the Karakoram can be expected to provide a large mass transport from the accumulation basins, with high precipitation, into the drier, lower areas. Glacier length seems not to be an indicator for surge sensitivity, as the longest glaciers in the Karakoram are not known to surge, while most surging glaciers in this area are medium-sized glaciers with a large elevation range (Reference HewittHewitt, 2007).
To date, glacier surges in the Karakoram have been observed either by repeated visits to the area and subsequent descriptions of the geometric changes of the glacier tongues and superficial moraines (Reference HewittHewitt, 2007), or the investigation of photographs from different dates (Reference Diolaiuti, Pecci and SmiragliaDiolaiuti and others, 2003). No information exists at present concerning the temporal evolution of the dynamic state of a surging glacier in the Karakoram.
North Gasherbrum Glacier
The region of the highest mountains around Baltoro Glacier, including the second-highest peak in the world, K2, was visited by European travellers as early as the 19th century. Members of a large expedition mapped this glacier and the surrounding mountains in detail in 1909 (Reference De FilippiDe Filippi, 1912). In contrast, the eastern part of this mountain group is much less known. To the east of Concordia, the main junction of Baltoro Glacier, the Gasherbrum and Broad Peak group comprises about a dozen peaks, with three reaching >8000 m a.s.l. To the north and east of these massifs is the drainage basin of North Gasherbrum Glacier, while the southern part of the Gasherbrum group contributes to the much smaller South Gasherbrum Glacier, a tributary of Baltoro Glacier. North Gasherbrum Glacier is fed by avalanches originating from the peaks of Gasherbrum II (8035 m a.s.l.) and Gasherbrum IV (7925 m a.s.l.) in the south and Broad Peak (8047 m a.s.l.) in the west of the main accumulation basin (Fig. 1). The slopes of the mountain walls surrounding the high-accumulation basins are steeper than 40° in many cases. The glacier is ∼23 km long, covers an area of ∼67 km2 and the glacier tongue is formed by four tributaries of different sizes. Besides the main trunk of the glacier, originating from a basin between Broad Peak, Gasherbrum IV and Gasherbrum II, only the lower, southern tributary seems to contribute a considerable amount of ice to the glacier. The accumulation basin of the main trunk is situated between 5400 and 6000 m a.s.l., with some parts reaching ∼6600 m a.s.l. However, the very steep slopes around the accumulation basins probably provide a considerable part of the total accumulation as avalanches. Ice in the ablation area of the glacier is mainly flowing eastward and later northeastward from ∼5000 m a.s.l. until the glacier terminus at 4230 m a.s.l. This part of the glacier is ∼15 km long, which results in a mean surface gradient of 0.05. The entire glacier shows an elevation difference of >3800 m from the highest peak to the glacier tongue, and thus the mean slope is ∼3× the slope of the glacier tongue. Due to the high fraction of avalanche snow in the accumulation basins, originating from elevations above 6500 m, it can be expected that the glacier is rather cold above the ablation zone. Only the lower tongue receives a large amount of meltwater, which will shift this part of the glacier into temperate conditions.
Satellite view (Landsat-7 Enhanced Thematic Mapper Plus (ETM+) image of August 2002) of North Gasherbrum Glacier in the Karakoram. The glacier flows towards the top right, where it terminates in a calving front at Shaksgam River. The flowline indicated on the image is used for the extraction of the longitudinal velocity profiles.
During a short visit in 2006 the lower part of the tongue was visually investigated. The glacier tongue shows supraglacial debris cover of ∼50%, which mainly originates from the confluence between the main trunk and a southern tributary 10 km from the snout, at ∼4700 m elevation. The northern half of the tongue is only partly covered by debris, but shows an extremely rough surface, with a very high density of ice sails, sharp ice pinnacles rising from the glacier surface up to 15–20 m. There is no indication of large velocity gradients between the main trunk and the inflowing tributaries. On the eastward-flowing lower glacier tongue, no shear zones can be detected between the individual glacier parts, indicating that the main trunk dominates the dynamical conditions of this part of the glacier. The snout of the glacier ends in the Shaksgam valley, a tributary of the Yarkand valley.
The Shaksgam valley has an area of 8200 km2, of which 32% is covered by glaciers and perennial snowfields. The long-term equilibrium-line altitude (ELA) is situated between 5100 and 5400 m elevation (Reference FengFeng, 1991). The main discharge from snow and glacier melt occurs between June and September, when 80% of the total discharge is generated. The snout position of North Gasherbrum Glacier is known to vary between ∼500 m from the right valley wall (Reference FengFeng, 1991) and touching the wall (Reference DesioDesio, 1930, Reference Desio1955). The water from these rivers drains via the Yarkand River into the Tarim basin, one of the largest inland depressions in the world. As the distance between the glacier terminus and the opposite valley wall is just a few hundred metres, even minor advances of the glacier could block the Shaksgam River and produce an ice-dammed lake. The glaciers in the upper Shaksgam valley are known to produce glacier-dammed lakes. So far the damming of a lake by North Gasherbrum Glacier has not been not observed, but the upstream Kyagar and Tram Kangri Glaciers have a long history of destructive glacier lake outburst floods (Reference FengFeng, 1991) which could affect almost two million people and widespread infrastructure downstream (Reference Hewitt and LiuHewitt and Liu, 2010). Such outburst floods from lateral glaciers blocking the main valley occur in both the Indus and Yarkand basins (Reference Feng, Liu, Xiangsong and YuchaoFeng and Liu, 1990), and >20 glaciers have been identified where this happened in the 20th century (Reference HewittHewitt, 1982; Reference ZhangZhang, 1992).
North Gasherbrum Glacier is classified as a surge-type glacier by Reference KotlyakovKotlyakov (1997). But so far no surge has been observed and there has not been any clear sign of former surge activity. This is in contrast to its southeastern neighbour, Urdok Glacier, which shows a series of looped moraines, originating from its northern glacier branch, and clearly indicates at least five surge events of this tributary during the past. While investigating glacier changes in the central Karakoram on the basis of remote-sensing imagery, it was observed that North Gasherbrum Glacier showed non-typical surface displacements over a certain time period. This observation led to a more detailed analysis of the temporal evolution of surface changes on this glacier; the neighbouring glaciers did not show any anomalies in the same time period.
The region east of the highest peaks in the main ridge is bounded by the Karakoram fault, a major geologic feature running for >750 km from the Pamirs to the Kailash region (Reference Searle and PhillipsSearle and Phillips, 2007). North Gasherbrum Glacier flows across this fault, and the middle part of the tongue (∼4500 m a.s.l.) is situated directly above the fault. In consequence the mountains in the accumulation zone of the glacier consist of Gasherbrum–Broad Peak quartz diorite, which intruded the limestones of the Shaksgam and Aghil formations. Thus the erosion rates in the higher regions can be expected to be low (Reference Searle and PhillipsSearle and Phillips, 2007), while the lower parts of the glacier might receive a considerable amount of debris from deformed limestone.
Glacier surges
Glacier surges are periodic advances of a glacier, often taking the form of rapid motion which leads to a sudden advance of the glacier snout. A typical example is Variegated Glacier, Alaska, whose advance in 1982–83 was well documented by Reference KambKamb and others (1985). Largely based on that work, the common perception of a surging glacier was that it always involved a rapid advance of the terminus, associated with high velocities. However, in their classic early survey, Reference Meier and PostMeier and Post (1969) advised that surges came in all shapes and sizes, a fact subsequently borne out by the detailed observations of Reference Frappé-Sénéclauze and ClarkeFrappé-Sénéclauze and Clarke (2007) on Trapridge Glacier, Yukon, Canada, and by Reference Murray, Dowdeswell, Drewry and FrearsonMurray and others (1998) on Bakaninbreen, Svalbard, both of whose most recent surges have been moderate. Surging glaciers occur throughout the world, although few have been studied in the Himalaya. In the European Alps, the last known surging glacier, Vernagtferner, lost its oscillatory behaviour after its 1900 advance due, presumably, to declining mass balance promoted by warming since the 1800s (Reference HoinkesHoinkes, 1969). In the Pamirs, where a large number of potential surge-type glaciers are identified, a very detailed observation exists for Medvezhiy glacier (Reference Kotlyakov, Osipova and TsvetkovKotlyakov and others, 1997). Here maximum velocities during the surge in 1989 were observed to be 70 cm d−1 ∼1 km upstream of the advancing glacier tongue.
Data and Image Analysis
Our analysis of the dynamic state of North Gasherbrum Glacier is based on Landsat-7 Enhanced Thematic Mapper Plus (ETM+) data from 2001–09. Landsat-7 provides multispectral images with different spatial resolutions, where the panchromatic band 8 provides the smallest pixel size of 15 m. Landsat-7 ETM+ data cover the period since 1999 and, due to a repeat coverage of 16 days, there is a good possibility of compiling multi-year time series of cloud-free imagery for a specific region. We used feature tracking (Reference Scambos, Dutkiewicz, Wilson and BindschadlerScambos and others, 1992; Reference Strozzi, Luckman, Murray, Wegmuller and WernerStrozzi and others, 2002) on a collection of repeated Landsat-7 ETM+ imagery to measure seasonal and interannual ice-flow rates for the period 2001–09. This method aims to find similar patterns (e.g. groups of crevasses advecting down-glacier, or specific supraglacial debris patterns) between two consecutive images. The principles of this method are based on a cross-correlation algorithm used on image patches of repeat-pass satellite image pairs (Reference Strozzi, Luckman, Murray, Wegmuller and WernerStrozzi and others, 2002). The method relies on first matching nonmoving features within the two consecutive images (such as rock outcrops), in order to find the geometric solution to co-register the image pair with subpixel accuracy. This was done by manually identifying tie points in the corresponding image pairs. To minimize any geometric distortions, only images of the same path/row were used in this study (path 148, row 35). North Gasherbrum Glacier is situated in the centre of this footprint, enabling us to use Landsat-7 ETM+ acquisitions, even after the scan-line correction failure in May 2003, providing full coverage over the entire glacier catchment (Fig. 1). Once the geometry matching had been achieved, patches of moving pixels (glacier ice, supraglacial debris) within the first (earlier) image were correlated within larger windows in the second (later) image to measure their displacement relative to the fixed features, giving the mean surface velocity (m d−1) between the two acquisition dates. We only used the high-resolution (15 m) band 8 for velocity extraction, and a total of 16 image pairs were processed (Table 1). The error in the derived surface displacement is primarily determined by the accuracy of the co-registration (<15 m) between the repeat-pass images and the temporal separation between the image acquisitions. Estimated errors for the analysed image pairs vary between ∼0.03 and ∼0.2 m d−1, depending on the time interval between acquisitions. (Due to a larger image separation for longer time-spans and a fixed location accuracy, the error decreases for an increasing observation period.)
Information about the Landsat-7 ETM+ acquisitions, the image pairs used and their periods
A large change in the dynamic state of the glacier was detected while inspecting a series of Landsat-7 scenes for general glacier evolution. Closer investigation revealed a sequential acceleration and slowdown of the glacier which is discussed here in more detail. The analysis focuses on a central flowline on the tongue of North Gasherbrum Glacier (Fig. 1). Along this flowline, surface displacements from the feature-tracking results were extracted. Examples for the surface displacement fields are shown in Figure 2. In general, valid displacements can be derived for the main glacier tongue below ∼5000 m elevation. The comparison of the quiescent phase with the main surge period shows that the lower tributaries are not involved in the surge activity. Usually successful tracking results provided surface displacements rather regularly along the lower part of the flowline, with only small gaps of several tens of metres. The calculated surface velocities could therefore be interpolated to provide continuous profiles. The surface geometry of the flowline shown in Figure 3 is based on the flowline location from Landsat-7 images and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation model (GDEM) (ASTER GDEM is a product of the Ministry of Economy, Trade and Industry, Japan (METI) and NASA) for elevation information. Along the profile length of ∼21 km, the glacier surface descends from ∼5800 to 4230 m a.s.l. The part in the accumulation zone is displayed in grey, while the ablation zone is shown in dark blue. The confluence of tributaries with the main glacier is shown as light blue lines at the bottom.
Surface displacements on North Gasherbrum Glacier resulting from the application of the feature-tracking algorithm on the Landsat-7 ETM+ scenes. Valid surface displacements are usually restricted to the main glacier tongue below 5000 m a.s.l. Here examples are given for (a) the quiescent stage, (b) the beginning of the surge, (c) the main phase of the surge and (d) the slowdown phase. Dates are dd/mm/yy.
Surface topography along the central flowline of North Gasherbrum Glacier (Fig. 1). The grey part of the curve indicates the accumulation zone; the grey lines above the distance axis indicate the inflow of tributary branches.
The advance, 2003–07
Figure 4 shows velocity profiles from the years 2003–07. The velocities are calculated in terms of the available images, but each image pair covers ∼1 year, so the results represent annual mean velocities. For 2003/04 we only have valid results for a small segment between ∼8 and 10 km along the profile. The remaining profiles show a steady growth in velocity from the minimum in 2002/03 (0.38 m d−1 or 139 m a−1 at 6.5 km distance) to the maximum at 2006/07 (1.1 m d−1 or 402 m a−1 at 18.2 km distance). The feature-tracking results terminate at increasing distances downstream (the snout is relatively stationary near 21 km), which we take to represent the boundary between active and relatively stagnant ice. The increasing velocity leads to an enhanced transfer of ice downstream, and the glacier surface rises by ∼30 m in its downstream reaches. The magnitude of this increase in surface elevation was estimated from the comparison of photographs of the glacier taken in 2006 and the extent of the glacier at the left lateral moraine in the satellite image of 2007.
The advance phase of the glacier, 2003–07. The figure plots the velocity (measured as annual average) as a function of distance downstream in the years 2003, 2004, 2005, 2006 and 2007. The velocity rises from a linearly decreasing function of distance in 2003 to a linearly increasing function of distance in 2006, with a shock wavefront at 15.5 km. By 2007 the wave has advanced a further 3 km.
The retreat, 2007–09, 2002
Between 2007 (when the annual velocity peaks) and 2009 the velocity decreases. Figure 5 shows this decline, as well as the profile from 2001/02, which, since it is higher than the subsequent 2002/03 profile, suggests that it is also part of a declining phase, and thus that the behaviour may well be periodic. The unknown number of years between the 2009 profile and a subsequent 2002-like profile suggests a period of at least 8 years and probably somewhat longer. Taking into account the active phase of the surge, the typical period is probably about 15 years.
The retreat phase of the glacier, 2007–09 and 2002. Velocity profiles are plotted as in Figure 4. Because the 2002 profile lies above the 2003 profile, it suggests that 2002 was still in a retreat phase, and may be similar to a profile beyond 2010. If we hypothesize that the glacier surges periodically, this suggests that the interval between the 2009 and 2002 profiles represents the quiescent phase, and thus that the surge period is liable to be ∼15 years.
Summer 2009
While the velocity profiles in Figures 4 and 5 are relatively smooth, the one shown in Figure 6 has noticeable shortwave variation on a length scale of 1 km. This profile represents the average velocity over the late summer months August and September and the transition to the beginning of winter in early December 2009. In the first half of the period, strong surface melt can be expected, while in the second half water at the base of the glacier should still be available from the preceding melt period. Thus, the sliding velocity can be expected to be large. A similar condition was also observed for the summer period (June–October) in 2008, but only a part of the profile could be analysed due to the quality of the feature-tracing results.
The summer velocity profile from August to December 2009. Note this short-term velocity profile has short-wavelength (1 km) variations, which are absent in the longer-term average profiles.
Frontal fluctuations
Glacier surges are usually thought to be connected with rapid advances of the glacier front. In the case of North Gasherbrum Glacier there is not much room for an advance. Between the glacier front and the opposite valley slope there was only a distance of 200–450 m in 2003, the year with the minimum velocity conditions. In 2007 this distance reduced to ∼30–40 m in the centre and still ∼400 m at the glacier margins. The maximum advance is observed in 2008, when the glacier at some locations seems almost to have made contact with the eastern valley slope. Compared to the total length of the glacier, this frontal advance of a maximum of ∼200 m is rather small. There are also indications of possible glacier fluctuations from earlier remote-sensing information. Figure 7 shows the glacier snout in 1971, 1980 (both from Corona imagery) and 2007. In 1971 the minimum distance was ∼80 m between the glacier and the eastern valley slope. The small channel along the centre of the glacier snout is ∼90 m wide. In 2007 a similar channel exists, with a mean width of ∼30–40 m. In 1980 the mean distance at the same location was ∼200 m. If a periodic fluctuation of the glacier is assumed, the state in 1971 is probably very close to a maximum, while conditions in 1980 are similar to those in 2002/03, which are assumed to show a minimum state. This suggests that the periodicity of the fluctuations is roughly a decade, if we assume that the behaviour is indeed periodic (we emphasize that there is no direct evidence for this).
Frontal positions of North Gasherbrum Glacier in 1971, 1980 and 2007 (the year with the maximum advance). The comparison is based on imagery of the Corona mission and Landsat-7 ETM+.
Interpretation of the Data
The velocity wave
Figure 4 represents a transition from a pre-surge velocity profile in 2002/03, which is more or less linearly declining from 9 to 16 km, to a surge velocity profile in 2005/06, which takes the form of a linearly increasing velocity from 10 to 16 km. Between 2006 and 2007, this profile advances as a wave, with the peak in the velocity advancing from 15.5 to 18.5 km, so the wave speed is ∼3 km a−1 . The declining part of the profile represents the diffusive structure of a shock wave with width ∼3 km.
Seasonal modulation
If we look closely at the 2006/07 profile in Figure 4, we can see that the average profile is modulated by a regular sequence of small-amplitude waves, with peaks at ∼14, 16.5, 18.5, 20 and 20.8 km, so the wavelength is in the region of 2–2.5 km. Once we have recognized these modulations, we can see similar modulations in 2005/06 (at 12, 14, 15.5, 16.5 and 17.5 km) and even in 2004/05 (at 10.5, 12, 13.5, 15 and 16 km). Assuming these modulated peaks correspond to each other, they represent small-amplitude travelling modulation waves with speeds in the range 1.5– 3.3 km a−1 .
Short waves
Figure 6, the only short time period profile, shows distinct variation on a short wavelength of ≲1 km, which is absent on the longer time mean velocity profiles. Assuming these features are real, they might either represent dynamic wave features associated with compressional waves, or they might represent variations associated with underlying bed-slope variations in the form of riegels.
Upstream slowdown
All the velocity profiles exhibited a marked slowdown in an abrupt step, located (for 2002) at 10 km. The location of this step is variable; in the years ordered 2003, 2004, …, 2009, 2002, it is located approximately at 9, 10, 10, 10, 12, 12, 11, 10 km; it thus appears to vary dynamically with the progression of the wave, advancing abruptly between 2006 and 2007 while the velocity wave propagates forward. One might suppose that the step is associated with one of the tributaries evident in Figures 1 and 3, located at ∼10.1 and ∼12.3 km along the flowline, and these do appear to be close to the range of the step, although the 2003 step is upstream of both. Alternatively, the slowdown could be caused by a steeper slope caused either by bed topography (but then it presumably should not move) or by drawdown of the ice during the surge (but then one expects upstream propagation of the step during the surge rather than an advance). Large variations in the surface slope are, however, not observed in the flowline geometry (Fig. 3).
A Simple Model
When large changes occur on a glacier, such as described here, there are three obvious ways to interpret them. (1) They could be a transient response (e.g. as a kinematic wave) to a prolonged change in external conditions (e.g. a prolonged series of high-accumulation years). (2) They could represent a threshold phenomenon of a stable but excitable system. (3) They could represent a self-sustained oscillation, as in a surging glacier. The latter two are the easiest to understand, and differ only in the occurrence of periodicity. It is therefore natural to consider the large disturbance to North Gasherbrum Glacier in this way.
In our interpretation of the data as that of a surging glacier, we focus on those aspects of a glacier which are thought to be instrumental in determining such behaviour. The oscillatory motion itself involves variations in sliding, which implies the presence of a subglacial drainage system, which can operate in different modes. The mechanism of such switching might be through freeze/thaw cycles (Reference MacAyealMacAyeal, 1993), hydraulic runaway (Reference Fowler and SchiaviFowler and Schiavi, 1998), Röthlisberger (R)-channel shutdown (Reference KambKamb and others, 1985) or some other mechanism. Distinguishing between different drainage mechanisms requires knowledge of the sliding behaviour of the glacier as well as its plumbing, neither of which are available without detailed fieldwork, such as that done on Variegated Glacier during its 1982 surge. The other principal external influence is the mass balance. It is, for example, likely that the demise of surging glaciers in the European Alps is due to declining accumulation rates. Statistical studies (e.g. Reference Clarke, Schmok, Ommanney and CollinsClarke and others, 1986) are of little critical use here, since they cannot unfold the complex interplay between the different parameters describing the ice flow, sliding and hydrology.
Having presented the observational data and interpreted them in terms of large-scale wave motion, we now turn to the problem of understanding why and how such behaviour can occur. In order to do this, we describe a model of the ice motion which explicitly involves a basal sliding law together with a subglacial drainage theory. These ingredients are essential to the description of the variations in velocity which are seen, but nothing is known about them, so our model is based on a reasonable inference of sliding and drainage characteristics. In particular, we use a description of drainage switching associated with the transition from a low-pressure channel system to a high-pressure linked-cavity system. Our comparison of the theory with the observational results then provides a context for understanding, and points the way to what future observations will be necessary.
Rapid variations in glacier motion are most simply ascribed to alterations in basal sliding. Rapid changes in shearing could only be due to rapid changes in temperature (impossible by thermal conduction, but conceivable through large surface melting and subsequent englacial refreezing), or by rapid change in depth, which itself requires rapid change in velocity. Thus, in practice the only viable alternative to variations in sliding is by surface melting. A simple calculation using a constant activation energy, Q +, for creep above −10°C suggests that a surface velocity, u, due to shearing will be enhanced by Δu given approximately by
where the excess temperature, ΔT, generated by refreezing of a surface melt of Δh ice equivalent, redistributed over a glacier of depth d i, is given by
where L is latent heat and cp is specific heat. Using values Q + = 115 kJ mol−1 (Reference Cuffey and PatersonCuffey and Paterson, 2010), R = 8.3 J mol−1 K−1, T 0 = 273 K and L/c p = 160 K, we find
To obtain the tenfold increase in velocity seen in the present data over a period of 4 years, we would thus need a melt of Δh ≈ d i/30 over each of the four years. If d i = 200 m, for example, this represents annual melting in excess of 6 m, which seems excessive. We therefore suppose that the velocity variations are due to sliding anomalies. If this is the case, it is sufficient to consider a model in which the velocity, u, downstream is entirely by sliding, and this we do.
The sliding law
The sliding law relates the basal stress, τ, to the velocity, u, and the effective pressure, N, which is the ice overburden pressure minus the water pressure in the subglacial drainage system. Reference LliboutryLliboutry (1968, Reference Lliboutry1979) was an early advocate of the effects of subglacial cavitation, which led him to pose general laws of the form
in which n is the exponent in Glen’s flow law. Lliboutry found that f could be a unimodal (one-humped) function. Reference IkenIken (1981) suggested a simple reasoning why this should be so. The result was subsequently confirmed in more elaborate calculations by Reference FowlerFowler (1986) and Reference SchoofSchoof (2005). (The calculations are for Newtonian ice, n = 1, but one can infer the form of Equation (4) on dimensional grounds, assuming the bed roughness is sufficiently large that regelation can be ignored.)
An issue arises concerning the decreasing part of the sliding function, f. Reference FowlerFowler (1987a) suggests that for real beds, with roughness at arbitrary large scales, f would continue to increase. Reference SchoofSchoof’s (2005) more precise calculations suggested that this could be an oversimplification, although he finally suggested a sliding law in which f increases towards a final maximum, thus corresponding to a frictional resistance, also advocated earlier by Lliboutry. In this paper we take the view that the presence of a real maximum in bed stress is associated with the catastrophic behaviour of hanging glaciers (Reference FlotronFlotron, 1977; Reference Pralong and FunkPralong and Funk, 2006; Reference Faillettaz, Pralong, Funk and DeichmannFaillettaz and others, 2008), and that where such behaviour does not occur, then either the bed is sufficiently rough that f increases, or else it effectively does so through the effects of side-wall friction (an extreme example of this is the behaviour of Ice Stream B, West Antarctica, where the driving stress exceeds the basal stress, and the balance is taken up by the lateral stresses; Reference Whillans and van der VeenWhillans and Van der Veen 1997). If we suppose a power law, then the sliding law takes the form
where c and r are constant and we would suppose r < 1/n. The form of this law receives some support from Reference Budd, Keage and BlundyBudd and others (1979), Reference BindschadlerBindschadler (1983) and Reference Iken and BindschadlerIken and Bindschadler (1986). The shear stress is adequately represented by the expression (Reference NyeNye, 1952)
where H is ice depth, x is downstream distance, ρ i is ice density, g is acceleration due to gravity and S is the mean valley slope; this expression neglects a further correction due to longitudinal stress gradients (Reference FowlerFowler, 1982).
Drainage
The classical theory of subglacial drainage is due to Reference RöthlisbergerRöthlisberger (1972), who described water flow through a channelized system by means of a set of equations which ultimately relate the effective pressure, N, to the water flux, Q. In the steady state, Reference FowlerFowler (1987b) showed that the solution could be adequately approximated, away from the glacier snout, by the algebraic relation
where n is the exponent in Glen’s flow law. This counterintuitive result (flow rate is less at higher water pressure) is responsible for the idea that such a channelized system will consist of one, or at best a few, channels, because for two slightly dissimilar neighbouring channels, the smaller (at higher pressure) will drain into the larger. This property, that ∂N/∂Q > 0, thus causes coarsening.
Reference LliboutryLliboutry (1968) envisaged drainage through a system of linked cavities and Reference KambKamb (1987) devised a theory for such a flow, which importantly showed that water pressure increased with volume flux, so ∂N/∂Q < 0 and the system is stable as a distributed system (as it must be, since cavities are distributed). In fact, Lliboutry’s earlier theory can be used to illuminate the two types of drainage system. Lliboutry identified a shadowing function, s, which we define to be the area fraction of the bed which is uncavitated. In keeping with the sliding theory discussed earlier, s is a monotonically decreasing function, s(Λ), with s(0) = 1. The effective pressure over the uncavitated part of the bed is then N/s, and in a linked-cavity system we suppose that flow from cavity to cavity goes through orifices which are simply little R channels. If there are n K such channels across the width of the glacier (essentially determined by the number of cavities), then the flow in each is Q/n K, so the analogue of Equation (7) is
Figure 8 illustrates the form of the drainage effective pressure given by Equation (8) relating effective pressure, N, to water flow, Q, through a field of linked cavities. To be specific, we define
so that C is a monotonically increasing (from zero at Λ = 0) function. We suppose that it saturates at large Λ, indicating that complete flotation does not occur. Note that when N is large, Λ is small, and Equation (8) reverts to Equation (7) with n K = 1. The functions used in Figure 8 are C(Λ) = (1−e −k Λ)/k and Λ = u/Nn , with k = 0.2, u = 1 and n = 3. For the curve labelled K n K = 200, while for the lower (R) curve n K = 1. The minimum of the curve (for a fixed value of n K) occurs at Λ = Λc ≈ 1/n when k is small. For Λ < Λc, to the right of the minimum, coarsening occurs, and a single (n K = 1) channel forms. To the left of the minimum, Λ > Λc, distributed drainage is stable, and this takes the form of linked cavities, K, for N > 0 (and n K ≫ 1). However, if N → 0, we suppose that some of the bed remains in contact with ice, so the shadowing function remains positive, and this allows a film flow, F, as the water flux decreases to zero. The three indicated components of the curve thus represent F: patchy film flow; K: linked cavities; R: R channel. Although very simple, this description has the merit of naturally representing the transition between a channel and a linked-cavity system which was inferred by observations on Variegated Glacier (Reference KambKamb and others, 1985), and also on Haut Glacier d′Arolla, Switzerland (Reference Nienow, Sharp and WillisNienow and others, 1998).
Drainage relation given by Equation (8). Starting at N = 0 the curve F increases sharply, reaches a maximum (not shown) and returns as the curve marked K.
Hydrological switching
Supposing Figure 8 is appropriate, it should be interpreted as follows. The minimum of the curve is at Λ = u/Nn = Λc, where Λc is constant. In terms of determining N, we must remember that Q is effectively prescribed, and u (and thus Λ) is determined by the sliding law. As argued by Reference FowlerFowler (1987b), the implicit dependence of N on u leads to a multivalued relationship, as shown in Figure 9. The values of N are presumed dependent on the water flux, as indicated in Figure 8, so N R is an increasing function of Q, while N K is a decreasing function of Q.
The transition between effective pressures in terms of the ice sliding velocity. The R-channel and Kamb (K) linked-cavity values depend on the water flux, Q, with the derivatives 
and 
.
The onset of a surge is associated with the gradual increase of depth during the quiescent phase. This leads to an increasing velocity and, if the accumulation rate is sufficiently large that u reaches the right-hand cusp in Figure 9, then the drainage system undergoes a transition to a linked-cavity system at lower effective pressure and the ice velocity increases. This is indicated in Figure 10, which shows how the multivalued drainage law leads to a multivalued relation between ice velocity and depth, even if the sliding law is single-valued, as in Equation (5).
The multivalued drainage law of Figure 9 leads to the multivalued relationship shown between u and H. The upper fast branch corresponds to the linked-cavity system, and the lower branch corresponds to the Röthlisberger system.
The data in Figure 4 suggest that the transition to fast flow occurs over a period of some years, and it is necessary to consider how fast we might expect the transition between systems to occur. Firstly, a scaling analysis of the Reference NyeNye (1976) version of Reference RöthlisbergerRöthlisberger’s (1972) model shows that the ice closure equation can be written in an approximate dimensional form (Reference FowlerFowler, 2009):
where f is a friction factor, ρ
w is the density of water, g is acceleration due to gravity, ρ
i is ice density, L is latent heat and K is a viscous closure rate coefficient. The timescale, t
0, of closure is approximately t
0 ∼ 1/KNn
, and is ∼1 day. Here, the cross-sectional area, S, is related to the water flux, Q, by S ∝ Q
3/4
. A further timescale occurs in the mass-conservation equation for water, but this is a factor 
smaller, where h
0 is the vertical descent over the water flow path. Typically ε ≪ 1, so the water flux is given quasi-statically by integrating the surface meltwater source term, M,
and the closure equation (Equation (10)) simply determines N. In the time-dependent oscillations which occur during jökulhlaups, N and S interact through the provision of an evolution equation for N at the channel entry. Normally this is at the portal from an ice-dammed or subglacial lake, but it could equally well describe the initial drainage from the surface through crevasses and moulins. At the beginning of the subglacial channel system, we would then prescribe an equation (cf. Reference FowlerFowler, 2009)
where A C is the effective crevasse/moulin surface area and M 0 is the surface meltwater supply rate to the channel inlet atrium. Choosing, somewhat imprecisely, A C = 100 m2, N ∼ 10 bar, Q ∼ 1 m3 s−1, this yields an inlet timescale:
Thus all the timescales associated with adjustment of the channel system are small compared with the annual timescale for meltwater supply variations. Consideration of the inter-cavity orifice transport in a linked-cavity system yields essentially similar results. A further timescale of relevance in the linked-cavity system is the adjustment time of the ice flow over cavities under changing effective pressures. This time is
where l C is the inter-cavity distance. For u = 100 m a−1 and l C = 20 m, this gives the longer timescale, t C ∼ 0.2 years. The timescale, t W, over which a local shutdown of the channel system will percolate across the width of a glacier might be expected to be somewhat longer than the ice-flow adjustment timescale, t C. Reference Hewitt and FowlerHewitt and Fowler (2008), following Reference WalderWalder (1986), also suggested a value t W ∼ 0.2 years.
We thus see that drainage transition is likely to occur on a timescale of weeks to months. This seems broadly consistent with observations of velocity and pressure variations over the course of a melt season (Reference KambKamb and others, 1985; Reference Nienow, Sharp and WillisNienow and others, 1998), and so we shall suppose that the transition time between branches in Figure 9 (and thus also in Fig. 10) is of this order. The fact that the annual velocity variation in Figure 4 occurs over a timescale of 3 years is not incommensurate with this, if we take into account that meltwater supply only occurs in the summer, so the drainage system has to rebuild itself each year. However, the main obvious cause of time variation in u arises from the fact that in a linked-cavity system u increases with water flux, Q, which itself increases with time during the melt season.
Surge model
The model of the surge of a glacier is based on that of Reference KambKamb and others (1985), as elaborated by Reference FowlerFowler (1987b). To the above relations between basal stress, drainage effective pressure and sliding velocity we add the equation of conservation of mass in the form
where the subscripts denote partial derivatives, a is the accumulation rate and we suppose for simplicity that all motion is by sliding. Our discussion of a surge follows that of Reference FowlerFowler (1987b, Reference Fowler1989) and is summarized here. As the glacier thickens during the pre-surge quiescent phase, H approaches H + in Figure 10, and when it reaches it, there is a transition of the drainage system to a distributed system on the part of the glacier (the reservoir region) where H > H −. This activated region shifts to the upper branch of Figure 10, and the resultant rapidly moving ice in the reservoir region slumps forward as a wave, and its height lowers as a consequence. This may or may not lead to an advance of the snout, depending on the volume of ice above the critical level, H −, in the reservoir area and the size of the downstream receiving area (Fig. 11).
Pre-surge (solid curve) and post-surge (dotted curve) surface profiles fortwo glaciers. In the upper profile, the difference between activation and deactivation elevations, ΔH = H + – H −, is sufficiently large that the surge leads to an advance of the snout, while in the lower, the snout is unaffected.
To this description we now add the effect of a water flux, Q, which increases downstream. Adopting Equation (6) in the form τ ≈H, then from Equation (5) the sliding law is
For Röthlisberger channels, we might suppose (Equation (7)) that the effective pressure, N, is determined by
while for the distributed drainage system, we might take for small (but not very small) values of N (Fig. 8) C ≈ Λ in Equation (9), (recall Λ = u/Nn ) and thus, from Equation (8),
and we approximate this by
From these, together with Equation (16), there follow the channel (R) and linked-cavity (K) sliding laws relating u to depth, H, and water flux, Q:
In the absence of specific information on water flow and sliding parameters, we assume generally that Q increases with downstream distance, x, and thus the R-channel sliding law has u decreasing with x, while the linked-cavity sliding law has u increasing with x. For illustrative purposes we take
This choice is partly motivated by Figure 4, where the 2002/03 and 2006/07 profiles suggest A R Hm = 0.1 × 104 m2 d−1 and A K Hm = 0.72 × 10−4 m2 d−1. We now interpret Figure 4 as follows. The reservoir area is between ∼10 and 15 km (central part of the ablation area) and between 2003 and 2006 the velocity rises from u R ∝ 1/x to the activated velocity, u K ∝ x, if the origin is taken at ∼4 km (roughly the position of the equilibrium line). This rise propagates somewhat further back, and also as soon as the velocity commences its ascent, a shock wave begins to propagate forward: already by 2005 it had advanced to ∼17 km.
Velocity wave
We can provide some information about the drawdown and the shock speed as follows. Based on Equation (21), H satisfies the approximate equation behind the surge front
during the surge (all the data are in the ablation zone, the equilibrium line being at ∼4 km). Observation of a net rise in elevation downstream (between 17 and 19 km) of ∼30 m suggests that
where as we surmise that the depth, d, of the glacier is ∼200 m. An estimate for the advance of the surge front, Δx f, can then be found by elementary geometrical considerations. If the reservoir area has length Δx and we assume a constant pre-surge surface slope, S, downstream of this, then we have the approximate relation
If we take Δx = 5 km, ΔH = 30 m, S = 0.04 (based on a 200 m drop over 5 km), then we calculate
This, obviously, is very crude, but it does provide an explanation for why the surge does not cause a real advance of the glacier.
Drawdown
We write Equation (22) in the form
If we put
(thus η is dimensionless and O(1)) then, since ΔH ≪ H −, the rate of drawdown is given by
where
assuming |a| ≪ GH −. From Figure 4, we can estimate G ≈ 2.63 × 10−2 a−1, so GH − ∼ 5 m a−1, compared with an ablation rate of 1–2 m a−1 . The timescale for drawdown of η by O(1), and thus the duration of the surge, is then
essentially as observed (the velocity increases from 2004 to 2007 and declines from 2007 to 2009).
Shock speed and width
The speed of the advancing shock front at x = x f, in which the velocity decreases from its surging value to essentially zero, is given by the jump condition
For values u − ∼ 1 m d−1 and ΔH/H − ∼ 0.15, this gives a wave speed of 2.4 km a−1, as observed. Moreover, the wave-front position accelerates, since by solving Equation (31) we have
and this acceleration can also be seen in Figure 4.
Shocks are discontinuities, but in reality are smoothed out by diffusive processes. Ordinary surface waves are smoothed by the term involving the deviation of the surface slope from the mean slope in Equation (6). In the present situation, the relaxation time, t W, for transition between differing drainage systems plays a key role. Suppose that the effective pressure satisfies an equation of the form:
In the shock, we put
where ΔX is the shock width; note that η and X are dimensionless. It follows from Equation (33) that the drainage transition occurs over the dimensionless range X = O(1) if we choose
where v
f = 
is the shock speed. If we take v
f = 2.4 km a−1 and t
W = 3 years, as suggested by Figure 4, then ΔX ∼ 7 km. This compares favourably to the apparent shock widths in Figure 4, which (by extrapolating the velocity in the front down to zero) take the apparent values 2, 3, 4.5 and 6 km in the years 2003, 2005, 2006 and 2007; the increase is associated with the increase in v
f as the front progresses.
To see that a shock structure does exist in the model, we write the mass-conservation equation in the form
where the function T(X) represents the transition from distributed to channelized drainage near the shock front. Ignoring the velocity in front of x f, we can take T to be a monotonically decreasing function of X, with T(−∞) = 1, T(∞) = 0. The slope-dependent term represents the normal surface diffusive effect present in Equation (6). Substituting Equation (34) into Equation (36) and retaining leading-order terms, we obtain an equation whose first integral is
with boundary conditions
We have used the definition of the shock speed given by Equation (31), and neglected ablation on the basis that (aΔX)/(v fΔH) ≪ 1, as is the case. The parameter δ is defined by
and for the values ΔH ∼ 30 m, S ∼ 0.04, ΔX ∼ 7 km, this gives δ ∼ 0.1, although it is 0.3 for the 2003 width of 2 km. For small δ, the solution is apparently η ≈ T(X), and this can be confirmed by consideration of Equation (37), where a graphical analysis of trajectories of w = η/T shows that a solution of Equations (37) and (38) exists for all positive δ.
Modulation waves
Finally, we turn to the seasonal modulation waves discussed above. Based on our inspection of Figure 4, we have analysed these data in more detail by looking at the variations of the velocity from the local spatial average. The results are shown in Figure 6, where four (or possibly five) clear amplitude peaks can be traced from 2005 to 2007. The waves move forward at an accelerating rate, just as the shock front does. The initial wavelength is ∼2.5 km.
The most obvious explanation for this variation is in the seasonal variation of water supply. Roughly speaking, the coefficient A K in Equation (21) (cf. Equation (20)) represents the meltwater supply. The simplest modification of Equation (21) to allow for this is to assume that A K is periodic, with a period of 1 year. Thus we reconsider the wave equation, Equation (26), in the form
where G is given by Equation (29) and we neglect ablation of order 1 m a−1 by comparison with GH ≈ 5 m a−1 . As before, we define η by Equation (27), and we define the dimensionless space-like scale
where l is a suitable length scale (e.g. the length of the ablation zone) Approximately, the surface perturbation, η, satisfies
where
We write
where 
is the time average of G and
gives the average drawdown of the surface (cf. Equation (28)). g′ is the time derivative of a (consequently) dimensionless periodic function, g, with period 1 year. Evidently φ should represent the surface perturbation profile that is shown in Figure 12.
The modulations of the velocity wave which is shown in Figure 4 are revealed in more detail here. The dates of the image acquisitions are identical with those in Figures 3 and 4, meaning that all three profiles represent annual time periods. A moving average for each of the years 2005, 2006 and 2007 is subtracted from the raw data, to reveal a sequence of waves. The lines indicate our interpretation of the movement of these. The motion reveals two phases: an initial transient while the velocity rises, followed by the subsequent propagation down-glacier at a speed of ∼2–3 km a−1. Since this inferred behaviour is so consistent with our theoretical discussion, it is worth pointing out that, historically, the data analysis was done first and the wave-propagation curves in this figure were added before the theoretical analysis was done.
The equation for φ is approximately (assuming small ν)
and a suitable boundary condition follows from the prescription that H = H− at the head of the reservoir region. In Figure 4, this is at ∼10 km, so if we take l = 10 km, then we put η = 0 and thus also φ = 0 at ξ = 0. The solution of Equation (46) is thus
At a fixed time, the surface is periodic in ξ with a period Δξ = G × 1 years ≈ 2.63 × 10−2, and thus the actual spatial period is Δx = (m + 1)xΔξ ≈ 0.26(m + 1) km, which corresponds to the observed period of 2.5 km if m ≈ 9 (if x = 10 km; m ≈ 5 if x = 15 km). This is a rather high value, but may be associated with our particular assumptions about the sliding and drainage laws. In particular, the exponent n in Equation (20) 2 is based on the sliding and drainage laws associated with viscous ice deformation over hard beds, whereas those associated with deformation of subglacial sediment may support a more plastic (high m) rheology (Reference KambKamb, 1991).
Additionally, the wave crests move at constant speed in the ξ plane, ξ = Gt, and thus x ∼ exp[(m + 1)Gt], consistent with the acceleration shown in Figure 12.
Conclusions
Glacier surges occur regularly in the Karakoram, even though most of these surges go unnoticed in remote areas. Nonetheless, such processes can pose a serious threat to humans by blocking valleys and impounding large amounts of water that can be suddenly released by a breach of the ice dam. North Gasherbrum Glacier clearly has the potential for blocking the Shaksgam valley, because even a small advance of a few hundred metres is sufficient. Therefore the observation that this glacier probably shows periodic fluctuations is somewhat disturbing.
Our study shows that the dynamic reaction of a glacier can be well analysed if a sufficient number of remote-sensing observations are available. Another prerequisite for such an analysis is the availability of an appropriate digital elevation model (DEM). For the rather smooth and large glacier tongue of North Gasherbrum Glacier the elevation information from the ASTER GDEM is sufficient for the analysis. In the case of smaller glaciers, the planned DEM based on the TanDEM-X mission is keenly awaited by the glaciological community.
We have provided an analytic interpretation of the remote-sensing imagery in terms of a postulated surge theory for North Gasherbrum Glacier, which includes the influence of basal water flux and the specific geometry on the changes in ice flux. Because very little is known of the glacier, it must be emphasized that this theory is hardly at the predictive stage. Instead, what we have tried to do is to interpret the data in terms of one particular theory of glacier dynamics, based on simple yet sophisticated concepts of glacier sliding with cavitation and subglacial hydrological switching. Using very simple quantitative arguments, we are able to explain in detail the form and propagation of the velocity wave, and even the modulative effects of annual meltwater supply variation. The only inference of our theory that may be at odds with expectation is the requirement of a drainage-switching timescale of 3 years, as opposed to the expected value of 0.2 years. This may be associated with the confounding effects of annual shutdown and restart, or it may be associated with a longer timescale associated with cross-glacier hydraulic transition. Equally, however, it may indicate that the premise of our theory is fundamentally flawed, and that the cause of the velocity fluctuations lies elsewhere. Even if sliding and drainage are responsible, the particular assumptions concerning the sliding law and hydraulic switching are very vague and unconstrained. At best, the present results only indicate a possible physical scenario for the observations.
It seems that, for the given amplitudes of velocity increase and the observed glacier geometry, massive advances of the glacier snout are not very likely. The ratio between reservoir and receiving area favours a rapid mass shift without a strong change in the snout position. Changes in the putative hydraulic system due to enhanced ice melt (if summer temperatures continue to rise), or changes in the ice geometry due to a higher mass flux from the accumulation area (if the higher winter precipitation, observed by Reference Archer and FowlerArcher and Fowler (2004) continues) could change the supposed surge characteristics of this glacier towards those of a glacier whose length oscillates periodically.
Acknowledgements
A.C.F. acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland mathematics initiative grant 06/MI/005. The free use of Landsat-7 and Corona data from the US Geological Survey for this analysis is greatly appreciated. C.M. and A.L. thank EvK2-CNR and the Geographic Institute ‘Ardito Desio’ of the University of Milan for generous support. We thank F. Ng and K. Hewitt for their thorough reviews and constructive suggestions.
