A. Measurements

The geometry of Variegated Glacier in June 1973 is displayed in Figure 1. Changes in thickness subsequent to June 1973 were measured by surveying surface markers spaced along the center line over most of the glacier length (Fig. 2). Repeated surveys of markers gave velocity averages over various time intervals (Fig. 3). This paper discusses primarily winter and summer seasonal intervals, as defined by the list of survey dates given in Table I. The extent of spatial coverage, spatial density of data points, and measurement time-scales is displayed, for the most part, by these figures and the table.

Less systematic measurements were made at other times, at off –center –line locations, across several transverse lines, and in the principal tributary (Fig. 1). These accessory data are available in the data report series (Bindschadler and others, unpublished [a]–[d]; Raymond and others, unpublished [a]‒[c]) and are not presented in detail here. We will however, have occasion to refer to them.

B. Treatment of the Data

Reference Center-Line and Datum Geometry

To establish distance along the glacier in a unique way, a reference center line (Fig. 1) was defined to follow closely the center-line markers established in June 1973. The kilometer distance cited in the text and figures refers to the distance from the glacier head along this line.

Datum geometries for the glacier and valley were established by assigning values of surface elevation, center-line bed elevation, cross-section area, surface width, and change of surface width with unit change in depth at each datum center-line point (Fig. 1). Surface-elevation values were based on elevations measured in June 1973 at M ¼ km intervals and on hand interpolation where reference points were not located precisely at survey points. Bed-elevation values were based on seismic reflection measurements, radio echo-soundings, and bore-hole depths as compiled from various sources by Kamb and Echelmeyer (1986, fig. 8) except between Km 0 and Km 3 where radio echo-soundings by MacQueen (unpublished) were used. Cross-section area, width, and its change with depth were measured in 1973 by seismic reflection at six transverse lines (Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others, 1977) and these were interpolated at intervening locations based on U.S.G.S. map Sheets. The longitudinal positions of the six transverse sections are labeled by capital letters in the figures displaying longitudinal profiles. Above Km 3, the surface and bed geometry are complex and are not necessarily well represented by center-line surface and bed elevations.

Fig. 1.

Geometry in June 1973. (a) Map shows center line with Km ticks indicating distance from head, (b) Longitudinal profile gives the center-line elevations of the surface and the bed. Bed elevation measured by seismic reflection (Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others, 1977) is represented by dashed line. Subsequent radio echo-sounding (MacQueen, unpublished) and bore-hole depths (personal communication from W.B. Kamb) have provided improved depth estimates on the upper glacier. The current “best bed” is shown by the solid line, (c) Shows longitudinal variation in width, change of width with height, and shape factor.

Table I

Data on seasonal variation time windows

C. History of Change in the Main Valley

Elevation Changes

A progressive trend of annual elevation change from one September to the next occurred over most of the glacier length (Fig. 4). Between 1973 and 1981, there was a progressive thickening of the upper glacier from about Km 2 to Km 10 with a maximum of about 60 m between Km 6 and Km 8, a corresponding thinning of more than 50 m near the terminus below Km 18, and a steepening of the intervening stretch of the glacier. The boundary between net thickening and thinning relative to 1973 is well-determined after 1975 and moved from Km 9.8 to Km 11.6 between 1975 and 1981, which corresponds to a down-glacier migration at an average speed of 0.3 km a⁻¹. The maximum elevation change appears as a bulge, which also showed evidence of a down-glacier propagation from 1974 to 1981 at a similar average speed.

The total elevation change after 1973 is unknown above Km 2, since measurements at fixed index locations were made only after 1976. However, from 1976 to 1981, elevation changes were smaller than the annual accumulation and unsystematic. Also, aerial photographs show that this zone was only mildly affected by the surges of 1964–65 and 1982–83. There is apparently no progressive evolution of this highest part of the glacier. This nearly steady-state condition may arise because this zone is separated from the main length of the glacier by a short complex ice fall between Km 2.0 and Km 2.7 (Fig. 1b).

Fig. 2.

Elevation change relative to datum profile of June 1973. Note change in vertical scale for different years. Data points are from surveyed locations and elevations relative to the optimally interpolated datum profile. Dashed curve is fifth-order polynomial fit to the data points. Solid curve is an optimally interpolated fit to the data points. (Data points and digitized curves are tabulated in Raymond and others (unpublished [d]).)

Fig. 3.

Horizontal velocity measured at the surface for various years (× — summer, late June or early July to early September; о — annual, September to September). Summer and annual velocity curves are determined from the corresponding data points using optimal interpolation. Winter u_w curves are found from the annual u_y and summer u_s curves by using the equation t_yU_y = t_s.u_s + (t_y ‒ t_s)u_w where t_y and t_s are the lengths of the annual and summer survey intervals. For annual and winter measurements, the year indicates the end of the measurement period. (Data points and digitized curves are tabulated in Raymond and others (unpublished [d]).)

D. Changes in the Tributary

