I. Introduction

It has been shown by Reference Campbell and RasmussenCampbell and Rasmussen (1969,1970) that a surge can be induced in a numerical model of a glacier by lowering the basal stress for some time. Physical mechanisms for the lowering of the basal stress of a glacier by the production of melt water have been discussed by Reference LliboutryLliboutry (1964, 1969) and Reference WeertmanWeertman (1962, 1969).

The question arises as to how a water lubrication mechanism for glacier sliding can be incorporated naturally into a numerical model. It is apparent from the high speeds and low basal stresses of the surging and the fast polar glaciers that sliding is the key to the problem.

A numerical model for glaciers in which sliding is small has been developed by Budd and Jenssen (in press) and found to give reasonable approximations to real glaciers over a wide range of sizes. To extend this model to surging glaciers, a realistic mathematical model for sliding must be developed.

A number of theoretical studies of sliding have been presented e.g. Reference WeertmanWeertman (1957, 1964, 1967), Reference LliboutryLliboutry (1958[a], [b], 1968[a], [b]), Reference NyeNye (1969, 1970), Reference KambKamb (1970). These theories have aimed at determining a direct relation between the basal stress and sliding velocity at a certain location of a glacier in terms of the roughness of the bed. Weertman and Lliboutry also deal with multi-valued relations between basal stress and velocity depending on the basal melt water. Reference ShumskiyShumskiy (1965) argued that a direct relation between the basal stress and the velocity at a point does not in general exist because the glacier as a whole governs what happens at each point. Reference MeierMeier (1968) and Reference HodgeHodge (1972) have shown that the Nisqually Glacier exhibits considerable independence between its basal stress and sliding velocity.

The model to be presented here introduces a simple theory of sliding which requires no direct relation between the basal stress and the sliding velocity at a point, since both are determined from the properties of the whole glacier.

2. Principles for a surging model

The non-sliding model of Budd and Jenssen (in press) is stable in that if the input, consisting of bedrock, balance, and the properties of the ice, are kept constant then the glacier eventually tends to a steady state from any initial configuration.

For a surging model the following principles are adopted:

• (a) With constant input consisting of the bedrock configuration, the accumulation-ablation balance distribution, and the properties of ice, a steady-state glacier does not result, but instead the glacier builds up slowly then surges rapidly, stagnates then builds up, and repeats the process periodically.

• (b) The same laws and ice properties apply to the first order to all temperate glaciers. The final state of a glacier, as surging or non-surging, would then depend on the glacier’s bed profile and balance profile.

• (c) Surging does not require as a necessary criterion any special features other than sufficient accumulation for the given bedrock profile. Other properties such as special bed smoothnesses, high geothermal flux, etc., are considered as secondary contributors to surging conditions.

• (d) Surging is a large-scale phenomenon which can be adequately represented at least to the first order by a numerical model using discrete values at spacings comparable to the ice thickness. It follows that small-scale bed features are not relevant except in so far as they contribute to bulk averages over larger scales.

The aim of this paper is to outline the theory for the simplest model of glacier flow which includes the class of periodically surging glaciers as a subgroup obeying the same laws.

Numerical results showing the surging sequences that result from such a model are given by Budd and McInnes (in press). These are discussed together with another example of a high-speed surging ice mass in Section 8. The use of the model for studying particular real surging glaciers will be presented subsequently.

3. System of equations

In order to highlight the essential features required for periodically surging, the model here will be kept as simple as possible, concentrating on reproducing the most important properties of real surges first. Greater generality and sophistication can be added later to give closer representation to reality as required.

Hence a two-dimensional model is chosen representing a vertical section along the central flow line of a glacier or any flow line of an ice cap. Parameterization for making the twodimensional profile appropriate for a given valley glacier is discussed by Reference Budd and JenssenBudd and Jenssen (in press). Additional complications such as special effects of tributaries or valley sides are omitted at this stage. Since very wide ice masses surge, such as Brúarjökull in Iceland Reference Thorarinsson(Thorarinsson, 1969), it may be expected that although valley sides are important to a particular glacier, they are not an essential feature of surges generally. The most important features of surging glaciers can be represented in two dimensions as changes of ice thickness and velocity with time along the length of the flow line.

Let Z be the ice thickness, V the average velocity through a column, and A the accumulation–ablation balance at a point distance x along the glacier. Cf. Figure (1).

