A Numerical Model of Temperate Glacier Flow Applied to the Quiescent Phase of a Surge-Type Glacier

Abstract

A time-dependent numerical model of temperate glacier flow without sliding is developed and applied to the quiescent phase of surge-type Variegated Glacier, Alaska. The model is based on a one-dimensional continuity equation but the transverse channel shape is explicitly included allowing the complex geometries of real glaciers to be modelled. Velocities and volume fluxes are calculated from the glacier geometry. Transverse stress is taken into account by shape factors which are fitted to measurements of geometry and velocity and are chosen to be insensitive to changes in geometry. Longitudinal stress gradients are taken into account by use of a large-scale surface slope. A Crank-Nicholson finite-difference approximation is used and it is unconditionally stable when a small contribution from the local slope is added to the average slope.

Model parameters are fitted to extensive data collected on Variegated Glacier in 1973 and 1974. Predictions of the model over a four year interval agree well with field measurements. Predictions of the current quiescent phase (1965–84) indicate depth increases in the upper glacier of more than 75 m with a twenty-fold increase in the volume flux. During this interval the base shear stress increases 40% in the upper glacier and decreases 20% in the lower glacier. During the mid to late quiescent phase, ice motion becomes more important than mass balance in the redistribution of mass over the central region of the glacier. If normal flow were to persist, the predicted steady-state profile would be an average of 100 m deeper and 41% more voluminous than in 1973.

The predicted base shear-stress gradient is never negative enough to satisfy Robin and Weertman’s (1973) condition for blockage of subglacial water flow. The annual rate of water production by dissipation of mechanical straining at the bed remains two orders of magnitude below that produced by summer surface melt. The predicted fractional increase in base stress during the quiescent phase is a maximum in the region believed to be the trigger zone of the surges.

Un modèle numérique de l’écoulement d’un glacier tempèrè et son application à la phase de repos d’un glacier à crue. Un modèle numérique à variable temporelle de l’écoulement d’un glacier tempéré sans glissement a été développé et appliqué à la phase de repos du glacier à crue Variegated Glacier, Alaska. Le modèle est basé sur une équation de continuité à une dimension mais la forme des profils en travers du chenal est explicitement prise en compte ce qui permet de reproduire les géométries complexes des glaciers réels. Les vitesses et les débits sont calculés à partir de la forme du glacier. La contrainte transversale est prise en compte par des facteurs de forme qui sont adaptés aux mesures de géométrie et de vitesse et sont choisis pour être insensibles aux changements de géométrie. Les gradients de contrainte longitudinale sont pris en compte en utilisant une grande échelle de pente superficielle. On utilise une approximation aux différences finies de Crank-Nicholson et elle est toujours stable lorsqu’une petite variation de pente locale vient s’ajouter à la pente moyenne.

Les paramètres du modèle sont ajustés grâce aux abondantes données recueillies sur Variegated Glacier en 1973 et 1974. Les prévisions tirées du modèle pour un intervalle de quatre ans concordent bien avec les mesures sur le terrain. Les prévisions pour l’actuelle phase de repos ( 1965–84) indiquent des accroissements d’épaisseur de la partie haute du glacier de plus de 75 m avec une multiplication par vingt du débit. Pendant cette période la contrainte de cisaillement à la base augmente de 40% dans le haut du glacier et décroît de 20% dans le bas. Pendant le milieu et la fin de la phase de repos le mouvement de la glace devient plus important que le bilan pour la redistribution des masses dans la région centrale du glacier. Si l’écoulement normal devait persister, le profil prévu au stade d’équilibre serait en moyenne de 100 m plus épais et de 41% plus volumineux que celui de 1973.

Le gradient prévu pour la contrainte de cisaillement à la base n’est jamais assez négatif pour satisfaire à la condition de Robin et Weertman (1973) pour le blocage de l’écoulement d’eau sous-glaciaire. Le taux annuel de production d’eau par dissipation d’énergie mécanique le long du lit reste de deux ordres de grandeur en dessous de celui dû à la fusion estivale en surface. L’accroissement partiel prévu des efforts à la base est maximum dans la région que l’on croit être celle du déclenchement des crues subites

Zusammenfassung

Ein numerisches Modell für das Fliessen temperierter Gletscher und seine Anwendung auf die Ruhephase eines ausbrechenden Gletschers. Für den Fluss temperierter Gletscher ohne Gleiten wird ein zeitabhängiges numerisches Modell entwickelt und auf die Ruhephase des ausbrechenden Variegated Glacier, Alaska, angewandt. Das Modell beruht auf einer eindimensionalen Kontinuitätsgleichung, enthält jedoch explizit das Querprofil des Gletscherbettes und erlaubt damit die Berücksichtigung der komplizierten Geometrie wirklicher Gletscher. Geschwindigkeiten und Massenfluss werden aus der Gletschergeometrie hergeleitet. Querspannungen werden durch Formparameter berücksichtigt, die sich aus Messungen der Geometrie und Geschwindigkeit ergeben und so gewählt werden, dass sie gegen Änderungen der Geometrie unempfindlich sind. Gradienten der Längsspannungen werden aus der Oberflächenneigung im grossen hergeleitet. Es wird eine Crank-Nicholson-Näherung mit finiten Differenzen angewandt, die ohne Einschränkung stabil ist. wenn kleine, lokale Neigungsänderungern zur mittleren Neigung hinzugefügt werden.