Changes in the one major tributary (Fig. 1) are not as well-documented as in the main valley. Based on up to six markers, a general description is possible. Elevation increased roughly 40 m from 1973 to 1978. This change was about the same as the elevation change in the main valley between Km 5 and Km 6.5 where the tributary enters (Fig. 4). From 1978 to 1981, there was no consistent trend in elevation change, which contrasts with the continuing progressive increase in the main valley (Fig. 4). Similarly, the speed increased from 1973 to 1978, but there was no consistent change from 1978 to 1981. There was also a distinct seasonal variation in speed similar to that observed in the main valley. In 1978 and 1979, the seasonal difference in velocity was about 0.3 m d⁻¹ compared to a winter velocity of 0.5 m d⁻¹. The history of the seasonal variation is not well known.

It appears that, by 1978, the tributary reached an approximate steady state in which the velocity had increased enough to drain the mass-balance input and stop the progressive build-up of ice mass there.

E. Prediction of Changes

Reference BindschadlerBindschadler (1982,unpublished) developed a numerical model to predict changes of Variegated Glacier between its surges. The model calculated changes in thickness with time from the longitudinal variation of ice flux and local width-averaged mass balance. Ice velocity and flux were calculated from ice deformation with no basal sliding. By fitting the model to the measured geometry and velocity in 1973 and driving it with the average mass-balance distribution observed from 1972 to 1976, he predicted changes in thickness and velocity in future years (Reference BindschadlerBindschadler, 1982, figs 12 and 14).

Fig. 4.

Elevation change from September 1973 to following Septembers, based on the optimally interpolated elevation-change profiles of Figure 2. Above Km 2, elevation changes were small and unsystematic from 1976 to 1981 and are assumed to be small relative to 1973.

Over the length Km 5 to Km 16 that Bindschadler examined, the shape of the predicted longitudinal profiles of elevation change and velocity agree remarkably well with the observed changes. The advance with time of the boundary betwen thickening and thinning, the growth of the bulge peaking between Km 6 and Km 7, and the increase in velocity over the upper part of the glacier are all predicted well. (Compare Figure 4 to Reference BindschadlerBindschadler (1982, fig. 12) and Figure 5b to Reference BindschadlerBindschadler (1982, fig. 14).) However, the size of the changes in both velocity and elevation predicted over a time interval comparable to 1973–81 are much smaller than what actually occurred. Thus, there is a quantitative discrepancy that needs to be examined.

The most fundamental discrepancy is associated with the under-prediction of velocity changes.Reference Bindschadler Bindschadler (1982) estimated deformational flow velocity from basal shear stress and ice thickness with rheological and geometrical parameters calibrated from measurements in 1973. He used a simple longitudinal averaging scheme to account approximately for longitudinal stress gradients; although this scheme had imperfections that could have led to erroneous local variations in velocity (Reference Kamb and EchelmeyerKamb and Echelmeyer, 1986), it should have been adequate to represent the large-scale flow pattern correctly. For lack of a generally accepted sliding law, Bindschadler was forced to assume no sliding in his model, which could be the source of the discrepancy. Whatever the cause, the under-prediction of velocity changes from geometry changes represents a serious problem for any attempt to model ice flow. In the remainder of the paper we shall be concerned with this discrepancy in velocity identified by the model results. However, we shall discuss it by examining directly the relationship between observed velocity and observed geometry changes.

The failure of the model to predict accurate elevation changes arises in part from the under-prediction of velocity changes mentioned above. In addition, Bindschadler (1982, unpublished) noted that observed surface mass balance and changes in ice thickness over the interval from 1973 to 1976 appeared to be inconsistent with the ice redistribution by the flow velocities that were actually measured. Ice thickness increased more rapidly than could be explained, especially above Km 7. Our own analysis of subsequent years confirms this inconsistency. It may be associated with problems that inevitably arise in the pragmatic application of a one-dimensional continuity model to a specific three-dimensional system described by incomplete data. It is necessary to account for differences between center-line and width-averaged thickness changes and mass balance, for input from ice-clad valley walls and all tributaries, and for possible changes in depth-averaged density.Reference Bindschadler Bindschadler (1982, unpublished) discussed the limitations of the data in assessing the possible impact of these sources of error in mass-balance accounting. We shall not discusss these problems here, but instead will limit our attention to the broader dynamical question of how observed velocity is related to observed geometry.

A. Background

The total surface speed of a glacier ucomes from two sources expressed as

(1)

Here u _d is a deformation contribution arising from internal shearing between the surface and bed and u _b is a sliding contribution.

The deformational contribution u _bis commonly estimated from

(2)

where h is the surface-normal ice thickness τ_b is the shear stress at the bed, and K = 2A/(n + 1) and n are parameters determined by the flow law of ice described as . Equation (2) arises from integration of shear strain-rate over depth, assuming that surface-parallel shear stress varies linearly with depth from 0 at the surface to τ_b at the bed and that it is the dominant contribution to effective stress and the corresponding effective viscosity in the non-linear flow law. (See Reference PatersonPaterson, 1981.)