Fig. 1.

A representation of a two-dimemional model Jot a glacier is shown, in part (i) a curve is given representing the net accumulation ablation A as a function of distance x along the glacier. In part (ii) b represents the given bedrock profile, Z the ice thickness, –α the surface slope, V the average total velocity. V _i the average internal deformational velocity, V _b, the sliding velocity and L the glacier length at time t = 2.

The variations of the ice thickness with lime t are then given by the equation of continuity

The balance curve A is given as input data together with the bedrock elevation profile b as a function of x along the glacier. Thus if the average velocity V can be expressed as a function of z and its derivatives at a point then this can be substituted in Equation (1) and the problem is reduced to one of kinematic waves. For sliding however no such relation between the velocity and the properties of the glacier at a point can be realistically assumed.

We separate the total average velocity V of a vertical column into a basal sliding component V _b and an average internal deformation component V _i as

For large-scale averages, if the longitudinal strain-rates are small and if a power law for flow applies, the average internal deformation velocity can be satisfactorily expressed in the form

where τ _b is the base stress and k and n are constants.

More general flow laws could be used here but it is shown by Reference Budd and JenssenBudd and Jenssen (in press) that this formula with n = 2 and k = 0.15 bar ⁻² year ⁻¹ gives a reasonable approximation to a wide range of glaciers of different sizes. For the case of large longitudinal strain-rates the octahedral stress could also be used. This offers no special difficulty numerically but represents an additional complication at this stage which is not a necessary requirement for the surging process. In fact since the base stress drops during the surge when the longitudinal stresses become high, a compensating simplification is also introduced by using for the average internal deformation the expression

where τ _c is the down-slope stress on the centre line given by

with ρ the ice density, g the gravitational acceleration, α the surface slope, and s a shape factor appropriate for the glacier cross-section, generally to be taken between 0.5 and 1.

Since

and s, ρ, g are constants, the internal velocity given by Equation (4) is determined at each point along the glacier by properties of the glacier at that position. For the basal velocity however we turn to the equation for longitudinal stress equilibrium for scales large compared to the ice thickness and for small surface slopes from Reference BuddBudd (1970):

where

is the average longitudinal stress deviator through the column at position x.

We define the average generalized viscosity η through the column

where

is the average longitudinal strain-rate through the column. For the simplest model η can be taken as constant. More generally a power or hyperbolic sine flow law could be used. In addition the octahedral stress could be used incorporating the horizontal shear stress as well as the longitudinal stress. However during surges, if the longitudinal strain-rates tend to dominate, it will be sufficient at this stage to neglect the effect of the horizontal shear on the longitudinal strain-rate.

It is emphasized that the flow properties of the ice subjected to longitudinal stresses during a surge are regarded as unknowns here to be determined from the deformation rates occurring during real surges. In fact we wish to study the behaviour of the model to a wide range of constant viscosity values in order to understand fully the role that the ice flow properties play in the mechanics of the surging before moving to a highly non-linear flow relation.

A further simplification is introduced by considering that the changes in the mean longitudinal strain-rate over the large scales we are dealing with are primarily associated with the changes in sliding velocity rather than the internal deformation, i.e.

Thus Equation (7) can then be written simply as

At this stage the basal stress τ _b has still to be determined. Two principles are used to derive it, viz:

• (a) The local lowering of basal stress is due to frictional lubrication.

• (b) The gross static equilibrium of an ice mass is preserved, i.e. complete break-off and avalanching are not considered here.

If the basal stress is lowered at some point of the glacier, then, by Equation (7). longitudinal stresses are introduced up-stream and down-stream. If the glacier is not to accelerate as a whole down the slope, then the basal stress must increase at some other part to prevent this. To preserve the gross equilibrium we require the total basal stress over the whole length L of the glacier to be given by

In order to take account of local lowering of basal stress due to frictional lubrication, we note that the rate of melt-water production is equal to the rate of energy dissipation due to the motion τ _c V. At this stage we do not wish to differentiate between water produced directly at the base due to the sliding (τ _b V _b) and that produced from the total energy dissipation τ _e V. This no doubt depends on the internal transport of water in the ice mass, which still requires further investigation. For the numerical model calculations it does not matter greatly which expression is used so we choose the total dissipation τ _c V for simplicity.

