Flow regime

Our attempt to create a foundation for the mechanics of glacier surging is based on an assessment of the glacier flow regime which develops during a surge. As discussed by Reference PatersonKamb and others (1985), and others (summarized by Reference Raymond and ColbeckRaymond (1980)), surging involves substantial longitudinal compression and extension of the ice in the glacier. We have deduced from the published data that the longitudinal deviatoric stresses involved in this compression and extension of the ice are comparable to or substantially greater than the conventional estimate for the basal shear stress in the same glacier. Reference McMeeking and JohnsonMcMeeking and Johnson, 1985 have analyzed in some detail two-dimensional glacier mechanics when such conditions prevail and, independently, Reference Shoemaker and MorlandShoemaker and Morland (1984) have modeled one particular flow regime where longitudinal deviatoric stresses are significant. These two papers show that the conventional treatment of glacier mechanics (see, for example, Reference Hutter and ReidelHutter, [c1983]) is inappropriate for the flow regime that develops during surging. This is so because the conventional analysis is only valid where the longitudinal deviatoric stress is very small compared to the basal shear stress. This shows that a model for the stresses and flow regime during surging must be chosen with some care so that a proper foundation for the mechanics of surging can be utilized.

As a justification for the comments in the previous paragraph, we consider now the flow regime of some glaciers during surging. In Table I, we have calculated a strain-rate for each case of surging except for Variegated Glacier (1983). This calculation was done by taking the maximum velocity of the snout and dividing it by the length of the surging part of the glacier. This must be a lower bound to the maximum strain-rates that actually arise because the maximum velocity of the snout may occur due to stretching of a smaller length than the whole surging segment. In the case of Variegated Glacier (1983), Reference PatersonKamb and others (1985) have measured a local compressive strain-rate which we have quoted in Table I. Next, we have estimated the longitudinal stress required to produce the surging strain-rates by referring to data presented by Reference Frost and AshbyFrost and Ashby (1983). The estimates for stress taken to be twice the deviatoric stress s_xx should be valid near the upper surface where neither shear stress nor normal pressure will be very significant. Finally, we have calculated an estimate for the basal shear stress under the usual assumption that the down-slope component of weight is Table I supported by shear stress at the base. For this, we have used the average slope given by Reference Raymond and ColbeckRaymond (1980) and a depth of 200 m. The key point here, which can be seen in Table I, is that the longitudinal stress is usually somewhat larger than the conventional estimate for basal shear stress, In one case, (Variegated Glacier, 1983), the longitudinal stress locally is an order of magnitude greater than the shear stress. As a result, we must consider the analysis of glacier mechanics for situations where longitudinal stress is comparable to or larger than the basal shear stress if we are to understand glacial surging. As noted, the analysis of this situation has been discussed in separate papers by Reference McMeeking and JohnsonMcMeeking, 1985, and Reference Morland and JohnsonJohnson and Morland (1980).

The role of longitudinal stress during surging has been recognized by others such as Reference BuddBudd (1975). Starting with an equation which relates longitudinal stress gradients to the basal shear stress and the glacier weight, he postulated a dynamical model for surging due to variations in the basal shear stress. The model involves global averages of the basal shear stress and the physics is not entirely clear, although surges are generated. In addition, Reference LliboutryLliboutry (1968) and others have discussed surging due to sliding laws like (a) in Figure 1 and longitudinal stress plays a part, but the details of a mathematical model are not brought out.

Fig. 1 Sliding law − relationship between basal shear stress and basal velocity.

Compression front

Another feature which has been observed in some glacier surges is that a compression front moves down the glacier. This was observed by Reference PatersonKamb and others (1985) for the 1983 Variegated Glacier surge and by Reference Dolgushin and OsipovaDolgushin and Osipova (1975) for Medvezhy Glacier in 1963, in the form of a bulge. Reference PatersonKamb and others stated that below the front in Variegated Glacier the ice was almost stagnant and above it the flow was rapid. There were very high compressive strain-rates just behind the bulge and a corresponding thickening of the glacier. In the rapidly moving segment of the ice, 95% of the surface velocity was due to sliding and the remaining 5% was due to shear flow. Kamb and others attributed the rapid sliding to highly lubricated conditions at the base. They measured basal bore-hole water levels there which fluctuated around the ice-overburden value and concluded that the ice is supported mainly by a hydraulic system. They deduced that this gives rise to a very low level of drag at the base of the glacier.

The propagating compression front seems to be a discontinuity in the sliding conditions at the base of the glacier. Ahead of the front, the base of the ice is either stuck or sliding very slowly over the foundation. Behind the front, the base of the ice is relatively free of shear drag and moving rapidly compared to the ice beyond the front. The difference in velocities gives rise to a compressive region adjacent to the slowly moving ice.

We are not in a position to say why there should be nearly stuck portions and more freely sliding portions of the glacier. Reference PatersonKamb and others (1985) proposed that, as the surge front propagates down Variegated Glacier, the basal hydraulic system is disrupted and alters from one which causes high drag at the base to one which causes low drag. They thought that the surge front causes the closure of drainage tunnels and that this would flood the base of the glacier. Kamb and others suggested that the tunnel closure would be due to the surging mass up-stream from the front.

If the mechanism suggested by Reference PatersonKamb and others (1985) for Variegated Glacier is correct, we would view its operation as follows. The surge front separates the high-drag region from a low-drag region. In the low-drag region the down-slope component of weight of the overlying ice exceeds the drag on the sole. As a consequence, some of the down-slope weight is supported by the longitudinal stress. This is the only way that the large compression and extension rates can be developed. The stress thrusts upon the region of the glacier subject to high drag and as a result some of the weight of the surging mass is being supported by the slowly moving segment of the glacier. This force is transmitted to the base as a shear stress within about one glacier thickness from the basal discontinuity. The shear stress at the base in this transition region will be much larger than elsewhere, giving rise to a shear-stress concentration.