Die Modellparameter wurden den extensiven Daten angepasst, die am Variegated Glacier 1973 und 1974 gewonnen wurden. Die Vorhersagen des Modells über einen Zeitraum von vier Jahren stimmen gut mit den Feldbeobachtungen überein. Vorhersagen für die laufende Ruhephase (1965–84) lassen Dickenzunahmen von mehr als 75 m im Oberteil des Gletschers erwarten, verbunden mit einer 20-fachen Zunahme des Volumenflusses. In dieser Phase wächst die Scherspannung am Untergrund im oberen Gletscher um 40% an und nimmt im unteren Gletscher um 20% ab. Um die Mitte und gegen das Ende der Ruhephase ist für die Umverteilung der Eismassen im Zentralbereich des Gletschers die Eisbewegung von grösserer Bedeutung als die Massenbilanz. Würde die Normalbewegung anhalten, so wäre das vorausberechnete stationäre Profil im Mittel 100 m dicker und 41% voluminöser als 1973.

Der berechnete Gradient der Scherspannung am Untergrund ist niemals so stark negativ, dass er die Bedingung von Robin und Weertman (1973) für die Unterbrechung des subglazialen Wasserflusses erfüllen würde. Die Jahresrate der Wasserproduktion durch Umwandlung mechanischer Spannungen am Gletscherbett bleibt zwei Grössenordnungen unter der Produktion von Oberflächenschmelzwasser im Sommer. Die vorausberechnete anteilige Zunahme der Spannung am Untergrund während der Ruhephase hat ihr Maximum in jenem Gebiet, das für die Auslösezone der Ausbrüche gehalten wird.

Introduction

One of the primary goals of glacier dynamics is to predict the response of glaciers to climate variations. Nye (1960, 1963[a], [b]) was the first to attempt this in a systematic, quantitative way. He made a perturbation analysis of a one-dimensional equation obtained by integrating the continuity equation over the transverse cross-section area. His analytical studies examined the response of a glacier to small departures from the steady state. Subsequent models have usually relied on computers to solve a more complex set of equations. Models developed by Campbell and Rasmussen (1969), Rasmussen and Campbell (1973), and Mahaffy (1976), solved a two-dimensional equation obtained by vertical integration of the continuity equation. This technique is useful in modelling two-dimensional surface features and ice sheets. Budd (1975) has developed a one-dimensional numerical model than can undergo cyclic variations simulating surge behavior through the effect of a basal shear-stress boundary condition sensitive to lubrication by basal melting.

The model developed in this paper is designed to follow the details of a flow field more carefully than has been achieved in past models. In a strict mathematical sense, the model is one-dimensional but the transverse shape of the channel is parameterized by shape factors to include its important effect on response behavior. The values of these shape factors are determined from the measured motion and corresponding geometry of the glacier at one time. The development of the model is such that these shape factors are insensitive to changes in the glacier surface and represent a considerable improvement over the parameterization used by Nye (1960) when the changes in geometry are large.

Longitudinal stress gradients are approximately accounted for by using a surface slope averaged over a distance many times the glacier depth in the calculation of the base shear stress (Budd, 1968). Because no theory of basal sliding yet proposed (Lliboutry, 1968; Weertman, 1964; Kamb, 1970) has been able accurately to predict sliding velocities from a known glacier geometry, sliding is ignored in this model.

To insure numerical stability of the finite difference equations, a small local slope contribution must be added to the larger-scale average slope. The resulting set of equations is unconditionally stable; no restrictions are placed on either the time increment or the grid interval. No smoothing of the elevation profile is ever required (Budd and Jenssen, [1975]).

The model is applied to the quiescent phase of the surge-type Variegated Glacier, Alaska. The surge phase of this glacier cannot be simulated with this model because sliding is neglected. The large magnitude of the geometry changes that have been measured during part of the quiescent phase provide a favorable test of the model. In addition, the relative importance of mass balance and ice motion in these geometry changes can be examined quantitatively. By modelling the entire quiescent phase, the temporal changes in various geophysical parameters can be predicted and results compared with some proposed ideas on surge release mechanisms.

The Numerical Model

Continuity equation

The fundamental equation of the model is the continuity equation integrated over the transverse cross-sectional area to achieve the one-dimensional form given by Nye (1963[b]).

(1)

where S is the cross-sectional area, Q the volume flux, ȧ the mass balance (ice-equivalent units), and W the surface width. The grid system chosen is Eulerian, defined by points X _i fixed in space (see Fig. 1). The one-dimensional axis is curvilinear following the central flowline of the glacier.

Fig. 1. Coordinate system. Glacier geometry is specified by bed elevations, Y_i and ice depths, H_i at coordinates X_i. Changes in ice depth are determined from continuity equation using calculated volume fluxes Q_i±½, mass balance a_i, and channel shape (not shown).

A Crank‒Nicholson finite-difference formulation was used to approximate Equation (1). It has the following form:

(2)

where the subscripts refer to the grid point and the superscripts refer to the time step (time = m Δt). The volume flux (and velocity) are not calculated at the grid points x _i , but rather at the mid-points of segments between grid points. x _{i± 1/2} (Fig. 1). The grid spacing can be variable but should not vary rapidly. Equation (1) is one-dimensional if S, Q, ȧ, and W depend on known constants and a single variable. The independent variable is vertical ice depth, which, along with the fixed bedrock configuration, defines the glacier geometry.

Glacier geometry and mass balance

The surface width is expressed as

(3)

where D _i and E _i are constants to be specified and H _i is the vertical ice depth. Equation 3 represents the linear combination of channels with a parabolic and V-shaped cross-section respectively. The cross-sectional area is related to ice depth by integrating Equation (3) to give

(4)

The case F _i ≠ 0 corresponds to a grid point where Equations (3) and (4) define a channel shape over a limited range of H _i . Although S _i (H _i = 0) ≠ 0 in this case, the fit will be better within the intended range of depths than if F _i = 0. Specifying the central bedrock elevation Y _i at each grid point completes the detailed description of the glacier channel.

The mass balance may be specified as a function of position, elevation, and time as

(5)

Flux distribution

The principal problem is to relate the volume flux to the depth distribution. Assuming the glacier deforms only by simple shear and the surface slope α is small and nearly parallel to the bed, the base stress is