Thus we seek a formula which gives a local reduction in τ _b with increasing τ _c V. Many such formulae of various complexity can be developed, but at this stage we look for the simplest. Hence we define the local lubrication-lowered stress as

where φ is an adjustable parameter called the friction lubrication factor. Here we take φ as a constant, but for increased complexity it could be made variable and also include more complex functions of τ _c V.

To preserve gross equilibrium we define average stresses over the glacier by

The basal stress τ _b is then taken

This equation has the property that gross equilibrium is preserved, i.e.

and that locally the stress difference (τ _c – τ _b) is greatest where τ _c V is greatest.

Having determined the basal stress distribution, the basal velocity can be obtained by integration of Equation (10), provided we can specify the boundary conditions.

Since the boundary conditions are not well known the following will be taken as working assumptions.

(i) The longitudinal strain-rate due to sliding tends to zero at the terminus as Z tends to zero, i.e.

This says that as the ice thickness becomes small near the front, any sliding that occurs tends to be like a block at the end point.

(ii) The horizontal velocity at the start is taken as zero, i.e.

For ice masses starting at ice divides this is certainly true, but for glaciers starting on steep slopes some modification may be necessary.

Hence we obtain for the longitudinal strain-rate due to sliding

The sliding velocity is then given by

Thus, having determined the sliding velocity, the total velocity is given by Equations (2) and (4). This is then used in the equation for basal stress and the system can be cycled for convergence. For a time-varying numerical model, the velocity is fed into the continuity equation (1) and the system is cycled by time-stepping. The time step must be chosen sufficiently small so that the glacier dynamics can keep up with the changes in ice thickness, or otherwise more frequent cycling for the velocities is required, cf. Table I.

Table. I.

Calculation routine

4. Operation of the Model

To begin with, the model can be run with the lubrication factor φ zero, then the sliding is eliminated and the glacier grows to steady state. For steady state

and the glacier reaches a length L determined by

The maximum flux Φ_m occurs at the firn line, say x = f, and is given by

For sufficiently small non-zero φ, τ _b is slightly less than τ _c, so a small amount of sliding occurs, depending on the value of the parameter η. Provided η is not too small, a new steady state is reached similar to that for zero φ except that, for the position of maximum flux, V is slightly greater and Z is slightly less than for the case of zero φ. As φ is increased the sliding velocity increases and, provided η is sufficiently small, a critical stage is reached for which φ τ _c V is large enough to lower the base stress appreciably which causes the sliding velocity to increase and through the feed-back of Equations (5^∗) to (9^∗) gives rise to a high velocity zone near the region of the glacier where τ _c V is a maximum. This fast zone with longitudinal extension up-stream and compression down-stream moves as a dynamic (non-kinematic) flux wave and spreads out down the glacier. The surge finishes as the fast velocity zone spreads through the terminus carrying mass with it and causing the rapid advance. After the surge, the reduced slope causes a reduction in the stress and the velocity which through Equations (3^∗) to 9^*) reduces the sliding so that the motion reverts to the former slow mode. The advanced front then gradually becomes reduced by ablation while the accumulation zone thickens until τ _c V is sufficiently large once again and the surge is repeated.

The speed of the surge depends on the values of η and φ, but also on the length of the glacier, the thickness, the bedrock slope, etc. If η is too large, the sliding is also reduced below the critical value to surge.

5. Critical values

In order to obtain an idea of the appropriate range of values to use for φ, we examine the function τ _c ^∗. For the case of zero sliding, if we substitute V from Equation (4) into (12) we obtain

This function of τ _c for constant φ, k and Z has a maximum when

The studies of ordinary and fast glaciers for numerical modelling by Reference Budd and JenssenBudd and Jenssen (in press) suggest that a maximum τ _c of about 2 bars is appropriate. Similarly with Z ≈ 500 m the maximum deformational velocity would be V ≈ 300 m/year. Thus to obtain a maximum of τ _c for τ _b ≈ 2 bars and V ≈ 300 m/year. we may expect a maximum for φ of

Since the sliding velocity increases the value of τ _c V, τ _c ^∗ drops even lower, but Equation (25) gives a starting point for the program and we can examine the range φ ≈ 10⁻⁵ – 10⁻².

The maximum value of τ _c V may be expected to develop near the firn line for which from Equations (5) and (22)

Since this is independent of ice thickness, it follows that, for this model, glaciers of various sizes may surge, provided the input,