Fig. 5.

Surface horizontal velocity in summer (a) and winter (b). Curves are from Figure 3 for locations below Km 3. Above Km 3, curves are based on hand interpolation of data values.

This section examines the effectiveness of Equation (2) to describe the relationship between changes in geometry and velocity observed on Variegated Glacier, and whether appropriate parameters K(or equivalently A) and n can be identified. To do this, we assume u _b«u _d so that u _b∼u, and followReference Budd and Jenssen Budd and Jensen (1975) by comparing depth-averaged shear strain-rate u _d/h with τ_b at different locations and times, which should define a line of slope non a log-log plot.

It would also be possible to express deformational velocity directly in terms of depth and surface slope without reference to basal shear stress using the assumptions above and additional assumptions that will be necessary to evaluate τ_b (e.g.Reference Kamb and Echelmeyer Kamb and Echelmeyer, 1986). This path would allow a direct comparison of velocity changes to changes in depth and surface slope. Echelmeyer (unpublished) has shown that a response factor is needed to account for changes in ice cross-section shape as an ice surface rises or falls in a fixed valley cross-profile. We avoid the need for explicit reference to the response factor by using Equation (2) and accounting for thickness, slope, and shape changes in the evaluation of τ_b as explained in section B below.

Fig. 6.

Variation of surface horizontal velocity versus time at two longitudinal positions measured on a seasonal time-scale. Curves show the variation in velocity predicted from internal deformation with no sliding, using Equation (6) and the geometry and velocity for winter 1973–74 as a reference state with n = 3 (solid curve) and n = 4 (dashed curve where deviation from solid curve is perceptible). Differential speeds measured over short intervals in bore holes are represented by points (personal communication from H.F. Engelhardt and W.B. Kamb).

B. History of Base Stress

Shear Stress

The base shear stress beneath the center line is commonly estimated from the glacier geometry through the equation

(3(a))

where ρg = 0.088 bar m⁻¹ is the weight per unit volume corresponding to ice density ρ (900 kg m⁻³), h is the local ice thickness normal to the upper surface at the center line, α is the surface slope, and f is a cross-section shape factor. Here, τ_sis referred to as the local slope stress. We calculate it using a slope averaged over 500 m; this corresponds to the resolution limit of our measurement spacing and is approximately 1–2 ice depths.

To account for longitudinal stress gradients, we used a longitudinal averaging based on the longitudinal coupling theory of Kamb and Echelmeyer (1986, Equation (3a)). It is expressed in terms of τ_sas

(3(b))

where lis a longitudinal coupling length that sets the effective length scale of the averaging. A constant value of l= 0.66 km (4l = 2.66 km) was used between Km 2 and 15. Above and below these limits, l= 0.25 km (4l = 1.00 km) was used because of lower ice depth. These coupling lengths are consistent with values found appropriate by Reference Kamb and EchelmeyerKamb and Echelmeyer (1986), who used their theory to analyze the distribution of velocity on Variegated Glacier during the years 1977–78. The base-stress distribution given by Equations (3a) and (3b) follows closely base stress estimated from Equation (3a) alone using a slope averaged over 2 km Reference Bindschadler, Harrison, Raymond and Crosson(Bindschadler and others, 1977), but it is somewhat smoother and, theoretically, more accurate.

Values of shape factor f were obtained from the numerical calculations tabulated by Reference NyeNye (1965, table IV), assuming a parabolic cross-section shape, width ratio determined from depth h and width w, and no basal slip. At cross-sections where section shape was known from seismic reflection measurements (B, C, D, E, EF, and F locations in longitudinal profile plots), f was also calculated from the actual shape (Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others, 1977, figs 3 and 7) using finite elements. The difference caused by deviation from parabolic shape was small (<6%) except at two sections near the tributary (EF ∼ 30%, F ∼ 11%).

The calculation of τ_b in the vicinity of the tributary between EF and F is, in any case, problematic, since the cross-section shape is not clearly defined at the confluence and the effect of the entering ice on the longitudinal force balance could be significant. Similarly, the surface and bed geometry above Km 3 are complex (MacQueen, unpublished) and Equations (3) are probably inadequate to capture the details of the distribution of τ_b there. In spite of these difficulties in estimating the full spatial distribution of τ_b we expect the broad features of the spatial pattern and the time trends to be correct.

The geometrical parameters α, h, w, and the resulting f all varied significantly with time and contributed to the shear-stress changes. Our interest is to compare velocity and stress. Therefore, the stress was computed from the average of the initial and final geometries of a survey interval, in order to approximate better the geometry responsible for the average velocity over that interval. The results for each winter interval are shown in Figure 7a.