(6)

and the surface velocity is

(7)

(Nye, 1952). Here ρ is ice density, g the gravitational acceleration, and A and n are parameters in a power-type ice flow law (Glen, 1955; Nye, 1957), u _b is the sliding velocity, and H is the ice depth normal to the surface. In this model, u _b is assumed negligible. Drag along the valley walls also helps support the glacier. Nye (1965) took this effect into account by modifying the base shear stress in Equation 6 by a multiplicative factor f. This factor is referred to here as the “velocity shape factor” and has values between zero and unity.

This theory neglects the effects of longitudinal stress gradients on the ice motion. Budd (1968) derived a complex expression for the basal shear stress that explicitly included the effect of longitudinal stress gradients, but not all terms of the expression could be evaluated from the glacier geometry alone. Budd argued that this expression is simplified somewhat when averaged over a longitudinal distance of five to ten times the depth but the evaluation of certain terms still cannot be done rigorously. Using a still longer averaging length (10 to 20 times the depth) the expression reduces to the form of Equation (6). Using field data from Variegated Glacier, Bindschadler (unpublished) concluded that the expression

(6′)

where the averaging of slope was over a distance of 8 to 16 times the depth, provided calculated values of surface velocity (Equation (7)) which were as accurate as those it was possible to obtain with the more complex expressions.

The volume flux can be calculated from

(8)

where Q is the volume flux, S is the cross-section area, and f ^* is the ratio of velocity averaged over S to the center line surface velocity, referred to here as the “flux shape factor”. Based on numerical calculations of cross-section flow by Nye (1965), f and f ^* (and thus their product f ^* f ⁿ in Equation (8)) are insensitive to changes in geometry. Once f and f ^* are determined, they may be assumed to remain constant in time without introducing significant errors into the calculations.

For the numerical calculations, the velocity and flux (Equations (7) and (8)) are calculated from the following set of finite difference equations:

(9)

(10)

(11)

(12)

where the large-scale surface slope is averaged from x _k to x _l. The cosine of the local slope (Equation (9)) appears in Equations (11) and (12) to convert the vertical depths H _i to surface-normal depths but is strictly valid only when surface and bed are parallel. Velocity and flux have the same sign as α but not, necessarily, α. Thus local regions of up-glacier slope without up-glacier flow are allowed.

Boundary conditions

Two boundary conditions must be specified: one at the head, the other at the terminus. Three general types have been developed for the model and are illustrated in Figure 2. The first is a constant flux condition where the volume flux into the first section (or out of the last section) of the glacier is held constant. The second is a deforming wedge at either the head or the terminus. The conditions within this wedge are that volume is conserved while the volume flux vanishes when the ice depth vanishes. The slope of the wedge and the terminus position are determined by the glacier dynamics and vary in time. The final possible boundary is an ice divide at the glacier head. In this case a second ice body is considered to be flowing in the opposite direction from the glacier head. The volume flux feeding this second glacier is a specified fraction k of the volume flux feeding the modeled glacier. A symmetric ice divide corresponds to k = ‒1. The details of each type of boundary are discussed in Bindschadler (unpublished).

Fig. 2. Boundary conditions developed for model. See text for description of each case.

Method of solution

Equations (2) through (5) and (9) through (12) represent the complete set of equations that are solved on a computer. Due to the non-linearity of the equation set, the depth profile at the unknown time step cannot be solved exactly and a Newton–Raphson predictor-corrector scheme (McCracken and Dorn, 1964) is used to solve Equation (2) by iteration to the desired accuracy. The practical accuracy of the solution is constrained by the accuracy of the field measurements.

The numerical stability of this set of equations has been examined by both theoretical analysis assuming n = 1 and experimentally with n ≥ 1. When the velocity and flux (Equations (11) and (12)) depend only on the local slope, i.e , the system is unconditionally stable. If the large-scale slope in Equation (11) is replaced by a constant, i.e. non-diffusive flow, the system is marginally stable; errors neither grow nor decay. However, with the use of the large-scale slope in Equation (11), the system becomes unconditionally unstable.

To use this system of equations, then, either the bedrock and surface profiles must be considerably smoothed initially and smoothed frequently during the computations (see Budd and Jensen, [1975]) or some other means of restoring stability must be found. By replacing the large scale slope in Equation (11) by a weighted linear combination of the unstable large-scale slope and stable local slope (Equation (9)), stability can be restored and we can still maintain use of a large-scale slope to account for longitudinal stress gradients. Thus, Equation (11) becomes

(11′)

where

(13)

The quantity 〈α _{i + ½} 〉, is referred to here as the “effective slope”. Theoretically determined, the allowed values of ϕ depend slightly on the length of the large-scale average but an upper bound of 0.8 is characteristic (see Bindschadler, unpublished, figure 5.4). This bound was verified by experimental runs of the model. Although the use of this effective slope is primarily to restore numerical stability, it is also more realistic to use a combination of averaging scales rather than to rely on a single averaging scale. Use of a single averaging scale results in the peculiar behavior that any surface feature (a sawtooth wave, for example) with a periodicity equal to the averaging length would produce a constant surface slope and thus tend to be permanent and not diffuse with time.

The numerical stability of this revised system of equations is unconditional; any time step or grid interval, however large, is permissible. The practical limitation is of course that the finite difference method only approximates the differential equation (Equation (1)). The difference between the true and numerical solutions, called the truncation error, can be shown to be O(Δt ², Δx ²). More details concerning stability and errors are given by Bindschadler (unpublished).

Comparison of model predictions and analytical theory