is sufficiently large. Glaciers of steeper slopes would be expected to reach the critical stage earlier, i.e. for smaller thicknesses. To apply this in practice however, it should be noted that the width needs to be parameterized.

To simulate the fast polar glaciers it is merely necessary to increase the accumulation zone flux to sufficiently high values to maintain the glacier flow in the fast mode, i.e. high velocities with relatively low base stresses.

Thus with this model a natural classification scheme for glaciers presents itself as follows:

• (i) “Ordinary glaciers” are those whose flux for their bed profile is low enough for the glacier to remain in the “slow mode”, i.e. with the base stress increasing with the velocity.

• (ii) “Fast glaciers” are those whose flux is sufficiently high to maintain the glacier in steady state in the “fast mode”, i.e. with a high velocity and a lowered base stress.

• (iii) “Surging glaciers” are those whose flux is sufficiently high for the glacier to reach the fast mode before reaching steady state but not sufficient to maintain it in the fast mode. Thus these glaciers build up slowly until they reach the fast mode, drain out rapidly in a surge then slowly build up again—periodically oscillating from one mode to the other.

In times of changing climate it is possible for surging or non-surging glaciers to pass from one type to the other by sufficiently large changes in the balance.

As an estimate of the order of magnitude of the critical flux from Equation (26), the studies of real glaciers for numerical modelling by Reference Budd and JenssenBudd and Jenssen (in press) and Reference Budd and AllisonBudd and Allison (in press) suggest that a value of sαVZ ≈ 500 m² year⁻¹ may separate the surging from the ordinary glaciers.

6. Special characteristics of surging glaciers

In order to test the hypothesis of a critical τ _b V we note from Equation (26) that for the critical zone separating surging from non-surging glaciers the surface slope would be inversely proportional to the flux rate φ times the shape factor s. To calculate the shape factor we really need the ice thickness, which is not available for most surging glaciers. However if we select surging glaciers over a wide range of surface slopes and flux rates then the precise value of the shape factor, which varies generally between 0.5 and 1, is not so important.

Hence in Table II a list of glaciers is shown with estimates of the following quantities:

Mean surface slope

,

Area of accumulation zone S,

Average net accumulation rate above the firn line Ā,

Total mass flux at the firn line Q = ĀS,

Width of the glacier near the firn line W,

Flux rate Φ = V _c Z where V _c is the centre-line surface velocity and Z is the centre-line thickness near the firn line.

Table II.

Approximate surface slopes and flux rates for some surging glaciers

The flux rate Φ is estimated from the total mass flux Q by the relation

where W is the mean width of the section and V is the mean velocity over the section.

Two further shape factors are introduced following Reference Budd and JenssenBudd and Jenssen (in press)

The flux rate per unit width’at the centre line is then given by

For the above calculations the product s ₂ s ₃ has been taken as ≈0.5.

From the values of Φ and

in Table II it is apparent that the inverse proportionality of Equation (26) does tend to apply.

Using the empirical graph of glacier surface slopes versus velocity as discussed by Reference Budd and JenssenBudd and Jenssen (in press) and Reference Budd and AllisonBudd and Allison (in press) these flux rates and surface slopes for surging glaciers can be compared to corresponding values for the ordinary glaciers and the fast polar glaciers, cf. e.g. Carbonnell and Bauer (1968).

The heavy dashed line in Figure 2 represents a curve of approximately constant α Φ which, for constant shape factor, would also be a constant τ _b curve. This curve tends to separate the bulk of ordinary glaciers from the fast polar glaciers. The plotted points for the surging glaciers tend to fall on the fast side of this curve. Thus for steady state these surging glaciers would tend to have high τ _c V values if they were to flow in the ordinary mode.

Fig. 2.

From a compilation of measurements on ordinary glaciers, ike mean centre-line ice thicknesses Z are shown as a function of surface slope α and centre-line surface velocity V. From these curves the centre-line flux φ = VZ obtained as shown by the full curves. Individual points for some fast polar glaciers are shown by the circles in the high-velocity regan, The crosses are obtained from the fluxes and slopes of some surging glaciers. These tend to be generally between the values for the ordinary glaciers and the fast polar glaciers. The heavy broken line represents a constant α φ curve which suggests this may define a transition zone between the ordinary and fast glaciers.