Above Km 3, τ_b decreased with time, or the changes were small and unsystematic. Between Km 3 and Km 6, τ_b increased progressively from 1973 to 1979. From 1979 to 1982 the increase stopped, because increase in depth was compensated by decrease in surface slope. Between Km 6 and Km 13 there was a progressive increase over the complete interval. The largest changes occurred between Km 6 and Km 9, where both thickness and slope increased substantially with time (Fig. 4). Below about Km 13, τ_b decreased slightly or remained roughly constant where a decrease in depth was nearly compensated by an increase in surface slope.

Normal Stress

The basal normal stress σ _b beneath the center line may be estimated from

(4)

Results for winter intervals are shown in Figure 7b. Corresponding effective normal-stress distributions are considered in section IV.B based on a theoretical calculation of water pressure.

C. Deformation of the Ice Column

Depth-averaged shear strain-rate is given by u _d/h To estimate it, we assume zero sliding velocity (u _b = 0), so that u _d = u. The errors from this assumption should be smallest in winter, when u and presumably u _b were at their seasonal minima. Therefore, we used the velocity measured in winter. Because the ice depth h is slowly varying over most of the glacier length (Fig. 1), the pattern of variation in shear strain-rate appears rather similar to the pattern for winter velocity (Fig. 5b). When estimated in this way, a maximum average shear strain-rate of 1.8 × 10⁻³ d⁻¹ occurred in the vicinity of Km 4 during the winter of 1980–81.

Fig. 7.

Evolution of center-line basal shear (a) and normal (b) stress estimated using Equations (3) and (4). and geometrical parameters described in the text.

For comparison, the pattern of longitudinal strain-rate (∂_u/∂_x) during the winter seasons can be inferred from Figure 5b. A maximum magnitude of longitudinal strain-rate of 0.15 ×10⁻³ d⁻¹ occurred in the vicinity of Km 8 in the winter of 1980–81. At all locations and winter seasons, the longitudinal strain-rate was less than 0.2 times the depth-averaged shear strain-rate. Under these circumstances, softening of the ice column by coupling of the longitudinal strain-rate or corresponding deviatoric stress component into the effective non-linear ice viscosity should have a negligible effect on ud(Reference Raymond and ColbeckRaymond, 1980). Therefore, Equation (2) would not be expected to be invalid because of longitudinal strain effects.

Fig. 8.

Depth-averaged shear strain-rate estimated as winter velocity (u_w) divided by ice thickness (h) versus base shear stress (Fig. 7a). Points are computed for every ¼ km and every year 1973–81 for which data are available. Line represents prediction from flow law recommended by Reference PatersonPaterson (1981) for temperate ice (n = 3. A = 5.3 × 10⁻¹⁵ k Pa⁻³ s⁻¹ = 0.167 bar⁻³ a⁻¹). Curves join data points from fixed locations indicated by Km positions. Circles indicate data points from winter 1973–74; stars indicate dala points from winter 1980–81. Data points from intervening winters are not distinguished by distinct symbols in order to maintain clarity where density of data points is high.

These ideas are tested in Figure 8 by plotting u/hversus τ_b using values calculated every integral ¼ km for winters 1973–74 through 1980–81. The distribution predicted from Equation (2) using the flow law recommended by Paterson (1981; T = 0°C, A - 0.167 bar⁻ⁿ a⁻¹, n = 3, K = 2.28 × 10⁻⁴ bar⁻³ d⁻¹) is also shown for comparison. The overall pattern of data points is reasonably consistent with Paterson’s recommended parameters. Notice that most points from the winter 1973–74 lie below the expected distribution, while the data points from the winter 1980—81 lie above it. Additionally, points for which τ_b ≥ 1 bar suggest a slope steeper than n = 3.

The curves in Figure 8 link data from different years at selected fixed locations representative of the upper and lower parts of the glacier. Such curves, for all locations, cut across the spatially determined pattern of data points with slopes that are greater than correspond to n = 3. Steep slopes are especially noticeable on the lower part of the glacier at all locations below Km 10, where τ_b and u remained roughly constant and h decreased, thus giving increased u/hat constant τ_b and a nearly vertical trace on the plot (e.g. curve for Km 14 in Figure 8). On the upper part of the glacier at locations between about Km 3 and Km 7, there appears to be an increase in slope beginning in 1977 (e.g. curve for Km 5 in Figure 8). Figure 8 shows that changes in velocity over time at a fixed location were distinctly anomalous in comparison to the spatial pattern.

D. Velocity Anomaly

The parameter Kin Equation (2) can be calibrated from measured values of u, τ_b, and hat a reference time and location, so that

(5)

This presumes that u _b = 0 for the reference measurement. The deformational velocity u _d at any other time and/or location is then predicted to be

(6)

The difference between measured surface velocity u and u _d predicted from Equation (6) defines a velocity anomaly u _a

(7)

In the next section we assume that any such unexplained velocity arises from basal sliding; this is implied by Equation (1).