This stress concentration will cause viscous flow of the ice to be more rapid there than elsewhere. We suggest that the enhanced stresses at the discontinuity wilt allow the tunnel closure in a rapid fashion as suggested by Reference PatersonKamb and others (1985). Certainly, the high stresses will perturb the equilibrium of the basal water tunnel system. The non-linear interaction in the creep response between the enhanced shearing flow and the stresses generated by the overburden pressure around the tunnels would seem to suggest the tunnels would tend to shrink by creep. They may close entirely but, in any case, some water would be squeezed out.

The ice dynamics of a model for this kind of transition at the surge front would determine the rate at which the tunnels close. This, in turn, would determine the rate at which the surge front would propagate down the glacier. However, we believe that the processes involved in the transition would be controled by the stress field developed around the shear-stress discontinuity at the base. Thus the events can be characterized by the stress concentration at the discontinuity. That, in turn, can be determined neglecting the details of tunnel closure or any other relatively small-scale process involved in the transition. As a consequence, we concern ourselves only with the gross details of the stress field near the discontinuity with the understanding that these gross features of the stress field would control the fine details in a small region near the actual discontinuity. The fine details would in turn determine the propagation rate.

Basal condition discontinuity

With the above understanding, we are in a position to consider the mechanics of a basal shear discontinuity whatever the reason for its existence. For example, the discussion above refers to Variegated Glacier which is a temperate glacier. Reference Clarke, Clarke, Collins and ThompsonClark and others (1984) have discovered that there is a discontinuity at the base of Trapridge Glacier. This discontinuity is stationary and marks a transition from ice which is at the melting temperature in the upper part of the glacier and ice further down the ice stream which is frozen to the underlying stratum. It would appear that the frozen segment does not allow any basal sliding while the glacier is quiescent. However, the temperate segment slides and as a result a bulge builds up against the frozen dam. A possible surge mechanism is the failure of the frozen base due to the stress concentration or singularity at the discontinuity. This would occur probably when the stress reaches a critical level. Once the frozen base is ruptured near the discontinuity, part of the weight of the sliding ice would be thrust on to the next segment of frozen base. This, in turn, is likely to rupture and as a result a surge front would propagate down the glacier. Sliding more freely, the glacier would extend. As before, the stresses around the discontinuity would control the detailed rupture process. For the purposes of modeling, we can consider the overall stresses and characterize failure in terms of some critical stress level.

The picture is complicated in the case of Trapridge Glacier, because it is known to rest on deformable debris (Reference Clarke, Clarke, Collins and ThompsonClarke and others, 1984). The sliding part of the glacier may be moving relative to the debris or the deformation of the debris may be causing the motion. Presumably, the segment which is stuck overlies debris which is frozen. The frozen segment of the base would fail either by the glacier shearing away from the debris or by the frozen debris failing or breaking away from unfrozen debris underneath. The propagation of the slip surface accompanying the failure in each possibility would still be controled, we think, by the stress singularity or concentration at the now propagating stuck-sliding discontinuity. The modeling of this case of a temperate-sub-polar interface as in Trapridge Glacier can be taken care of by the same overall model as would be relevant to the purely temperate case of Variegated Glacier. We will discuss only an isothermal glacier model with the understanding that the sliding discontinuity at the base incorporated into our model may actually be a feature of non-isothermal conditions as in Trapridge Glacier. The model is meant to describe the general surge features of both an isothermal and a non-isothermal glacier, although certain detailed phenomena which distinguish these two cases are omitted. In particular, the treatment of ice rheology will be deficient as to the effects of temperature variations in the glacier, but we do not believe that plays a major role in controlling surging.

Hydrology

Basal hydrology is a significant feature of glacier behavior and in some cases of surging such as Variegated Glacier (Reference PatersonKamb and others, 1985). We will not treat the hydrology in our model, but it should be understood that water may be at the root of features that we describe in a phenomenological way. Thus the transition from nearly stuck to low-drag conditions at the base may be due to differences in the way that water invades the sole of the ice. One example of this would be the model of Reference PatersonKamb and others (1985) for the closure of drainage tunnels and the redistribution of water at the base into disconnected cavities. We will not attempt to characterize the drag at the base by modeling such features as described by Kamb and others. Instead, we will use sliding laws with features that would seem to have the right character for a variety of possible phenomena, including those driven by hydrological effects. Similarly, we are unable to characterize features of a purely hydrological nature but which may play essential parts in the surge cycle. For example, if the stuck-sliding discontinuity is acting as a dam to basal water or as a severe constriction as envisaged by Reference PatersonKamb and others (1985), then the restrained water would be released when the discontinuity nears the snout. As it is likely that this sort of phenomenon is responsible for the termination of surging in some cases, it is clearly important. However, we can only model this behavior in our approach by asserting that surging ceases at some point. Furthermore, the conditions for termination of a surge may be due to any one of a number of reasons such as that given above. One possibility is that the supply of melt water dries up seasonally and lubrication is lost. Another is that the glacier stretches and thins to such an extent that the driving force for the fast sliding falls below the critical level for sustenance and a return to slow sliding or nearly stuck conditions occurs. Due to our neglect of such issues, our surging model is quite incomplete but nevertheless brings out several salient features of the mechanics.

Basal sliding law

Hydrological, thermal, and basal deformation effects such as those mentioned above will be modeled in a sliding law chosen to display characteristics of many of the phenomena discussed above. These laws have shapes shown in Figure 1 where the shear drag on the base is plotted against the basal sliding velocity. Reference LliboutryLliboutry (1968) suggested that surging arose as a result of double-valued relationships between sliding velocity and drag of the kind marked (a) in Figure 1. Lliboutry based his law on an analysis of sliding with cavitation.