In spite of the scarcity of data for surging glaciers this empirical analysis suggests that the frictional lubrication basis for a surging model may at least match reality in the first approximation.

7. Matching real surges

In the present model the magnitude of the surge depends primarily on the amount by which the flux, for the given slope and cross-section, is greater than the critical value. Thus the size of the surge is altered by changes in the climate through the balance curve.

Similarly the period between surges is altered by changes in the balance associated with climate changes. The length of the period depends on the time taken for the glacier to build up to the critical size.

The duration of a surge for the model is dependent on the speed of the surge and the length of the glacier. Thus the values of η and φ in the model control the duration of the surge through their control of the surge speed. The durations of real surges have considerable variation, e.g. Kuliah 12 km in 2 months (Reference DesioDesio, 1954), Brúarjökull 8 km in 4 months (Reference ThorarinssonThorarinsson, 1969), Muldrow 6.6 km in about 2 years (Reference PostPost, 1960), Walsh 11.5 km in about 6 years (Reference PostPost, 1966; Reference Meier and PostMeier and Post, 1969).

Although the duration of a surge in the model of a particular glacier can be controlled by the η and φ values, it may not be possible with the present model to match the wide range of durations given above with the same values of η and φ for all glaciers. Two obvious extensions can be made to the model which will effect the surge duration. First a non-linear longitudinal stress-strain-rate relation would simulate a wide range of effective viscosities for surging, in that those which develop large longitudinal stresses would have lower viscosity and surge faster than those with smaller longitudinal stresses. Secondly the effect of the valley sides, especially on glaciers wide relative to their depth, can be important to the shear-stress distribution around the boundary of the cross-section. A simple partitioning scheme for separating vertical and horizontal shear stresses as discussed by Reference Budd and RadokBudd and Radok (1971, p. 32) could be useful for controlling the duration of the surge. This problem will be treated further in the context of using the model in matching particular real surges.

8. Sample results

There are a large number of different types of surging ice masses which could be examined by the present model. Just two examples are presented here to show the scope of the model in representing two particularly interesting types of different surging ice mass. The first is an example of a typical large valley glacier (such as Walsh, Muldrow or Steele) which moves at the rate of several km/year over a number of years. The second is an example of a high-speed surging ice mass (such as Brúarjökull) which moves an order of magnitude faster at the rate of several tens of km/year (over 100 m/d) for a period of a few months.

The same model is used in each case with just different values of the input, viz. the bedrock profile and balance profile (or two-dimensional parameterized influx profile) characteristic of the particular ice mass. In each case the model is examined for a wide range of the parameters viscosity η and lubrication factor φ.

For a particular run with the given input, the parameters η and φ are set, and the ice is grown from zero thickness and run for a peried of from 700 to 1000 years, by which time the “final state” of the ice mass is apparent.

The input bedrock and balance (or influx) curves are shown for the two cases in Figures 3 and 5. In each case for small values of φ and reasonable η the ice masses grow up to a steady state, as shown, e.g. in Figure 3b. As the lubrication factor φ is increased, the maximum value of the sliding velocity increases. For any value of η, not too large, a certain stage is reached for which the ice mass surges and continues to do so periodically until a final uniformly repeating surging state is reached as shown in Figures 4 and 6.

Fig. 3.

The input data far the model consists of the net accumulation/ablation balance curve (a) as function of distante along the glacier, and the glacier bedrock profile (b, lower curve). For the non-surging case the glacier builds up to a limiting steady stale profile (b, upper curve). For the surging case the glacier develops periodically a steep profile (c, 1) which rapidly transforms to an extended flatter profile (c, 2) during the surge.

Fig. 4.

The growth of the model towards steady state for an ordinary glacier (r) and a surging type glacier (í) is illustrated by the variation with ime of the glacier’s length (a), the maximum velocity (b), and the position of the maximum velocity (c) The ordinary glacier gradually ends towards a constant steady state whereas the surging glacier develops a periodically oscillating state characterised by a long period of slow build-up with retreat followed by a rapid advance at high speed. In the course of a surge cycle the position of maximum velocity travels down the glacier as a flux wave with increasing speed from the accumulation zone to the region of the terminus.

Fig. 5.