To apply Equations (5)–(1), we took the currently accepted value of n = 3 (Reference HookeHooke, 1981; Paterson 1981) and calculated K at each location from the values of u τ_b and h measured there over the winter 1973–74. This procedure assumes that u _b = 0 at all locations in the winter 1973–74. The resulting values of K have a mean of 1.5 ×10⁻⁴ bar⁻³ d⁻¹ and show a spatial variation up to about 30% from the mean with above-average values between Km 5 and Km 9, near average between Km 9 and Km 11, and below average below Km 11.

The mean of K is smaller than the value of K = 2.3 × 10⁻⁴ bar⁻³ d⁻¹ predicted from Paterson’s preferred parameters. This difference is also evident in Figure 8 by comparing the data points for the winter 1973–74 with the line predicted from Paterson’s parameters. Thus, the glacier was apparently stiffer than what Paterson identified as typical temperate ice. Any actual sliding motion not taken into account by our no-slip assumption would make this difference even larger. That assumption, therefore, appears to be reasonable. In terms of the experimental flow law found byReference Duval Duval (1977), the mean value of K would be consistent with tertiary creep of temperate ice with a water content of about 0.1%.

The spatial variation of K corresponds to a deviation of data points for the winter 1973–74 from any single straight line of slope 3 in Figure 8. Such a spatial variation in apparent flow law is also evident in the analyses of the velocity in 1973–74 by Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others (1977, fig. 7), Reference BindschadlerBindschadler (1982, fig. 9), andReference Raymond and Harrison Raymond and Harrison (1987, fig. 2). It could arise from: (i) basal sliding, (ii) data errors in surface velocity, elevation, ice depth, and cross-section shape, (iii) inadequacy of Equations (2) and (3), or (iv) a real longitudinal variation in ice stiffness. We assume that (ii) is the dominant contribution and that the same errors influence Equations (2) and (3) in a similar way at all times so that the apparent value of K does not change with time.

The variation of u _dwith time predicted from Equation (6) at two fixed locations is shown by curves on Figure 6. Figure 9 shows u _a given by Equation (7) for the summer and winter seasons. These anomalies represent the part of the seasonal velocity history that is not explained by geometry changes figured by Equations (2) and (3), and calibrated by Equation (5) in the reference interval winter 1973–74. It is significant that by 1981 u _a may be 50% or more of u at locations between Km 3 and Km 9. The implications are discussed in section IV. First we discuss how u _a relates to two other analyses of velocity on Variegated Glacier.

E. Discussion of Deformation Calculations

The prediction of u _d expressed by Equations (2)-(6) is essentially the same as used by Kamb and Echelmeyer (1986) to examine the distribution of velocity measured on Variegated Glacier over the interval July 1977-July 1978. Our approach differs from theirs primarily in the calibration of K by choice of the reference values of u, τ_b, h and in Equation (6). Our choice emphasizes the development of anomalous velocity with time relative to the winter 1973–74; any peculiarities of the spatial pattern in the winter 1973–74 are folded into the analysis as a spatially variable K. Reference Kamb and EchelmeyerKamb and Echelmeyer (1986) were concerned with the spatial pattern of velocity at only one time interval (1977–78) and calibrated K at a single reference location (Km 9.5). This gave a value of K = 2.4 × 10⁻⁴ bar⁻³ d⁻¹ in good agreement with the Paterson value. With this spatially constant value of K, they predicted u _d(x), which matched the measured u(x) very well over most of the glacier length; the only significant differences occurred up-glacier from Km 5, which they attributed to localized sliding at rates up to 0.2 m d⁻¹ (Kamb and Echelmeyer, 1986, fig. 9).

Because of the differences in calibration of K, these two approaches are not directly comparable. Both predict anomalies localized on the upper part of the glacier. However, the spatial effect identified by Kamb and Echelmeyer was more narrowly confined and locally larger than the anomaly u _a identified here (Fig. 9), which extended down to Km 8 and below. This difference arises primarily because the spatial pattern of velocity in the winter of 1973–74 contained similar features to that of 1977–78 (Fig. 5b), so that the spatial effect identified by Kamb and Echelmeyer was calibrated away by our spatially varying K.

It is also of interest to note that the value of K used by Kamb and Echelmeyer to fit the distribution of velocity in 1977–78 predicts a velocity in 1973–74 that is larger than measured everywhere between Km 3 and 15, which corresponds to a distinctly negative velocity anomaly over all of the well-measured length of the glacier. The negative anomaly is equivalent to the earlier recognition that in Figure 8 almost all of the data points from the winter 1973–74 lie below the line corresponding to the value of K given by Paterson’s preferred flow law, which is nearly equal to that used by Kamb and Echelmeyer. It is also related to the behavior illustrated in Figure 8 that velocity changes more rapidly at a given location than would be predicted from the spatial pattern. Such a negative anomaly could not reasonably be explained by basal sliding.