In addition to the above possibility, we propose that the suggestion by Reference PatersonKamb and others (1985) that there is a disruption of the basal drainage system near the surge front would lead to a sliding law like that marked (b) in Figure 1. Our thoughts are that the large stresses associated with sliding at velocity u ¹ _b are sufficient to close drainage tunnels and this leads to a drop in the drag and a corresponding increase in the sliding velocity. As long as the new drainage system is stable, the low-drag conditions will prevail. As mentioned previously, loss of basal water would lead to a return to high-drag conditions.

Another possibility, perhaps relevant to a partially frozen base, is that the maximum marked 1 in Figure 1 is the stress at which the base fails. Perhaps this would occur without any prior sliding (i.e. u ¹ _b = 0) and thereafter the failed base is capable of sustaining relatively low stresses while slip endures. The failure involved may occur in the ice, at the ice–base interface, or actually in the base in the case where the underlying stratum is deformable debris. To justify our assertion that a base of deformable debris may fail with the characteristics of curve (b) in Figure 1, we note the work of Reference RoscoeRoscoe (1970) on the failure of dense sand and similar granular materials such as glacial till.

We feel that sliding laws of the type shown in Figure 1 can be justified for a number of possible phenomena, although we have not attempted to characterize fully any particular case. Of course, the relevant behavior must be identified from field work before any particular case can be invoked. However, we have attempted to show that several possibilities fall within the purview of the sliding laws shown in Figure 1. There may be effects neglected in the sliding behavior and not present in the laws we consider, but we believe that the idealizations we have made retain the nature of the behavior relevant to several known surging glaciers. The absence of explicit thermal or hydro-logical effects from the models for surging is an idealization and a simplification helpful for the analysis. However, we will proceed with the understanding that the parameters of our sliding law may vary due to thermal or hydrological effects.

If we have omitted anything essential in the sliding law, then our surging model is invalid. It is more likely that we have simply omitted non-surging fluctuations which occur on a regular predictable basis, perhaps seasonal. Then all that occurs due to the fluctuations is that the critical condition to be met happens at one time of the year or cycle but cannot be met at other times of the year or cycle.

The overall concept is that slowly moving segments of the glacier are those sliding with a velocity up to u¹ _b (Fig. 1). This slow motion can occur up to drag levels which are high and approach τm. Rapidly sliding parts of the glacier are those moving faster at the base than u¹ _b. The drag involved would fall below τM at least for speeds greater than u¹ _b by only a limited amount. Our study of the mechanics involved in glacier surging has led us to believe that the initiation of surging is bound up with the transition from speeds below u¹ _b to those greater than u¹ _b. In passing, we note that we believe that phenomena bound up with sliding laws of the type shown in Figure 1 are essential for surging. If the sliding law is monotonic and never becomes otherwise but experiences seasonal fluctuations, then the motion of the glacier will be unsteady. As discussed by Reference FowlerFowler (1982), kinematic waves will be generated on the glacier, but it seems unlikely that any unstable surging behavior can result.

We have outlined how the down-slope component of the weight of the rapidly sliding region would thrust on to the more slowly moving segment, raising the basal shear stresses there. We envisage that these high stresses cause the mechanisms of transition to operate and so the rapidly sliding region would extend into previously slowly moving segments of the glacier. The peak drag T_m (see Fig. 1) would determine the thrust required for transition mechanisms to operate and so that would specify the conditions under which the sliding region would spread. Of Figure 2 course, this parameter would vary from point to point in the glacier. In addition, we would expect that there would be seasonal fluctuations due to melt-water flow or other effects.

Fig. 2. Schematic of a non-sliding glacier with a rapidly sliding segment AB.