Forati idealized smooth analogue to a flow line of the Brúarjökull, the influx curve shown in (a) represents an accumulation! ablation balance curve modified to take some account of the varying width. A corresponding smoothed bedrock profile is shown in (b) and (c). The resultant surface profile of the ice mass after about 700 years is shown for a retreat phase at 10 year intervals in (b) and for a rapid advance phase in (c) at intervals of o.1, 0.1, o.15 years from the start of the advance.

Fig. 6.

The Change in the length of the idealized Brúarjäkull model is shown in (a) for the glacier grown from zero thickness with influx and bedrock profiles given in Figure 5. After about 700 years the ice mass exhibits regular periodic surging. In (b) the magnitude of the maximum sliding velocity is shown on a logarithmic scale. The sliding velocity is very low for most of the period, gradually increasing to about 50 m/year after which the surge sets in and the velocity reaches 50 km/year for a few months before dropping back to nearly zero to repeat the cycle.

For a given bedrock and balance profile the critical value of the lubrication factor for surging depends on the value of the viscosity η. The distribution of sliding velocities over the domain of lubrication factor and viscosity is shown in Figure 7. Here the surge–non-surge boundary is shown by the broken line which lies generally between the 10² and 10³, m/year velocity isopleths. For low viscosities this boundary is very sharp but for high viscosities it broadens until the viscosity becomes too high for surging altogether. A variation of the input data for the balance or bedrock profiles results in a general shift in the sliding velocities with respect to the axes, but the general pattern is preserved. The data for Figure 7 come from several hundred different runs growing the ice mass from zero to a final state using the same input but varying the values of φ and η. The output is very large and similar diagrams can be constructed for other features of the ice mass such as maximum thickness, surge duration, surge period, length of surge advance, etc., to show how these vary over the domain of φ and η.

Fig. 7.

The variation in the maximum sliding speed of the idealized Brúarjökull model of Figure 5, over the domain of different values for the viscosity η and lubrication factor φ, shown in m/year. The dashed line shows the division between the steady-state, to the left, and the surging results, to the right. The maximum speed for the steady-slate results usually occurs early in the development so that final sliding speeds much greater than 100 m/year tend to develop surging.

In this way it is possible to select some features of a measured glacier surge which can be matched by the model using certain values of η and φ.

For temperate glaciers the value of η depends on the stress. Even at a stress of one bar, field and laboratory values seem to be scattered over an order of magnitude. Thus it is necessary to consider a wide range of η values. Typical values should lie between 10⁻¹ and 10 bar year. For cold ice masses the viscosity could be typically one to two orders of magnitude higher and Figure 7 shows that under similar conditions the cold ice mass would surge with lower velocities than the corresponding temperate ones.

For the valley glacier example of Figure 3, a smooth bedrock and balance profile was selected similar to the measured values on Lednik Fedchenko but somewhat lower as one would expect without tributaries. An increase in φ to 10⁻⁴ bar⁻¹ m⁻¹ year⁻¹ results in a surging glacier with maximum speed of over 4 km per year, advance duration 3 years, surge period 88 years, and advance distance of 10 km, as illustrated in Figure 4. These surge features are not unlike those measured, e.g. on Walsh or Steele Glaciers. The change in the thickness profile is also of a similar magnitude as shown in Figure 3.

The progress of the surge is characterized by the movement of a broad wave down the glacier as shown by the position of the maximum velocity shown in Figure 4c, 2. A similar wave can be seen as a travelling thickness perturbation for the second example in Figure 5. This wave starts building up in the accumulation zone right after the surge and moves down the glacier with increasing speed towards the terminus. The rapid advance occurs as the wave passes through the terminus.

With the same values of η and φ, the surging glacier reverts to a non-surging glacier by a reduction of the balance curve. Similarly the non-surging glacier can be changed to a surging one by an increase in the balance curve.

Variations of the bedrock slope also shift the pattern of sliding velocities and the surge-non-surge boundary over the domain of η and φ.

Finally, even with this first simple model it seems to be possible to match measured glacier surges in many respects by an appropriate selection of the values of η and φ.

For the very high-velocity example, a smooth bedrock profile similar to a flow line for Brúarjökull is used as shown in Figure 5. The influx curve is a modified balance curve chosen to include some effects of the varying width of the ice mass along the flow line, although no attempt to give an exact match has been made at this stage. The resultant sliding velocity distribution over the domain of φ and η is that shown in Figure 7. For values to the left of the dashed line a steady state is reached while those to the right result in regular periodic surging as shown by the velocity and length changes in Figure 6. The steady-state profiles of ice thickness tend to be generally thicker than the measured profile, whereas the surging profiles of ice thickness as shown in Figure 5 are somewhat thinner. This is a result of higher average velocities for surging ice masses. It is expected that closer matching can be achieved by using a more representative influx curve.