To establish that the numerical calculations could accurately simulate glacier behavior, various test cases were considered. In all test cases a straight parabolic channel (D = 57.7 m^1/2, E = 0, F = 0) was used with a uniform bed slope of 5° and a uniform grid spacing of 200 m. Glen’s (1955) values of the flow parameters (A = 0.148 bar⁻ⁿ year⁻¹ and n = 4.2) were assumed. Both f and f ^* were set to 0.55; values typical of valley glaciers (Nye, 1965; Bindschadler and others, 1977).

Steady-state test glacier

In the first set of tests an altitude-dependent mass balance, a symmetric ice divide at the head, and a wedge at the terminus were specified (see Fig. 2). No large-scale surface slope was used (i.e. ϕ = 0 in Equation (13)). The bed slope was increased locally near the head (see Fig. 3). A glacier of initially arbitrary shape was allowed to adjust to the steady-state configuration shown in Figure 3. In this configuration, volume was shown to be conserved to within 1% over most of the glacier and errors in the velocity determination arising from the linear interpolations between grids (see Equations (11’) and (12)) were also 1% or less. With a different grid the same steady-state configuration was achieved, confirming that this configuration was determined by the mass balance and channel geometry rather than the grid.

Fig. 3. Steady-state configuration for test glacier model with diffusion and ϕ = 0. Upper plot is longitudinal profiles of centre-line bed and ice-surface elevations (vertical exaggeration, 5:1). Lower plot shows longitudinal profiles of depth, volume flux, and center-line surface velocity.

Perturbations from steady state

This steady-state glacier was then used to study the dynamic behavior of the model further. Two perturbations of the steady state were examined: addition of a uniform (1 m) slab, and a uniform mass balance increase of +0.1 m year⁻¹. The test glacier’s response in each experiment was compared to the analytical treatments of Nye (1963[b]). Figure 4 shows the results of each experiment. Although precise quantitative comparison with the analytical treatment was not possible (the numerical model presented more complex situations), qualitatively the two methods agreed very well.

Fig. 4. Temporal response of test glacier (Fig. 3) to: (a) addition of one meter slab at t = 0 years, and (b) uniform mass-balance increase of 0.1 m year⁻¹ (ice equivalent). Depths are shown as deviations from equilibrium depths (Fig. 3). EQL marks equilibrium-line position.

A more quantitative comparison was possible with the second set of tests where only a section of an ideal glacier was considered. The section was of uniform thickness (300 m), the mass balance was zero everywhere, and the volume flux into, out of, and at each point within the section was the same. Thus the section was in a steady state. No large-scale slope averaging was used (ϕ = 0). A small perturbation was added to the surface and its behavior was followed in time. Two cases were considered: non-diffusive flow where the effective surface slope is constant in time (the variations of flux depending only on the depth variations); and diffusive flow where the effect of slope variations is included. In both cases waves were expected.

The analytical prediction of wave velocity is easy to make based on the work of Nye (1965). For the parabolic case, Equation (12) becomes

(14)

where k represents a combination of constants. If the surface slope is constant, the kinematic wave velocity is

(15)

By rewriting Equation (12),

and, recalling Equations (3) and (4),

the kinematic wave velocity is

(16)

for the values used here (n = 4.2, f ^* = 0.55). By a different method Nye (1965) reached a similar result. To have ignored the effect of width in the above analysis would have led to a predicted kinematic wave velocity of (n+1)u _s = 5.2 u _s, more than 100% too large. This illustrates the importance of using a realistic two-dimensional channel cross-section.

The results of this set of numerical calculations are shown in Figures 5 and 6. The perturbation is a Gaussian wave one meter high (0.3% of the total depth) with a 1200 m half-width. Figure 5 presents the non-diffusive case and corresponds exactly to the above analytical analysis. Here the wave travels along the surface at a velocity of 117 m year⁻¹ without diffusing. The surface velocity of the ice is 48 m year⁻¹; the ratio is 2.46 as predicted above. The instability appearing behind the perturbation in Figure 5 is a result of the marginal stability of the equations when there is no slope variation in the determination of flux. Small positive errors in the depth form regions of higher velocity which become progressively deeper and faster. This unstable behavior was always observed whenever there was no diffusion used.

Fig. 5. Motion of Gaussian hump on uniform slab without diffusion. Depths are expressed as deviations from equilibrium depth (300 m). Successive profiles are vertically offset. Dashed line shows velocity of hump maximum is 117 ± 2 m year⁻¹.

Fig. 6. Motion of Gaussian hump on uniform slab with diffusion. Depths are expressed as deviations from equilibrium depth (300 m). Successive profiles are vertically offset. Dashed line shows velocity of hump maximum is 120 ± 2 m year⁻¹.

When diffusion was introduced, i.e. non-constant surface slope, Figure 6 shows the marked difference in behavior of the wave. The velocity of the wave crest is nearly unchanged (120 m year⁻¹ in this case) but the perturbation has diffused rapidly. In addition, the rate of diffusion agrees with diffusion theory. Considered separately from the wave velocity, the initial perturbation should have the following profile:

at a later time (Bindschadler, unpublished, p. 132). Thus after twenty years, the wave should be 0.26 m high with a half-width of 4630 m: from Figure 6 a height of 0.26 m and a half-width of 4600 m were measured. An important property of glacier behavior becomes strikingly clear when comparing Figures 5 and 6: although diffusion does not markedly alter the kinematic wave velocity, it is responsible for transmitting the margins of the perturbation down-glacier (and possibly even up-glacier) at speeds much higher than the kinematic wave velocity. Thus for the dynamic experiment shown in Figure 4a, the duration of the unstable response of the lower glacier is greatly curtailed because the leading edge of the wave (whose arrival restores stability; Nye, 1960) travels down-glacier at a speed much higher than the kinematic wave velocity. This also has the effect of dampening the intensity of the whipcrack-like behavior of the terminus described by Nye (1960).

Application of model to Variegated Glacier

Background