Like Kamb and Echelmeyer (1986),Reference Raymond and Harrison Raymond and Harrison (1987) examined the predictions of Equations (2)-(6) with spatially constant K. They showed that the mean value of K identified here and n = 3 produce the best least-squares fit of prediction to observations for the winter 1973–74 and that no time-independent, spatially constant value of K could match acceptably the full data set from all of the years, which is illustrated here by Figure 8.

This discussion underscores some of the ambiguity of predicting ice deformation and the importance of an accurate definition of the ice-flow law and a means to determine ice-structure parameters that affect it.

A. Relationship of Velocity Anomaly to Basal Sliding

If Equation (6) gives a correct prediction of u _d, then Equations (1) and (7) imply

(8)

where u _b is the basal sliding velocity in the center of the glacier and u _a is the deformational velocity anomaly defined by Equations (2), (3), (5), (6), and (7). Figure 9 would then display the longitudinal variation of center-line sliding velocity. The validity of Equation (8) is obviously subject to uncertainty from various sources.

The main uncertainty concerns the proper values of the flow-law parameters, which is evident in the discussion in section III.Е. In addition, ice stiffness could possibly change with time, in view of the effects of ice fabric and water content on the flow law (Reference DuvalDuval, 1977, 1981), and conceivable changes that could occur in the time-dependent flow environment of a surge cycle. We cannot entirely discount the possibility that softening of the ice contributed to the growing velocity anomaly u _a between 1973 and 1981.

It is also clear that there is a limitation on the spatial resolution with which u _b can be deduced. For example, if u _b varies longitudinally, longitudinal stress gradients are set up which damp the resulting variation in u at the surface (Reference Balise and RaymondBalise and Raymond, 1985). The damping is large for velocity variations on a longitudinal scale of 10 ice thicknesses or less (3–4 km on Variegated Glacier), so that features at this scale are not likely to be captured realistically by the estimate of u _bfrom Equations (6), (7), and (8). On Variegated Glacier, transverse velocity profiles (Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others, 1977) and the absence of cracking at the margins indicated sliding at the edges was negligible; therefore, if sliding occurred, it was more rapid in the center than at the edges. The resulting stress redistribution in comparison to the no-slip case causesf, τ_b, and u _d to decrease beneath the center line (Raymond, 1980) and Equations (3), (6), (7), and (8) would underestimate u _b there. Similarly, smaller-scale lateral variations in u _b would be attenuated at the surface. Therefore, the estimated value of u _b can at best represent a fairly large-scale spatial average over both length and width.

These various problems are not necessarily resolved by measurements of differential motion u _d in bore holes. Such measurements were made by Kamb and Engelhardt over intervals of several weeks during the early parts of each of the summers 1978, 1979, and 1980 in the vicinity of Km 7 on Variegated Glacier. Preliminary analysis of the tilts (personal communication from H.F. Engelhardt) indicates values of u _d plotted in Figure 6. A detailed analysis of measurement error has not been made but, based on the measurement interval of 3–4 weeks, an accuracy of approximately 0.05 m d⁻¹ or better is likely.

These bore-hole measurements clearly show that during the actual intervals of measurement (early summer) a large fraction of u _a is caused by sliding. On the other hand, the bore-hole measurements indicate a dramatic increase in u _d between 1978 and 1979 that is inconsistent with Equation (6) and also appear to suggest that during the winters of 1978–79 and 1979–80 u _a was not caused by sliding but rather by anomalous deformation. However, this interpretation is not straightforward. The holes made in different years were not at exactly the same locations and therefore may have sampled different local variations at the bed. The measurement intervals were short and included the occurrence of mini-surges (Reference KambKamb and Engelhardt, 1987) during which there could be substantial perturbation of the internal deformation pattern (Reference Balise and RaymondBalise and Raymond, 1985). Therefore, there is considerable uncertainty about how u _d and the corresponding u _b measured in these bore holes relate to the average over large areas of the bed and the long time intervals, which Equations (2)-(9) attempt to predict.

Fig. 9.

Velocity anomaly computed with n = 3 from Equations (6) and (7) for the summer (a) and winter (b) seasons.

From this discussion and with recognition of the seasonal change in velocity (Figs 5 and 6), it is fairly certain that the distribution of u _a determined for summer is correctly interpreted as being dominated by sliding. The interpretation of u _a determined for winter intervals is not clear. Since the spatial patterns of u _a in summer and winter have essentially the same profile shapes (with the exception of the summer 1981) and differ primarily in magnitude (Fig. 9), it is at least reasonable to assume u _a truly represents sliding during winter as well.

Even though there are interpretational problems, it is nevertheless important to investigate whether the estimated distributions of u _b.can be explained by proposed sliding laws (e.g.Reference Lliboutry Lliboutry, 1968, 1979;Reference Budd, Keage and Blundy Budd and others, 1979; Reference FowlerFowler, 1986) which predict u _b in terms of basal shear stress and water pressure. Toward this goal we will focus primarily on the winter distributions, in recognition that during summer the sliding velocity and underlying physical controls are tied to melt cycles dominated by short time-scale transients Reference Harrison, Raymond and MacKeith(Harrison and others, 1986; Humphrey and others, 1986; Raymond and Malone, 1986; Kamb and Engelhardt, 1987), which are not completely documented by data or easily approachable theoretically.