Nevertheless the example given has many features similar to those associated with the measured surges of Brúarjökull. If a realistic value of η of I bar year is chosen, then a value of φ = 1.7 × 10⁻⁴ gives a maximum surge speed of 51 km/year and a surge duration of 0.3 year. These results are illustrated in Figures 5 and 6. The observed features of the 1963 surge included an average speed of 5 m/h (44 km/year) over a duration of less than 3 months (1969). The advance derived from the model is 7 km, a little less than that measured. The change in thickness also seems compatible with the measured changes. The period of the surging is about 64 years compared to 73 years between 1890 and 1963. It is expected that the period increases slightly with the ice thickness here.

The changes in the thickness profile for 10 year intervals through a retreat phase is shown in Figure 5b. The dashed line shows the most advanced position just after the surge identical to that of Figure 5C, which shows the changes in the profiles during an advance for the times 0, 0.1, 0.2, 0.35 years from the commencement of the advance. These illustrate the rapid movement of the thickness wave to the terminus and beyond within just a few months.

The pattern of thickness changes with time also have a remarkable resemblance in general form to the few data available from measured glacier surges, including the bulged advancing front and the straightening retreating profile. Figure 6 shows the change in length and sliding velocity with time. The logarithmic scale on the velocities illustrates the dramatic changes during and after the surge. The gradual approach to regular periodic surging over 700 years is clearly apparent.

Finally it also seems possible on this simple model to match these high-speed surging ice masses by appropriate values of the parameters. The difference in parameters from one ice mass to another will depend on a number of factors such as the effect of valley sides, the general bed roughness, and the geothermal flux. For zones with unusually large geothermal flux rates (say q in mechanical units), it may be appropriate to include this as an extra parameter contributing to the melt production rate by writing

This would tend to give surges more readily in regions of high geothermal flux. Variations in the appropriate η values from one ice mass to another will be dependent on the stress level as well as the temperature. It is expected that the use of suitable non-linear flow laws and octahedral stresses will clarify this.

The output for the programme includes profiles of other features of the ice mass such as base stress, longitudinal stress, longitudinal strain-rate, τ _c, τ _c ^∗, V, τ _b V, VZ, etc. The changes of these features during surges will be left to a separate paper in which a study of the matching of a variety of ice masses such as those of Table II will be presented.

9. Conclusions

The model presented here has been simplified to include just those features of an ice mass required to represent the essential physical processes of surging. These include the variations of ice thickness, stress, strain-rate and velocity along the central flow-line as a function of time. The information required for input include the bedrock and balance profiles, the flow properties of the ice, and the frictional lubrication factor.

Thus, with suitable two-dimensional parameterization, real glaciers, both surging and non-surging, can be modelled directly. A measure of the accuracy of the model can be obtained by a comparison of the calculated velocities and ice thicknesses with those of the real glaciers.

The model has considerable scope for improvement by increased generality and sophistication without changing the basic principles upon which it depends. The obvious extensions include generalization to three dimensions, use of more complex flow relations and octahedral stresses, and the development of more complex relations for base-stress lubrication. Even so, the present model has considerable flexibility in being able to match real surges in magnitude, speed, duration, and period.

By using the model to match real glaciers, important conclusions can be drawn on otherwise difficult problems such as the magnitude and distribution of the longitudinal stresses, the magnitude of the base stress, the sliding velocity, and the propagation of dynamic waves in glaciers during surging and non-surging conditions. The general problem of glacier sliding can be approached from a phenomenonological viewpoint with this model. Surging glaciers need not then be regarded as presenting any new difficulties, but rather as throwing considerable light on the processes of glacier sliding. This may then be used to improve the model to provide closer agreement with reality.

Thus the model presented here may be regarded as the simplest model to include periodically surging glaciers as a natural sub-group of glaciers generally.

Acknowledgements

I am indebted to Mr Barry M. Mcinnes of the Meteorology Department of the University of Melbourne for obtaining the large number of computer runs required for preparing the diagrams. The advice of Dr D. Jenssen of the same Department on the model programming has also been much appreciated.