Consider there to be a segment AB of the glacier in Figure 2. In this segment, the local peak drag τ _m would be exceeded by the down-slope component of the weight of the ice overhead. As a result, the down-slope component of the weight of AB would no longer be borne completely by the basal shear stress and the excess would have to be taken up by the longitudinal stress. Simultaneously, the segment would commence rapid sliding (the governing equations require large velocities when the longitudinal stress is significant; (see Reference McMeeking and JohnsonMcMeeking and Johnson, 1985). In turn, the longitudinal stress would thrust on to neighboring parts of the glacier outside AB and would tend to increase the shear stress at the base there. If the resulting basal shear stress is below the local τ _m, then the segment AB would be stable. If the shear stress at the base of the neighboring segments tends to exceed the local τ _m, then the region of rapid sliding would spread. This would involve a traveling compression front as discussed by Reference PatersonKamb and others (1985) and Reference PatersonPaterson (1981), or a traveling tensile front. As the sliding region spreads, increased weight would be thrust on to the slowly sliding parts, and so spreading would tend to continue.

The phenomenon of spreading would occur whether a law like (a) or (b) in Figure 1 were involved. In case (a) at a later stage, some of the rapidly sliding ice could be doing so with the down-slope weight borne by shear stresses at the base, that is, to the right of point 2 in Figure 1. This corresponds to what is called “continuously fast” sliding (Reference Raymond and ColbeckRaymond, 1980). On the other hand, if the law is like (b), some down-slope weight would always be supported by longitudinal stresses in the rapidly sliding region.

Idealization of the sliding law

Since we believe that the processes in a developing surge are bound up with behavior near the peak drag in laws like those shown in Figure 1, we will concern ourselves with modeling the sliding law around point 1 in Figure 1. At sliding speeds less u¹ _b than the down-slope component of the weight of the glacier is supported by the basal drag and the basal sliding speeds are relatively low. A substantial part of the ice flux is due to shearing flow. We will model this situation as one in which there is a no-slip condition at the base of the glacier. This modeling is possible because the no-slip condition gives rise to an ice flux solely due to shearing flow and a leading order basal shear stress which supports all of the down-slope component of glacier weight.

At sliding speeds greater than u¹ _b in Figure 1, the down-slope component of glacier weight will generally exceed the basal drag and most of the ice flux will be due to the basal velocity with very little due to shearing. We will model this situation as one of constant drag T_S at any sliding speed. This will give rise to an ice flux dominated by the basal velocity and some of the down-slope component of weight will be shifted from the basal drag to the longitudinal stress. This model of either no-slip or constant drag at the base will give rise to a very abrupt discontinuity at the surge front rather than a gradual one which would actually arise for laws like those shown in Figure 1. The nature of processes occurring at the gradual discontinuity will be preserved in the modeling and analysis of the abrupt discontinuity. At the gradual discontinuity there will be the large stresses of a stress concentration. At the abrupt discontinuity there will be the large stresses of a stress singularity. The kinematics of the transition in the abrupt discontinuity will be similar to the kinematics of the gradual transition.

It follows that the mechanics of such abrupt discontinuities are relevant to surging. The analysis of flow at such a boundary discontinuity has been considered by Hutter and Olunloyo (1980) but in the context of Newtonian viscosity. The transition involved was from no-slip to shear-free conditions. They found that a singularity in stress arises at the discontinuity as must be the case. In fact, the stress at such a discontinuity in any creeping solid or fluid must be singular. Rice and Rosengren (1968) and Reference HutchinsonHutchinson (1968[b]) (together denoted HRR) have analyzed such singularities for power-law materials including the special case of linear viscosity. Although the HRR analysis does not consider gravity-driven flows, we will show below that the analysis of the singular part of the stress near the discontinuity is the same whether the flow is gravity driven or not. The basic conclusion to be drawn is that there will be very large stresses developed at the discontinuity.

The infinite stresses will, of course, exceed any finite peak drag τ _m of the more realistic sliding laws shown in Figure 1. Our point of view is that the actual stresses near the surge front arising from laws like those in Figure 1 can be considered to be a perturbation on the singularity solution we have outlined. That is, the processes of transition from high drag to low drag would take place in a relatively small zone near the discontinuity and would smooth out the singularity. However, away from the transition zone, the stresses would be relatively unperturbed from those of the singularity solution and as a result characterized by parameters of that solution. These outer stresses in turn would determine the behavior in the transition zone, meaning that the overall transition would be characterized by the singularity solution. A critical value of a parameter which characterizes the singular solution would then control the process. This parameter which will be discussed in more detail later will be referred to as J. A somewhat similar situation is accepted as a rationale for the validity of non-linear fracture mechanics in structural analysis. We propose adopting the same approach here and so we will consider in our modeling and analysis only the abrupt transition and the resulting singularity.

The abrupt transition is used here as a model for the gradual transition arising from laws like those shown in Figure 1. Since we believe that the transition as described by Reference PatersonKamb and others (1985) for Variegated Glacier would give rise to a law like (b) in Figure 1, the abrupt discontinuity is therefore an approximate model for that case. However, the abrupt discontinuity and the resulting singularity may be more precise when used as a model for the frozen-unfrozen divide in glaciers such as Trapridge Glacier (Reference Clarke, Clarke, Collins and ThompsonClarke and others, 1984). An abrupt discontinuity of this type would give rise to a singular stress. In a relatively small region around the discontinuity, non-linear behavior of the ice, perhaps due to internal damage, would be activated. This would have the effect of eliminating the singularity. However, the effects would be present only very close to the discontinuity. The stresses farther afield would still be those of the singularity which would characterize all aspects of the discontinuity. Thus, the singularity model is one that we view as being relevant to basal discontinuities of a variety of types and we will use the model in the rest of the paper without reference to the character of the discontinuity.

We believe that it is the stress concentration at the surge front which causes the discontinuity to move down the glacier. The high stresses activate the process of transition at the discontinuity. The tunnel-closing mechanism suggested by Kamb and others for Variegated Glacier would be one example of a process of transition. For the case of Trapridge Glacier, we think that the frozen base must be ruptured somehow for the discontinuity to move and the processes of transition would be the rupturing mechanism. If the stress concentration is not sufficiently high or the high stresses are not spread over a large enough region, then the transition processes would not reach their critical extent and the discontinuity would not move. If, however, the stresses are high enough and spread out over a sufficiently large volume of ice, then the transition processes would change into those of lower drag. This would move the region of stress concentration and start to cause transition in the neighboring part of the glacier having high drag or a stuck base. If the critical conditions for transition can continue to be met, then the discontinuity will continue to move. Since we will be modeling the discontinuity as abrupt, the critical conditions will be phrased in terms of the parameter J which characterizes the singular stresses there. When J reaches a critical level, the discontinuity will move.

The effect would be like that which causes a slip surface or shear crack to propagate through over-consolidated clay and leads to the failure of embankments. A model for this, including the role of stress concentration or singularity in extending the failure, has been proposed by Palmer and Rice (1973). We will use certain aspects of their model and, in doing so, develop a framework for the mechanics of the propagating surge front.

Surging model

We are now in a position to give an overview of our surging model. From observations of Variegated Glacier by Reference PatersonKamb and others (1985), we see that the upper segment of the glacier was sliding rapidly for a couple of years before the surge year. This may have occurred because the critical stress for slow sliding in a law like those shown in Figure 1 was exceeded. The fast-sliding, low-drag region may have been confined originally to a small segment which later spread. It is possible that rapid sliding may occur locally in small sections of the glacier every season or it may be a permanent feature in certain areas. Our viewpoint is that it is very possible that such surging nuclei are common in many glaciers for a variety of reasons. In most years they are stabilized against growth by the prevailing state of the glacier. It should be remembered that, after a surge, the glacier is in a process of rebuilding and what may be stable one year may not be so in following years. Certainly, in Trapridge Glacier, the sliding region is present permanently in the temperate segment Reference Clarke, Clarke, Collins and Thompson(Clark and others, 1984) and the discontinuity at the temperate-sub-polar interface is stable for many years.

Rebuilding proceeds and eventually the critical conditions are met for the propagation of the discontinuity. The discontinuity may move slowly at first and then accelerate, as seems to be the case in Variegated Glacier (Reference PatersonKamb and others (1985)). Such time effects may be due to the natural response of the viscous flow of the ice or may be due to temporal influences of the hydrology. We do not take these effects into consideration but concern ourselves with the question of whether the pre-existing surge nuclei can be stable. When they start to propagate, we show that in common situations and according to our assumptions, the rapid-sliding region will tend to grow and capture the whole glacier.

Thus a surge is initiated and will continue until the rapid sliding is terminated. This may be due to the loss of lubricating ground water, as seems to be the case in Variegated Glacier (Reference PatersonKamb and others (1985)), or it may be because the glacier thins to such an extent that the driving shear stress is insufficient to sustain the state of sliding. This would correspond to the basal shear stress falling below the minimum or the asymptote beyond point 1 in Figure 1. The glacier would now return to its quiescent state and rebuilding would commence. Now, any surge nucleus would be stable for several years. This would be so because the stresses generated at the discontinuity involved would be below the critical level due to the new shape of the glacier. Thus, a full-scale new surge would not be initiated for several more years.

Models similar to the one outlined above have been used previously, such as in the work of Reference LliboutryLliboutry (1968). Similarly, Ramb and others (1985) have recognized the importance of sliding law non-linearities in their assessment of the surging of Variegated Glacier. However, we have considered the theoretical implications of the rapid sliding, the extending-compressing flow, and the discontinuities between slow- and fast-sliding regions. These aspects of the mechanics we have developed and incorporated into our model.

We now discuss the foundations of the mechanics which we believe must be understood for a complete modeling of surges. We will first discuss the stresses and velocities in the stuck parts (modeling slowly moving high-drag regions) and the sliding parts (modeling rapidly sliding relatively low-drag segments.) Following that we will consider the spreading of the sliding region and then a model surge, making comparisons with Medvezhy Glacier.

Non-sliding part

The leading order estimate of the shear stress in a stuck or non-sliding part of a glacier is

where ρ is the uniform ice density, g is the gravitational constant, and y _s(x) is the y position of the upper surface of the ice. The other stresses are

The longitudinal deviatoric stress s_xx is negligible. As discussed by Reference HutterHutter (1981) and others, the first correction to Equations (2) and (3) in a regular perturbation scheme would have magnitudes 0(δ) times the magnitudes in Equations (2) and (3), where δ is roughly the glacier thickness divided by the length and therefore small compared to unity. The velocity in a non-sliding part is

where h(x) = y _s(x) − y _b(X) is the glacier thickness. The velocity component v is 0(δ) times the magnitude of u because h′(x), y′_s(x), and y′_b(x) are all 0(δ) compared to unity.

In an aside, it should be noted that there is a boundary layer near the free surface where s_xx is more significant than σ_xy. This feature has been analyzed by Reference Johnson and McMeekingJohnson and McMeeking (1984) and its influence on the analysis of longitudinal deviatoric stress has been discussed further by Reference McMeeking and JohnsonMcMeeking and Johnson, 1985. However, the layer can be detected only when one seeks higher-order corrections to Equations (2), (3), and (4).

Sliding part

The solution for this segment has been studied by Reference McMeeking and JohnsonMcMeeking and Johnson, 1985. In this part there must be the possibility of substantial extensional and compressional flows. This will give rise to rapid sliding away from the basal slip discontinuity. In addition, the longitudinal stress must be large enough that it bears some of the down-slope component of the weight. The down-slope component of the weight is the component parallel to the base and the amount transferred to the longitudinal stress is the difference from the basal drag, i.e. ρghsin α – τ _s. In this part, the basal drag is not capable of adjusting itself to equate with the down-slope component of the weight. It is this fact which causes the longitudinal stresses to be large and determines the flow regime.

Reference McMeeking and JohnsonMcMeeking and Johnson, 1985 have found that the flow in this situation is to leading order a plug flow and so the longitudinal stress to leading order is independent of y,

In these circumstances, the shear stress can be found by integration of the momentum balance subject to the condition that the basal value equals the drag τ_S , so

As usual, the stress σ_yy balances the weight

but the longitudinal stress is an order of magnitude larger than σ _xy and σ _yy The boundary condition at y = y _s which requires the surface to be shear-free provides a differential equation for σ _xx ,

and we will simplify this by taking y _b constant. Thus

where h _A is h at x _A as in Figure 2.

The velocities are such that

and u will generally be much larger than the velocity in the non-sliding part. The drag acts against the direction of sliding. Lastly, the shape of the glacier, i.e. the ice thickness h = ys – yb, is determined from mass conservation or, equivalently, the kinematic boundary condition on the velocity at the surface. We reiterate that, if the stresses are not as outlined above, then there can be no large compressive flows and as a result no compressive front as observed in surges.

At the intersection of a sliding part and a non-sliding part there will be a region straddling the intersection, approximately one glacier thickness in length, in which neither the solution in the sliding nor in the non-sliding part is valid. Strictly, we should solve the governing equations there and then match that solution to the sliding region solution on one side and the non-sliding region solution on the other side. Instead, we will match the two sections in an ad hoc way. The leading order velocity in a sliding region is generally several orders of magnitude larger than the velocity in a stuck part. Consequently, it is appropriate to set the velocity in the sliding part equal to zero at the join. Matching the velocities at higher orders will account for the relatively weak shearing flow in the stuck part. We have not performed this higher-order matching but we believe it can be completed successfully. This means that in the sliding region

If the sliding region is contained between two non-sliding regions as in Figure 2, then

from which we can compute σxx (x _A) in Equation (9). Alternatively, if the sliding region has extended all the way to the snout at l then we can take, as an approximation, σ_xx = 0 there. Then

which means that the excess of the down-slope weight over τ _s(l – x) is supported by the longitudinal stress at x.

Generally, there will be a large longitudinal stress at the join between the sliding and non-sliding parts. This stress will be transmitted to the bed as a shear stress in the stuck part within approximately one glacier thickness from the join. In fact, the join between the sliding and non-sliding parts is a singularity for shear stress in the stuck region.

Stresses at the discontinuity

The singularity at the boundary-condition discontinuity is exactly of the type that arises in crack-tip fracture mechanics for structural materials Reference Rice and Rosengren(Rice and Rosengren, 1968; Reference HutchinsonHutchinson, 1968[b]). This singularity, known as the HRR singularity, dominates the near field and can be expressed as

where r and θ are polar coordinates from the point of discontinuity as shown in Figure 4 and the parameter J to be discussed below scales the magnitude of the stress field. The function has been computed for pure shear with n = 3 Figure 4 by Reference HutchinsonHutchinson (1968[a]). The pure shear crack-tip calculations of Hutchinson provide the correct near-tip field because anti-symmetry and incompressibility lead to stationary velocities along the axis of anti-symmetry ahead of the crack. Reference HutchinsonHutchinson (1968[a]) was considering the case of non-linear elasticity rather than viscous flow but his result can be carried over by interpreting displacements as velocities, etc. The parameter J which characterizes the stress field near the discontinuity is the J integral defined as a generalization of the form used by Reference Palmer, Rice and SerPalmer and Rice (1973)

where b is the uniform (gravitational) body force per unit volume,

where ∊ is the strain-rate and n is a unit vector on the contour ⌈ from D to E to 0 as shown in Figure 4. When the base of the glacier is a straight line, the integral is path-independent.

Fig. 4. Boundary-condition discontinuity with contour Γ for J-integral.

The stresses and strain-rates near the discontinuity are characterized by the single parameter J. It follows that any process of rupture or break-down that takes place very near the discontinuity will be controlled by J. We propose that the criterion for spread of the sliding region is that J equals or exceeds a critical value J _c. A similar criterion has been proposed for the growth of cracks in creeping metals at high temperature (Reference Landes and BegleyLandes and Begley, 1976). Similarly, the deepening of crevasses may be controlled by J and this will be discussed separately by McMeeking in future work. The path-independence of J will be demonstrated in that paper.

At the discontinuity, there will be infinite stresses, or in reality very high stresses. Such infinite or very high stresses alone are not considered to be sufficient to operate the transition mechanism and cause the discontinuity to move. This can be understood if one considers discontinuities being acted upon by very low thrust. The stresses will still be infinite or very high near the discontinuity but one concludes that the discontinuity is not likely to move. Instead, we deduce that it is necessary for high stresses to be present over a critically large volume of ice for transition to take place. As J increases, the high stresses near the singular point will spread over a larger volume of ice and eventually cover a sufficiently large volume. Depending on the situation, the processes involved may be the rupturing of ice or base (frozen base) or the creep closure of tunnels as suggested by Reference PatersonKamb and others (1985) for Variegated Glacier. Other mechanisms may operate in other cases. Thus, the critical value of J would be related to the processes involved and would depend on how large a region of large stresses is required for the transition to occur and for the discontinuity to move. We will use the parameter in a phenomenological manner and not attempt to relate the critical value of J to the processes involved.

This conludes our discussion of the foundations for the mechanics of surging. Further discussion of the issues of mechanics involved can be found in Reference McMeeking and JohnsonMcMeeking and Johnson (1985) for the stress analysis and in Reference HutchinsonHutchinson (1968[a], Reference Hutchinson[b]), Reference Rice and RosengrenRice and Rosengren (1968), and the wealth of literature on non-linear fracture mechanics. We will now proceed to a description of our model for surging.

Growth of the sliding region

During most of the time the majority of the glacier base will be non-sliding and the shear stress at the base according to Equation (2) will be

As discussed in the introduction, a surge nucleus may be a permanent, seasonal, sporadic or incipient feature of certain segments of the glacier. We think of a surge nucleus as a limited sliding region adjacent to or contained by stuck parts. Trapridge Glacier has a permanent nucleus in the form of a temperate region (sliding) up-stream from a sub-polar region (frozen stuck). Variegated Glacier may have a seasonal or sporadic nucleus. Seasonal and sporadic nuclei may occur because of annual cycles of melt-water flow. Incipient nuclei may be present in the sense that for some years conditions to form the nucleus may not be met because of the glacier state. Then, after some rebuilding, the basal shear stress is sufficiently high locally to cause more rapid sliding and a limited surge nucleus is formed. We have no stronger evidence for our suggestion of surge nuclei. But we would claim that the inhomogeneity of glaciers is such that seasonal, sporadic or incipient rapid-sliding regions are entirely possible. They may be very limited, but conditions for the initiation of surging would be those that determine whether the nucleus spreads rather than those that determine whether they exist in a given case. If nuclei are not present, our surge model is without validity. That situation would be remedied by considering more detailed aspects of what causes low-velocity sliding to accelerate at least locally.

Given the existence of a surge nucleus, we now consider the question of whether it will spread. Consider the nucleus to be the segment AB in Figure 2. This region is sliding rapidly and the longitudinal stress there is given by Equation (9) subject to boundary condition in Equation (13). Part of the weight of AB is no longer borne by the basal shear stress because the down-slope component exceeds the net drag. This excess has now been thrust on to the neighboring non-sliding parts and this will cause the basal shear-stress singularities as discussed previously. The discontinuity will move and the sliding region will spread if J at the discontinuity exceeds the critical level.

It is likely that a sliding region like AB would be quite short to begin with, i.e. less than one glacier thickness. The model would be characterized better if calculations of J for short sliding regions were available. Calculating J for the edges of sliding regions shorter than or comparable to h is quite difficult and would require numerical methods. Only when the sliding region has spread to somewhat longer than h can a leading order estimate be obtained easily. We will assume that such a surge nucleus exists.

The details of the calculation of J for the discontinuities at either end of a sliding region like AB are given in the Appendix. When AB is longer than h, the value of J at the discontinuity at B is

to leading order and there will be an equivalent value at A. This nucleus will spread if J _A or J _B equals or exceeds the critical value J _c and will be stable otherwise. As the region spreads, either σ _xx (x _A) or σ _xx (x _B), or both will increase in magnitude where A and B always define the limits of the current sliding region. If the critical value J _c is fairly uniform in the glacier, the lower extremity B will move down to the snout of the glacier at x = l. The upper limit A may travel to an ice reservoir or be blocked by some basal feature at which J _c is high. A third possibility is that the large longitudinal stresses will break off the surging part. We consider the final position of the upper limit to be at x = 0.

The rate of extension of the sliding region would probably be controled by the time necessary to incubate damage, e.g. tunnel closing around the boundary-condition discontinuity. However, it is clear from Equation (9) that the lower discontinuity at B would move down the glacier as a compression front as observed in Variegated Glacier by Reference PatersonKamb and others (1985). The upper discontinuity would move up the glacier as a tensile front. Of course, if the surge nucleus is originally near the top of the glacier, then the moving tensile front would not exist for long. However, a tensile region is explicit in our treatment of the model. If the ice is incapable of sustaining the tension that develops, then the surging part would break off. The resulting static imbalance would lead to unsteady motion of the glacier until surging terminates. We have not analyzed this possibility, although it may be important. Presumably, the drag forces involved would still be sufficiently high to keep the velocity of the surging parts to relatively low levels.

All or parts of the glacier are now sliding rapidly and stretching of the ice mass occurs. As discussed previously, we introduce into our model a surge termination in a phenomenological way. The surge may terminate due to loss of lubricating water or because of seasonal changes. The termination may occur first near the top and move down the glacier, as seems to be the case in Variegated Glacier or it may occur simultaneously throughout the glacier. The glacier returns to its quiescent state and the surge nucleus is once more stable. This will be so if the relevant surge nucleus is in a region of the glacier which has thinned during surging – generally, in the upper part of the surging segment. This can be seen from Equation (18) where the thickness enters directly. In addition, the total weight of the surge nucleus would be less than before and so the longitudinal stress σ _xx at the discontinuity would be lower than before. The surge nucleus would be stable until the height returns approximately to pre-surging levels due to rebuilding of the glacier. Note also that the inverse viscosity parameter B enters Equation (18) making it possible that seasonal fluctuations of ice rheology play a role in determining the value of J at a discontinuity.

Evolution of the surging glacier

The shape change during surging should be calculated while the surge fronts move down the glacier and until surging terminates. However, for simplicity, we will calculate the evolution by assuming that no shape change occurs while the surge fronts are moving through the glacier. Once a full surge part is established, the shape change of stretching and compression, and thinning and thickening occurs. This modeling is deficient, especially for Variegated Glacier where a contained surge region moves down through the glacier and surging terminates shortly after the surge region nears the snout. Our assumptions are equivalent to saying that the surge fronts move more rapidly than motions controled entirely by ice rheology which may be the case in certain glaciers.

At this stage in our model a lengthy segment of the glacier is sliding rapidly unconstrained at the snout. The complete solution during this period requires the solution of Equation (8) for the stress s σ _xx along with the determination of the velocity field and the glacier shape from Equation (11) and the equation for mass conservation. For ease of calculation, it appears to be most convenient to cast the equation for mass conservation in terms of the glacier volume between a postion x and the snout l,

In terms of V, conservation of mass takes the form

where we have neglected the generally small effect of accumulation and ablation during the surge period. Lastly, noting that ∂V/∂x = –h, we can write Equation (20) as

For simplicity, we will assume that the bed of the glacier has a uniform slope, i.e. y ¹ _b vanishes, and then the four equations for the four unknowns V, h, σ xx, and u which fully determine the motion and stress field for the surging glacier are

For illustration purposes, consider the case when the basal shear stress is such that it supports a fraction k of the down-slope component of the glacier weight, i.e. τ_s = kρghsin α. In this case. Equation (24) gives

and therefore

In addition, by dividing Equation (22) by V and then by differentiating it with respect to x, we can write Equation (22) as

Now, combining Equations (23), (27), and (28), we obtain the following equation for the stress

which is a non-linear kinematic wave equation. Using the initial particle position X and time as independent variables in place of the spatial variable x and time, this equation is generally written as an evolution equation along the characteristic curves, i.e.

on curves given by:

In this case, the characteristic curves determined from Equation (31) are simply the trajectories of the material particles x = x(X, t) where X is the initial position of a material particle. This is equivalent to using a Lagrangian formulation as opposed to an Eulerian formulation.

The solution to Equation (30) for the stress σ _xx of a material particle in terms of the initial stress field σ₀ is

where σ₀ is given in terms of the initial glacier shape by

and l ₀ is the initial snout position. The stress at a fixed point x in space is found from Equation (32) using the relation between x and X obtained by integrating Equation (31)

Here u is determined from the stress–strain-rate relation in Equation (23),

To conveniently use Equation (35), however, we must first change variables from x and t to X and t as follows

so that we have from Equation (35)

This can now be integrated to determine u(X, t) which in turn can be substituted into Equation (34) to determine the position x of a particle which started at a location X at t = 0. Using Equation (36), we find for Equation (34)

By interchanging the order of integration, this can be simplified to

where we have used Equation (33) and we define

In Equation (37) T is a dimension less time given by

Lastly, from Equations (25), (26), and (32) we find the ice thickness to be

A comparison between theory and observation was made using the data given by Reference Dolgushin and OsipovaDolgushin and Osipova, 1975 for the surge of Medvezhy Glacier which occurred in 1963. From the figures presented by Dolgushin and Osipova, data for the shape of the surging glacier were extracted using a Hewlett-Packard 7470A plotter and digitizing sight. The figures were enlarged to improve the resolution of the data. Using the initial configuration of Medvezhy Glacier in 1963, immediately prior to the surge, we computed the evolution of the glacier shape from our model. Results for three equal time intervals are shown in Figure 5. In addition, in Figure 6 we compared the observed post-surge shape to the shape calculated in our model when the distance the snout had advanced was coincident with that of the actual glacier. The agreement is quite remarkable. Some of the discrepancies can be attributed to the fact that we have assumed a uniform bed slope which is not quite accurate for Medvezhy Glacier and also to the fact that we have used a very simple form for the basal shear stress τ_s.

Fig. 5. Calculated surging motion of Medvezhy Glacier starting with the observed initial configuration. Glacier profiles are shown for four values of T: 0, — ; 0.58, .... ; 1.16, — — ; 1.74, − − − ; the glacier-bed axis and ice thickness are non-dimensionalized by the initial surge length l₀ and initial thickness at x = 0, respectively.

Fig. 6. Comparison of the final glacier profile between theory (− – −) and observation (....) for the surge of Medvezhy Glacier. The initial configuration is also shown (—). The glacier-bed axis and ice thickness are non-dimensionalized by the initial surge length l₀ and initial thickness at x = 0, respectively.

One feature to note concerning the present solution is the singularity in the stress field which occurs at a finite time (see Equation (32)). This is a fairly common characteristic for non-linear equations of the type considered here, i.e. Equation (29). Furthermore, note that the velocity field is determined from the stress field, and, consequently, it is also singular. This can be seen in Figure 7 where we have plotted the position of the snout of Medvezhy Glacier versus time. The singularity in the stress field occurs at T = 1/(H³) _max. and corresponds to the glacier thickness disappearing at some point as can be seen from Equation (38). The singularity in the velocity field actually precedes this by a small amount and is found to occur at T = (1 –e⁻³)/ (H³) max’. The singularity occurs partly because of the lack of inertial effects in our model. The physical significance of the approach in time to this singularity might be interpreted as follows. Since there is no limiting stress level or velocity, the velocity and longitudinal stress of a surging glacier will steadily increase and will eventually reach a point in time when the basal water system which is responsible for maintaining the high sliding velocities will be disrupted by the large strain-rates present. This will result in a drag greater than the down-slope component of the weight and the basal velocity will return to a level typical of the non-surge period when the glacier is essentially stuck to the bed. A new surge will not commence until the glacier, which has been thinned by the surge, has had sufficient time to build up once again and re-develop its basal water network. Ultimately, the down-slope component of the weight will again exceed critical levels in surge nuclei and a new surge will be initiated.

Fig. 7. The calculated position of the snout as a function of time for Medvezhy Glacier. The position of the snout is non-dimensionalized by the initial snout position l₀ and

is a dimensionless time.

From the present analysis, we can estimate the ratio of the recovery time (i.e. the time it takes the glacier to be rebuilt) to the surge duration time. We assume that the recovery time is the time it takes for the inflow at x = 0 plus the accumulation between x = 0 and l to equal the volume depleted during the surge. The volume depleted by the surge is the volume of ice which passed the original snout position l ₀. From our calculations for Medvezhy Glacier, we find this to be equal to Q _s h ₀ l ₀ where h ₀ is the maximum ice thickness and Q _s = 0.1911. For simplicity we assume that the net accumulation or ablation between x = 0 and l ₀ during the recovery period is equal to Q _A h ₀ l ₀ where Q _A is a constant. In the lower part of the glacier where the surge took place we expect spring and summer melting to exceed winter snowfall, so that the material carried beyond x = l ₀ by the surge would be melted away rather quickly. In any case, the mass balance defining the recovery period is

Consequently, the time to replace the volume, i.e. the recovery time, is

where an overbar on h ⁵ represents a time average in the time interval t = 0 to Δt _r.

On the other hand, the surge duration can be estimated from our calculation for Medvezhy Glacier. We found that the position of the glacier snout coincided with the observed position when the time T = 1.74, i.e.

where Δt _s is the surge duration time. Therefore, the ratio Δt _r/Δt _s is found to be

where γ = [(h /h ₀)⁵] is the time average over the time interval t = 0 to Δt _r Although it is difficult to use Equation (41) directly, it reveals the relevant physical parameters which control the surge cycle. Furthermore, we can use Equation (41) to estimate the parameter k which is a measure of the basal shear stress. For Medvezhy Glacier, l ₀ = 6000 m, h ₀ = 200 m, and the observed duration times are: Δt _r = 12 years; Δt _S = 2 months. In addition, from our calculations we found Qs = 0.1911 and, if we further assume that accumulation is responsible for replacing half of the surged volume, i.e. Q _A = 1/2Q _S, and we take γ; = 1, then we have the following estimate for k.

Recalling the definition of k, this indicates that during the surge of Medvezhy Glacier 88% of the down-slope component of the weight was being carried by the basal shear stress, the remainder being supported by the longitudinal stress.

Examining Equation (41) more generally, we see that as k → 1 (i.e. the basal shear stress supporting all of the weight), Δt _r /Δt _s → 0 or →t _s ∞, and there is no noticeable surge since the motion is occurring so slowly. Furthermore, Δt _r/Δt _S is a maximum and the surge is the most dramatic (since Δt _s is a minimum) when we have perfect slip, i.e K→ 0.