Variegated Glacier has a well documented surge history with a period of roughly 20 years (Post, 1969). The last surge occurred some time between August 1964 and August 1965 so the next surge is expected in about 1984. Through the quiescent phase this glacier undergoes large changes in its geometry making it a favorable test case for the numerical model. In addition, extensive data collection over several years (Bindschadler and others, 1977) provides the opportunity to parametrize the model accurately and to compare the model predictions with field measurements. Since this glacier is temperate (Bindschadler and others, 1976) and sliding is thought to contribute less than 5% of the total annual motion during 1973 (Bindschadler and others, 1978). the model can be applied.

Parameterization

The longitudinal line of markers established in 1973 (Bindschadler and others, 1977, fig. 1) provided a convenient center-line axis and grid spacing (Δx = 250m). In the accumulation basin some additional grid points were added where stakes were absent. In all, 76 grid points were used over a distance of 18.8 km. The 1973-74 seismic reflection survey provided ice depths every 500 m from Section F (at 5.1 km) to 1.2 km below Section B (Fig. 7). Detailed transverse profiles were completed at the six sections denoted F, EF, E, D, C, and B in Figure 7 (see Bindschadler and others, 1977, fig. 3). These data, along with maps of the exposed valley walls, allowed the parameters Y _i (bedrock elevation), D _i , E _i , and F _i (Equations (3) and (4)) to be determined. Figure 8 summarizes the mass-balance data used in the model.

Fig. 7. Variegated Glacier. Vertical ice depth measured from surface at center-line and longitudinal profiles of surface elevation (June 1973) and nominal bed elevation at center-line. Dashed portions indicate where data were interpolated.

Fig. 8. Profiles of width-averaged mass balance of Variegated Glacier from 1972 through 1976. Dashed portions indicate where data were interpolated.

Velocity shape factors

The values of the velocity shape factors f _i were estimated by two methods. The first used the numerical calculations of Nye (1965) for rectilinear ice flow through a parabolic section assuming no sliding and no longitudinal stress gradients. His tabulated values of f versus the half-width to depth ratio for n = 3 were interpolated to estimate f at each midpoint of the grid from measured values of width and depth. The second method used Equation (11′) to match the geometry to the velocity. The September 1973 elevation profile was used along with the annual velocity from September 1973 to September 1974 (estimated in the upper glacier and the terminal lobe). The averaging distance was two kilometers and ϕ = 0.8 for the effective slope (Equation (13)). Figure 9 compares the profiles of f _i determined by the two methods. Two profiles from the latter method are included for two different choices of the flow law parameters (n = 4.2, A = 0.148 bar⁻ⁿ year⁻¹, Glen. 1955, and n = 3, A = 0.0817 bar⁻ⁿ year⁻¹, Nye, 1953). A value of f greater than unity is unrealistic and indicates that either the velocity or the depth was poorly estimated above 3.1 km and below 16 km. A similar comparison of these two methods between Sections B and F using only the two kilometer average of the surface slope (ϕ = 1) and Glen’s (1955) values of A and n appears in Bindschadler and others (1977, fig. 7). These two figures show that the high-frequency variation of f values from the second method in Figure 9 is caused by the inclusion of the more rapidly varying local surface slope in Equation (13). Both analyses show similar discrepancies on the larger scale (especially from Section C to Section E) between the f values estimated from Nye (1965) and Equation (11′). This difference is attributed to the presence of longitudinal stress gradients which are ignored in the first method but which are accounted for in the second, although not perfectly, by the longitudinal averaging of the surface slope. Additional factors which might have contributed to the discrepancies are the non-parabolic shape of the channel and non-zero sliding rates.

Fig. 9. Velocity shape factor, f₁: ____________, from Nye (1965) assuming a parabolic cross section and n = 3; __________, required by measured surface velocity (Equation (11’)) assuming two flow laws (A = 0.148 bar⁻ⁿ year⁻¹, n = 4.2 and A = 0.087 bar⁻ⁿ year⁻¹, n = 3).

Accurate values of f are important because of the non-linearity of Equation (11). To support the assumption used in this model that f is time independent, Bindschandler (unpublished, fig. 7.8) calculated another profile of f values based on the measured geometry and motion of Variegated Glacier in 1976 and showed that the variation from the 1973 values was small.

Flux shape factors

The final parameter to be estimated is the profile of flux shape factor f ^*. Three different methods were used. The first was to use the values of f ^* versus the half-width to depth ratio (and interpolations between them) calculated by Nye (1965). The second was to estimate Q by using the seismicaliy measured cross-section and the measured surface velocity averaged over the width, which Nye (1965) showed was nearly equal to the velocity averaged over the section area, and to apply Equation (8) to calculate f ^*. This method was possible only at the six lettered sections. The third method used field measurements of elevation change and mass balance (over the 1973–74 balance year) and Equation (1) to obtain a longitudinal profile of ∂Q/∂x which was then integrated along the distance of the glacier. Because neither the elevation change nor the mass balance was well known at either end of the glacier, a large offset error in the calculated volume flux was expected. The longitudinal variation of Q over the central half of the glacier fits the estimates of Q at the individual sections (second method above) exceptionally well, especially after a slight rotation (−5 × 10⁵ m³ year⁻¹ km⁻¹) of the Q profile. The physical meaning of such a rotation was that a systematic error of either W(ΔH/Δt) or ȧW existed.

Figure 10 presents the estimates of f ^* by all of these methods. The first method gives nearly constant values for the central part of the glacier and is less than the individual section values estimated by the second method except for Section E. The differences between these two methods may be a result of non-zero longitudinal stress gradients and non-parabolic cross sections. The third method gives a much larger variation of f ^* which agrees well with values for the individual section. The large discrepancy at Section B is not a major concern because the surface velocity and flux are so low in this region. The discontinuity at 6.1 km is a result of the entrance of a tributary to the mainstream and is discussed below. A final observation is that in the region between 7.6 and 11.6 km, the trend of decreasing f ^* with distance corresponds to a region of positive longitudinal stress gradients and when the longitudinal stress gradients turn negative, the trend in f ^* reverses also (See Bindschadler, unpublished, fig. 3.1). Whether this is coincidence or not is unknown. A freehand f ^* profile following that produced by the third method but matching the values at the six sections was the profile used in the predictions of the next section.

