Introduction
One of the primary goals of glacier dynamics is to predict the response of glaciers to climate variations. Reference NyeNye (1960, Reference Nye1963[a], Reference Nye[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 Reference Campbell and RasmussenCampbell and Rasmussen (1969), Reference Rasmussen and CampbellRasmussen and Campbell (1973), and Reference MahaffyMahaffy (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. Reference BuddBudd (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 Reference NyeNye (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 (Reference BuddBudd, 1968). Because no theory of basal sliding yet proposed (Reference LliboutryLliboutry, 1968; Reference WeertmanWeertman, 1964; Reference KambKamb, 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 (Reference Budd and JenssenBudd 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 Reference NyeNye (1963[b]).
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, Yi and ice depths, Hi at coordinates Xi. Changes in ice depth are determined from continuity equation using calculated volume fluxes Qi±½, mass balance ai, and channel shape (not shown).
A Crank‒Nicholson finite-difference formulation was used to approximate Equation (1). It has the following form:
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
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
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
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
and the surface velocity is
(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
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
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 n 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:
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
where
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 2, Δx 2). More details concerning stability and errors are given by Reference BindschadlerBindschadler (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 m1/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−n year−1 and n = 4.2) were assumed. Both f and f * were set to 0.55; values typical of valley glaciers (Reference NyeNye, 1965; Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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−1. 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.
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
where k represents a combination of constants. If the surface slope is constant, the kinematic wave velocity is
By rewriting Equation (12),
and, recalling Equations (3) and (4),
the kinematic wave velocity is
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−1 without diffusing. The surface velocity of the ice is 48 m year−1; 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−1.
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−1.
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−1 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).
Application of model to Variegated Glacier
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−n year−1, Reference GlenGlen. 1955, and n = 3, A = 0.0817 bar−n year−1, 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, f1: ____________, 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−n year−1, n = 4.2 and A = 0.087 bar−n year−1, 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 × 105 m3 year−1 km−1) 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 m2 year−1 in Equation (2) was chosen for the calculations. This corresponds to a 0.01 m year−1 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).
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 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 × 108 m3. 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−1.
Base shear stress
From each depth profile the corresponding profile of base shear stress can be calculated from
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−1). 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−1 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 × 109 m3. 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
The units of F are meters year−1; 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 Reference Robin and WeertmanRobin 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 (Reference KambKamb, 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 Reference Robin and WeertmanRobin and Weertman (1973, Equation (13)), the condition for this to occur is
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 Reference KambKamb (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 Reference KambKamb (1970, table 2), 0.007, and the corresponding value of G, 1.5 (Reference KambKamb, 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 Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 Reference BuddBudd (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−2. This location agrees in general with the position of the trigger zone suspected by Reference Bindschadler, Bindschadler, Harrison, Raymond and CrossonBindschadler 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 (Reference Bindschadler, Bindschadler, Raymond and HarrisonBindschadler 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.
Acknowledgements
This work comprises the major portion of the author’s doctoral dissertation, written while attending the Geophysics Program, University of Washington, Seattle. Financial support was provided by NSF Grants No. DES 72-01629 A01 and EAR 76-22463 A01. Frequent interactions with Dr Charles Raymond contributed to this research and his many editorial suggestions improved the text. I also thank an anonymous reviewer for many helpful comments.