B. Control of Sliding Variations

Basal Water Pressure

Measurements of basal water pressure P _w are available only from a few locations and short time intervals during summers (Kamb and others, 1985; Kamb and Engelhardt, 1987). In order to extend this information in space and time, it is necessary to calculate the pressure.

The change in patterns of short-period velocity fluctuations and the termination of the major surge pulses in early summer suggest that by the latter part of each melt season the water flow had melted out a tunnel along the base of the glacier through which water flowed to the terminus stream (Kamb and others, 1985;Reference Harrison, Raymond and MacKeith Harrison and others, 1986; Humphrey and others, 1986). We shall assume that such a tunnel persisted during most of the winter and controlled the water pressure over wide areas of the bed. During the winter, the discharge was about 1 m³s⁻¹ at the terminus and transients in the flow should have been small. It is reasonable to estimate P _w in the tunnel, assuming steady state (Reference Raymond, Harrison and SenearRöthlisberger, 1972;Reference Spring and HutterSpring and Hutter, 1981). We assume that the tunnel is fed by water released from storage in the glacier and the ground, flowing to the tunnel through a distributed system of basal passageways under small lateral pressure gradients. Therefore, the steady winter discharge was assumed to vary spatially along the glacier length in proportion to the ice area up-glacier. This includes a step change in discharge at the tributary entrance to account for water from the tributary basin.

To calculate the distribution of P _w for each winter, the geometry (Figs 1 and 4) and overburden stress (Fig. 7b) were based on the average of the previous fall and following spring values. Ice flow-law parameters for tunnel-closure rate and tunnel-wall roughness were identical to those used by Reference Raymond, Harrison and SenearRöthlisberger (1972) to match pressure measurements on Gornergletscher. These particular flow-law parameters correspond to much softer ice properties than appropriate for shearing in the main body of a glacier as discusssed in section III. Whatever the cause of this inconsistency, for example, special conditions of stress or ice properties near the bed (Lliboutry, 1983), it promotes caution in applying the theory. Nevertheless, predictions from these tunnel-closure parameters do also appear to be consistent with the limited basal water-pressure measurements that are available from Variegated Glacier during June and July of 1978, 1979, and 1980 at locations near Km 7 (Kamb and Engelhardt, 1987). Water levels measured in the later part of July that most likely reflected the pressure in a basal tunnel and tunnel pressures predicted using the mean summer discharge at the terminus of about 16m³ s⁻¹ found by Humphrey and others (1986) both correspond to 14–15 bar below ice overburden. We assume the same tunnel parameters are also reasonable for winter.