Fig. 10. Flux shape factor, f _i ^*:-------, from Nye (1965) assuming a parabolic cross section; □, from width averaged surface velocity at measured cross sections, ---, inferred from elevation changes and net balance (Equations (2) and (8));______, idealized estimate, see text.

Fig. 11. Comparison of elevation change from June 1973 datum to September 1976 for measured data (•), fitted polynomial (----), and model prediction (________).

The entrance of a majority tributary alters the otherwise one-dimensional character of the glacier causing a discontinuity in the one-dimensional flux profile. Lacking sufficient information on the flow of the tributary, a simple parametrization of this tributary was sought. Because the accumulation basins of the mainstream and tributary are adjacent and cover a similar range of altitudes and because the tributary also participates in the surge, the tributary volume flux was simply scaled to the mainstream volume flux just above the confluence. Based on 1973 data the tributary flux was estimated to be 44% of the mainstream flux just up-stream of the confluence. Initial runs of the model showed that this value did not produce any large irregularities in the predicted depth changes in the vicinity of the confluence. This ratio was assumed constant for all model runs.

Boundary conditions

A constant flux at the head was specified. This is justified by the fact that the upper glacier does not participate in the surge. Thus, considering the large depth changes which occur over the rest of the glacier, this upper region can be assumed to be in steady state. The lower portion of Variegated Glacier consists of a lobe of stagnant, debris-covered ice. Both the velocity and mass balance reach negligible values well up-glacier from the terminus and the specific boundary condition used does not affect the solution over any significant part of the glacier. A maximum allowable residual of 10 m² year⁻¹ in Equation (2) was chosen for the calculations. This corresponds to a 0.01 m year⁻¹ uncertainty in the mass balance over a one kilometer width. A time step of 0.1 year was used.

Test of model prediction and data

To test the accuracy of this parameterization, the model was run for three years with an initial configuration corresponding to the measured September 1973 elevation profile. The resulting predictions of depth change were compared with those actually measured in September 1976 (Figure 11). The agreement is quite good. On the lower glacier the disagreement may be due to a poorly parameterized mass balance which overestimates the insulating effect of the extensive rock cover in this region. This is shallow, stagnant ice and has little effect on the behavior of the model. Near Section F the scatter of the data is too large to locate the position of maximum depth change. A more active tributary would raise the predicted maximum. The more general systematic error in Figure 11 is a reflection of a systematic error between field measurements of surface elevation change and mass balance (especially in 1974–75) (see Bindschadler, unpublished, fig. 2.11).

Predicted adjustment of Variegated Glacier

The model is now used to make predictions about possible future, and past, configurations of Variegated Glacier.

1973 to 1984

Of major concern is the configuration of this glacier just before the next surge. The model was run for eleven years, from 1973 to 1984, the best estimate of the time of the next surge. Figures 12, 13, and 14 show the predicted condition of the glacier at these times. Only one 1984 profile is shown in each figure, yet a number of models using the various profiles of f, f ^* and flow parameters discussed in the last section were run. The profiles shown in Figures 12, 13, and 14 used values of f fit to the measured surface velocity (Fig. 9) and the estimated f ^* values given by the heavy solid line in Figure 10, but the differences in the predicted depth changes using different profiles of f and f ^* were small (see Bindschadler, unpublished, fig. 8.1). In Figure 12, the reference surface is as measured in June 1973 while the initial condition for the model is the September 1973 surface. This demonstrates the large magnitude of depth changes expected over the last half of the quiescent phase relative to one summer’s ablation. From Figure 13 it can be seen that the dominant tendency is to smooth the flux profile with time. The glacier responds quickly in regions of large flux gradients (Equation (1)). Depending on the longitudinal variation of the channel cross-section and f ^* (see Equation (12)) this smoothing may or may not smooth the surface velocity (Fig. 14).

Fig. 12. Predicted depth deviation profiles from June 1973 datum for the times −8(1965),0(1973), + 11(1984), +21, +51, and +150 years since 1973.

Fig. 13. Predicted volume flux profiles for depth profiles shown in Figure 12.

Fig. 14. Predicted center-line velocity profiles for depth profiles shown in Figure 12.

1965 to 1973

Although the model cannot be run backwards in time (the diffusion term in Equation (11) would produce instabilities), the glacier profile in 1965 just after the last surge can still be estimated. By beginning with a 1965 profile the model predictions for 1973 can be compared with the actual measurements. By adjusting the 1965 estimate after each run to compensate for differences between the 1973 predictions and measurements and maintain smoothness, an estimate of the broad-scale features of the 1965 profile can be obtained (Fig. 12).

Over the region where the bedrock elevation is known, this 1965 estimate produced a maximum error of –5 m in the 1973 prediction with most points within ± 2 m of the measured 1973 profile. Figures 13 and 14 include the corresponding profiles of volume flux and surface velocity. Figure 14 shows that the ice in the reservoir area becomes nearly stagnant immediately following a surge.

Figures 12 through 14 illustrate the predicted changes of Variegated Glacier during a full quiescent phase. A maximum depth increase of 80 m is predicted in the reservoir area and a depth decrease of 50 m in the receiving area. The region of zero depth change is near Section E, not far from the mass-balance equilibrium line. An earlier estimate of this boundary for the 1964–65 surge based on aerial photographs placed it much lower on the glacier; somewhere below Section C (Bindschadler and others, 1977).

