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 km² 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.

Fig. 1. 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 km², 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 km upstream of the advancing glacier tongue.

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⁻¹ or 139 m a⁻¹ at 6.5 km distance) to the maximum at 2006/07 (1.1 m d⁻¹ or 402 m a⁻¹ 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.

Fig. 4. 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.

Fig. 5. 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.

Fig. 6. 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).

Fig. 7. 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+.

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⁻¹ . 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⁻¹ .

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).

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

(4)

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

(5)

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)

(6)

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

(7)

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

(8)

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

(9)

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/Nⁿ , 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).

Fig. 8. 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/Nⁿ = Λ_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.

Fig. 9. 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).

Fig. 10. 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):

(10)

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 ₀, of closure is approximately t ₀ ∼ 1/KNⁿ , 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 ₀ 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,

(11)

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)

(12)

where A _C is the effective crevasse/moulin surface area and M ₀ is the surface meltwater supply rate to the channel inlet atrium. Choosing, somewhat imprecisely, A _C = 100 m², N ∼ 10 bar, Q ∼ 1 m³ s⁻¹, this yields an inlet timescale:

(13)

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

(14)

where l _C is the inter-cavity distance. For u = 100 m a⁻¹ 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

(15)

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).

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

(16)

For Röthlisberger channels, we might suppose (Equation (7)) that the effective pressure, N, is determined by

(17)

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/Nⁿ ) and thus, from Equation (8),

(18)

and we approximate this by

(19)

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:

(20)

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

(21)

This choice is partly motivated by Figure 4, where the 2002/03 and 2006/07 profiles suggest A _R H^m = 0.1 × 10⁴ m₂ d⁻¹ and A _K H^m = 0.72 × 10⁻⁴ m² d⁻¹. 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

(22)

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

(23)

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

(24)

If we take Δx = 5 km, ΔH = 30 m, S = 0.04 (based on a 200 m drop over 5 km), then we calculate

(25)

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

(26)

If we put

(27)

(thus η is dimensionless and O(1)) then, since ΔH ≪ H ₋, the rate of drawdown is given by

(28)

where

(29)

assuming |a| ≪ GH ₋. From Figure 4, we can estimate G ≈ 2.63 × 10⁻² a⁻¹, so GH ₋ ∼ 5 m a⁻¹, compared with an ablation rate of 1–2 m a⁻¹ . The timescale for drawdown of η by O(1), and thus the duration of the surge, is then

(30)

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

(31)

For values u ₋ ∼ 1 m d⁻¹ and ΔH/H ₋ ∼ 0.15, this gives a wave speed of 2.4 km a⁻¹, as observed. Moreover, the wave-front position accelerates, since by solving Equation (31) we have

(32)

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:

(33)

In the shock, we put

(34)

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

(35)

where v _f = is the shock speed. If we take v _f = 2.4 km a⁻¹ 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

(36)

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

(37)

with boundary conditions

(38)

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

(39)

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

(40)

where G is given by Equation (29) and we neglect ablation of order 1 m a⁻¹ by comparison with GH ≈ 5 m a⁻¹ . As before, we define η by Equation (27), and we define the dimensionless space-like scale

(41)

where l is a suitable length scale (e.g. the length of the ablation zone) Approximately, the surface perturbation, η, satisfies

(42)

where

(43)

We write

(44)

where is the time average of G and

(45)

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.

Fig. 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⁻¹. 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 ν)

(46)

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

(47)

At a fixed time, the surface is periodic in ξ with a period Δξ = G × 1 years ≈ 2.63 × 10⁻², 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) ₂ 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.