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 Reference NyeNye (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)

(Reference NyeNye, 1952). Here ρ is ice density, g the gravitational acceleration, and A and n are parameters in a power-type ice flow law (Reference GlenGlen, 1955; Reference NyeNye, 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. Reference NyeNye (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. Reference BuddBudd (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, Reference BindschadlerBindschadler (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 Reference NyeNye (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 Reference BindschadlerBindschadler (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 (Reference McCracken and DornMcCracken 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 Reference Budd and JenssenBudd 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 Reference BindschadlerBindschadler, 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 Reference BindschadlerBindschadler (unpublished).

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 Reference NyeNye (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 Reference NyeNye (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 Reference NyeNye (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 (Reference BindschadlerBindschadler, 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; Reference NyeNye, 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 Reference NyeNye (1960).

Background

Variegated Glacier has a well documented surge history with a period of roughly 20 years (Reference PostPost, 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 (Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 (Reference Bindschadler, Bindschadler., Harrison, Raymond and GantetBindschadler and others, 1976) and sliding is thought to contribute less than 5% of the total annual motion during 1973 (Reference Bindschadler, Bindschadler, Raymond and HarrisonBindschadler and others, 1978). the model can be applied.

Parameterization

The longitudinal line of markers established in 1973 (Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 Reference NyeNye (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⁻¹, Reference GlenGlen. 1955, and n = 3, A = 0.0817 bar⁻ⁿ year⁻¹, Reference NyeNye, 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 Reference GlenGlen’s (1955) values of A and n appears in Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 Reference NyeNye (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 Reference NyeNye (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 Reference NyeNye (1965). The second was to estimate Q by using the seismicaliy measured cross-section and the measured surface velocity averaged over the width, which Reference NyeNye (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 Reference BindschadlerBindschadler, 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 Reference NyeNye (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 Reference BindschadlerBindschadler, unpublished, fig. 2.11).

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 Reference BindschadlerBindschadler, 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 (Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 Reference NyeNye’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 (Reference Robin and WeertmanRobin 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 , −|å|.