If the surge expected in 1984 transformed the 1984 profile to the 1965 estimate, the average elevation gain below Section E would be about 30 m. The approximate gain in volume in this region would be 2.5 × 10⁸ m³. If this volume flowed through section E in a one-year surge, the flux would be a fifty-fold increase from the 1973 value. Since the actual surge probably occurs in less than one year, surge velocities probably reach 100 times the 1973 value or about 360 m year⁻¹.

Base shear stress

From each depth profile the corresponding profile of base shear stress can be calculated from

(17)

where the slope–depth product is averaged over two kilometers to partially account for longitudinal stress gradients. Using values of f from Nye’s (1965) numerical calculations (which assumed no basal slip, no longitudinal stress gradients, and a parabolic bed) shown in Figure 9, values of τ_b were calculated at each center-line position where the depth was known (about every 500 m). This profile is smoother than the one calculated from the velocity profile (which also includes the local slope variations). The choice is not crucial since the following analysts will concentrate on changes of τ_b with time. The smoothed profiles of τ_b corresponding to the predicted geometries for a number of dates are shown in Figure 15. During the surge cycle a large, almost consant increase in τ_b occurs above Section E. Below 12 km τ_b decreases by, again, an almost constant amount. In between the temporal variation is more complex.

Fig. 15. Smoothed base shear stress (Equation (17), two kilometer average) for depth profiles shown in Figure 12.

Steady state

This model cannot simulate the surge itself due to the absence of any sliding contribution to the velocity. During the quiescent phase the glacier attempts to reach a steady-state configuration but this process is interrupted by a surge. However, the model can simulate this approach to steady state assuming this unstable behavior were not to occur. The specific steady-state configuration depends on the mass balance. In this model mass balance is a function of distance and not altitude, so the computed steady state may be slightly under-estimated. Using the same parametrization as before, the model was run for 150 years by which time the depths were changing by a negligible amount (less than ±0.01 m year⁻¹). The profiles of depth change, volume flux, surface velocity, and base shear stress for 21, 51, and 150 years after 1973 are included in Figures 12, 13, 14, and 15. The model predicts that the steady-state configuration would be much thicker than that at any time during the quiescent phase. The maximum depth would be 520 m (Section E). The steady-state volume flux (Fig. 13) equals the balance flux everywhere by definition. For a periodically surging glacier, this equality must hold when integrated over a complete surge cycle. Thus the balance flux must be higher than the volume flux throughout the quiescent phase but less than the fluxes realized during the surge. The predicted steady-state velocity profile bears little resemblance to the velocity profile at any time during the quiescent phase (Fig. 14). The double peak (at 4.8 and 6.8 km) is still present but the maximum velocity of 310 m year⁻¹ would occur much further down-glacier. The predicted base shear stress (Fig. 15) reaches a maximum of nearly two bars; maintaining the same general spatial variation as in 1973 but about 0.7 bar higher.

The volume of the glacier in 1973 was calculated as 3.43 × 10⁹ m³. The 1965 estimate was 1% lower. In 1984 the volume had increased less than 3%. These volume changes are negligible considering the accuracy limits of the model (about 1%) and possible errors in the assumed mass balance. The steady-state volume, on the other hand, is 41% greater than in 1973. This increase in volume begins slowly because the integral of mass balance over the glacier is close to zero and the ice is relatively stagnant everywhere. Eventually, the ice in the accumulation region becomes thick enough to transport an increasing volume of ice to the lower glacier. This transition from stagnant to active ice is examined in more detail below.

How a non-negligible sliding rate affects these predictions cannot be examined quantitatively. Additional depth increases (decreases) would occur when the longitudinal gradient of sliding velocity has a high negative (positive) value. The indication from field measurements is that, in the later stages of the quiescent phase, sliding is becoming more important, yet because the inferred longitudinal variation of sliding is similar to that of the velocity due to internal deformation, no large deviations in the general behavior of the predictions made here are expected.

Relative importance of flow and mass balance

The temporal changes in depth that have been predicted for Variegated Glacier are a result of two processes: mass balance and ice motion. For a surge-type glacier the relative importance of each process varies through the surge cycle. During the surge, ice motion is most important, presumably through the action of fast sliding, while immediately following a surge, mass balance has been assumed to dominate (Robin and Weertman, 1973). To examine the relative importance of mass balance and motion quantitatively, a parameter F is denned as

(18)

The units of F are meters year⁻¹; the lower bound is −|ȧ| and is independent of time. A positive (negative) F signifies a region where flow of the ice is more (less) important than mass balance in determining the change of ice depth.

Figure 16 presents profiles of F for five different times between 1965 and steady state (F ≡ 0). The profile of −|å| is included to illustrate the size of |(1/W) (ΔQ/Δx)| everywhere. The succession of profiles shows the formation of a single F-wave that grows and spreads out as it travels down the glacier. Within this “wave”, flow is the more important factor in depth changes. In 1965 mass balance dominates. The two positive peaks in 1965 are due to the sharp boundaries of the reservoir area (Fig. 13) and are probably fictitious due to the sensitivity of Q to errors in either depth or surface slope. By 1973, flow has become more important than mass balance only in the vicinity of Section EF. In the final stages of the quiescent phase the region of flow-dominated depth adjustments has grown considerably, extending from Section EF to Section D, and within it F has increased in magnitude. Assuming no surge were to occur and that sliding remain negligible, the growth and movement of this F-wave continues: Behind the wave the relative importance of flow decreases and eventually steady state is achieved. There is no evidence that more than one wave ever forms. 150 years after 1973, F is near zero everywhere except near the terminus.

Fig. 16. Profiles of F parameter (Equation (18)) for depth profiles shown in Figure 12: -------,−8(1965);---, 0(1973);_______, +11;(1984);_________, +21; , +150 years since 1973, and , −|å|.