Results for water pressure P _w in the form of effective normal stress (σ _b – P _w) are shown in Figure 11. There is a one-to-one correlation between the maxima and minima in the spatial patterns displayed in Figure 11 when compared to the distribution of ((σ _b – P _w) calculated byReference Bindschadler Bindschadler (1983, fig. 1) for Variegated Glacier in 1973 using assumptions similar to ours. Our results give somewhat lower mean levels of (σ _b – P _w), primarily because of a difference in assumed ice flow-law parameters. We shall focus on the relative changes in space and time that are insensitive to this assumption and are controlled by surface-slope and ice-thickness profiles. Figure 11 shows that above Km 14 variations in space and time of (σ _b – P _w) were less than 10% of the mean. This can be compared with the fractional changes in shear stress, which were up to about 50% (Fig. 7a). From this, one might expect the sliding variations to be controlled by τ_b and displayed by Figure 10. Nevertheless, the variations in (σ _b – P _w) may be significant, especially if u _bis very sensitive to (σ _b – P _w).

The small time-dependence of (σ _b – P _w) showed a gradual increase at locations down-glacier from Km 6 and a decrease at locations at Km 6 and above (Fig. 11). This boundary occurred at the peak in elevation change (Fig. 4). Down-glacier from the peak, the surface slope increased with time. Theoretically, as time passsed, the water discharge was driven out of the glacier with a higher head gradient, greater energy loss and conduit-melting rate, and a correspondingly faster creep-closure rate, which allowed a larger pressure difference between the flowing water and the ice. Above the peak in thickness change, the surface slope decreased with time and caused opposite effects. This direct relationship between surface slope and effective normal stress changes is a qualitative feature that we might expect to find for any model of basal water flow in which water pressure opposes closure of water-flow passageways. For example, it is predicted for water flow through linked cavities by Kamb (1987). These temporal trends in (σ _b – P _w) may Partially explain the difference in apparent sensitivity of u _b to change in τ_b above and below Km 6 (Fig. 10b). Below Km 6, the effect on σ _b by the increase in τ_b was opposed by increasing (σ _b – P _w). Above Km 6, the effect on u _b by increase in τ_b was augumented by a decrease in (σ _b – P _w).

Above Km 10 there was a gradual, although not monotonie, increase in (σ _b – P _w) with distance up-glacier. This trend is related to the up-glacier increase in surface slope (Fig. 1) and is explained qualitatively by the same reasoning as in the foregoing paragraph. Again, we may expect this qualitative pattern to be consistent with other models of basal water flow. This spatial trend may be a partial explanation of the lower u _b for given τ_b at the up-glacier locations.

These qualitative patterns in space and time indicate an inverse relationship between u _b and (σ _b – P _w).

Fig. 10.

a. Basal sliding velocity versus basal shear stress estimated at 0.25 km intervals on center line for summer and winter seasons 1973—81 (based on Figures 7a and 9).

b. Basal sliding velocity vesus basal shear stress for winter seasons 1973–81 (based on Figures 7a and 9). Curves link data points with same location.

Discussion of Sliding Behavior

Raymond and Harrison (1987) have used the data presented here in a search for a quantitative sliding law u _b = f(τ _b, σ _b − P _w) based on least-squares fit assuming several sets of possible flow-law parameters K and n. All of these choices yielded sliding laws such that u _b increased with τ_b and decreased with (σ _b – P _w), thus indicating that our qualitative conclusions appear to be robust and insensitive to details of assumptions about K and nHowever, none of the best-fit sliding laws were good fits to the data in the sense of significantly reducing the residuals in the overall data set below the level achieved by ice deformation alone with an intermediate ice stiffness, for example, Paterson’s preferred value for temperate ice. Said another way, if predicted sliding u _b were subtracted from measured surface velocity uto find corrected u _d in Figure 8, the scatter would remain essentially the same.

Fig. 11.

Effective normal stress computed using steady-state flow of water in a “Röthlisberger” tunnel, geometrical parameters in Figure 1, stress from Figure 8b, and water flux at the terminus of 1 m³ s⁻¹. Flow-law parameters (B and n in ε = (τ/B)ⁿ) and conduit Manning’s roughness (nʹ) were taken as B = 317 bar s^1/3, n = 3. Nʹ = 0.1 m^−1/3 s. The value of B corresponds to A = 3 × 10⁻¹⁴ K Pa⁻³s⁻¹ = 1 bar⁻³ a⁻¹ in the flow law έ = Aτⁿ.

It is possible that the relationship we seek is obscured by errors in the data and theory used to derive u _b, τ_b, and (σ _b – P _w). Certainly, we must suspect the calculation of basal water pressure. Perhaps the assumed basal tunnel did not exist over the full glacier length or the complete winter intervals (Humphrey and others, 1986) or the pressure on the adjacent bed was not tied to the pressure in the tunnel. If so, other water-flow models would be needed to calculate transport through linked cavities (Kamb, 1987), basal debris (Reference ClarkeClarke, 1987), or a combination that accounts for erosion and deposition of debris in the passageways. On the other hand, we expect the relationship between sliding rate and base stress to be affected by bed roughness or debris distribution, which could vary from place to place and with time. Whether the problem is caused by failure of the tunnel theory or non-existence of a spatially independent sliding law u _b = f(τ _b, σ _b − P _w), the resolution will apparently require knowledge of parameters related to detailed features of the bed structure. If this is true, it is effectively impossible to predict sliding rate without some calibration of the effects of the unknown bed structure.

C. Relationship to Surge Initiation

Despite the lack of success in discovering a “sliding law” that holds at all locations and times, the qualitative features are interesting. Of particular interest is the extreme sensitivity of u _b to small changes in profile above Km 6 in the few years before 1982. Such sensitivity is similar to a failure behavior that may have been related to the eventual surge initiation in 1982.

The theoretical ratio I = τ _b/(σ _b − P _w) is plotted in Figure 12.Reference Bindschadler Bindschadler (1983) referred to this ratio as the separation index, because as it rises, the probability of growth of basal cavities on the down-glacier sides of rigid bedrock bumps increases (Lliboutry, 1968; Kamb, 1970; Iken, 1981). More generally, we may expect that rising I increases the potential for rapid sliding independent of the actual basal structure. In view of the uncertainty in calculating basal water pressure and the corresponding effective normal stress, it is probably not realistic to ascribe much meaning to the actual numerical values of I in Figure 12. The spatial pattern, however, may be realistic. With this in mind, sliding rate appears to increase rapidly when I rises above a threshold that was reached in the few years before the surge between Km 4 and Km 7 (Figs 9b and 12). This location corresponds to the up-glacier flank of the growing topographic bulge centered between Km 6 and Km 7. Although the increase in τ_b (Fig. 7a) contributed strongly to the progressive rise in I, the drop in (σ _b – P _w). (Fig. 11) was the dominant influence just before the surge. This is of interest with regard to the location of the first rapid surge motion, which was predominant near Km 4. The surge initiation is discussed in more detail in a separate paper in preparation.