Discussion of surge release mechanisms

Quantitative measures of the base shear stress, and its variation in time and space, permit a limited number of surge release theories to be tested. One theory which can be tested is that by Robin and Weertman (1973). Briefly, this theory uses the idea that if the glacier is supported by normal stresses on an uneven bed, the lowest normal pressures occur on the down-glacier side of bumps in the bedrock and it is in these low-pressure regions that the subglacial water collects, draining from one low-pressure pocket to another. According to sliding theories, the magnitude of the pressure deviation across a bump is directly proportional to the base shear stress and inversely proportional to the bed roughness (Kamb, 1970, equation (123)). Thus, in a region where a negative base-stress gradient exists (i.e. base stress decreases down-glacier), the magnitude of the pressure deviation around each bump will also decrease down-glacier. This will cause a contribution to the pressure gradient between water pockets which decreases up-glacier, opposing the normal down-glacier flow of water. If large enough, this gradient could balance the gravitational forces thus blocking the drainage of water. Water accumulating up-glacier would eventually lubricate the bed enough to initiate a surge. According to Robin and Weertman (1973, Equation (13)), the condition for this to occur is

(19)

where ρ _i and ρ _w are ice and water density, g is gravitational acceleration, α and β are surface and bed slope respectively, and ξ and G are parameters used by Kamb (1970) to describe the bed roughness and effect on τ of the local variations in effective viscosity associated with the bumps. To obtain a lower bound for the right-hand side of Equation (19) the lowest value of ξ given by Kamb (1970, table 2), 0.007, and the corresponding value of G, 1.5 (Kamb, 1970, fig. 4), were used. Values for this expression were calculated every 500 m from the measured 1973 and predicted 1984 geometries along with the corresponding base-stress gradients (see Figure 17). The inequality required by Robin and Weertman is never satisfied anywhere on the glacier. More significantly, the changes from 1973 to 1984 are small for both terms of Equation (19). This strongly suggests that at no time during the later quiescent phase does this type of “pressure dam” occur.

Fig. 17. Longitudinal profiles of calculated negative base-stress gradient in 1973 and 1984 compared to −∆τ/∆x required to block flow of basal water (Equation 19).

Robin and Weertman hypothesized that mass balance dominated the changing geometry during the quiescent phase and thus τ would increase (decrease) in the reservoir (receiving) area with a zone of growing negative base-stress gradient in between, which they believe would be the trigger zone. The results of the model predict that flow plays a major role in the redistribution of mass in the mid-to-late quiescent phase (Fig. 16). This prevents the growth of a large base-stress gradient. Further, there seems to be no reason to expect the development of the largest negative base-stress gradient where the change in base stress is zero as do Robin and Weertman. This region of zero change in base-stress gradient (around 11 to 12 km) is much lower than the trigger zone suspected by Bindschadler and others (1977), near section EF. Nevertheless, it is conceivable that the negative base-stress effect may be a very local effect. The resolution of Figure 17 is not particularly good due to the local variability of depth (and f). It remains possible that locally large gradients exist but arc beyond the resolution of either Figure 17 or the model.

Another model of surge release has been proposed by Budd (1975). He hypothesizes that water produced by the straining of the basal ice layers and sliding lowers the base stress locally and allows an increase in the sliding rate. The efficiency of this lubricating process is hypothesized to depend on a single parameter, the “lubrication factor,” and if it is large enough both sliding rate and water production can increase in a positive feedback loop to produce surges. In this theory, the critical regions of a glacier are those where τ _b u (rate of energy dissipation) is high. Since base stress varies much less than velocity, the critical areas correspond to regions of high velocity. The maximum velocity predicted in 1984 is in the broad region between Sections EF and F where the maximum energy dissipated is 0.6W m⁻². This location agrees in general with the position of the trigger zone suspected by Bindschadler and others (1977). This theory, however, ignores the possible additional lubrication effects of water entering the glacier from the upper surface by melt water and rainfall. During summer this additional water far exceeds the volume produced internally and at the bed and is believed to be responsible for the increased velocity of the glacier during the summer months (Bindschadler and others. 1978).

Additional support for expecting the surge to begin near the vicinity of Section EF comes from examination of the fractional increase in base stress over the glacier during the complete quiescent phase. This fractional change was calculated using the 1973 values of base stress as the reference, i.e.

where t was either 1965 or 1984. Figure 18 shows that just down-glacier of Section EF the base stress increased dramatically during the quiescent phase. The non-linearity of the flow law (Equation (6)) amplifies the corresponding increase in deformation. Without a more reliable sliding law the effect of this increase upon sliding can only be assumed, but it is likely to be large as well.

Fig. 18. Fractional change in base shear stress (Fig. 15) from 1973 to 1965 and from 1973 to 1984.

Summary

The model presented here has enabled a number of important insights into glacier behavior. Studies of test glaciers confirm the sensitive response behavior of the terminus to small perturbations of the glacier (Fig. 4) and illustrate the modifications to the motion of kinematic waves resulting from diffusion (which quickly spreads the wave out; see Figs 5 and 6) and the valley walls (which slow the wave down: see Equation (16)). The application of the model to Variegated Glacier demonstrates how one can use actual field measurements to determine the model parameters accurately (Figs 9 and 10). The fact that this parameterization is insensitive to change in the glacier geometry emphasizes the advantage of this particular model. The large changes that occur allow a careful study of how a real glacier far out of balance with its climate responds (Fig. 16) and provide predictions which can be compared with field measurements (Figs 12, 13, and 14). Finally, the surge-type character of Variegated Glacier permits the testing of some theories of surge release using the predictions of the model (Fig. 17) as well as the calculation of various dynamic parameters, such as base shear stress (Fig. 15), and the temporal changes of these parameters (Fig. 18) which may be important in determining surge behavior.