Introduction

In what Reference ClarkeClarke (1987) has described as the first example of numerical modeling in glaciology, Reference NyeNye (1965) obtained finite-difference solutions for steady-state flow in uniform, inclined channels of various geometries using Glen’s nonlinear flow law for ice deformation and a boundary condition of zero or uniform sliding. Although Nye’s predictions of cross-sectional flow patterns were consistent with data for glacier-surface velocities, subsequent observations of the internal distribution of flow in a glacier section (Reference RaymondRaymond, 1971; Fig. 1) were significantly different from Nye’s predictions, primarily due to the existence of marked lateral variations in basal sliding not accounted for in Nye’s original boundary-condition assumption. At Athabasca Glacier, Reference RaymondRaymond (1971) measured basal sliding rates up to about 40 m a⁻¹ at the center of the glacier (equivalent to about 80% of the surface velocity at the center line), but observed decreasing basal sliding rates away from the center line, dropping to less than a few metres per year at the margins. As Raymond pointed out, this observation has important implications for the distribution of velocity and stress in the cross-section, and the “theoretical analysis of flow in valley glaciers may be in systematic error until such lateral variation in sliding velocity is taken into account.” (Reference RaymondRaymond, 1971, p. 76).

Raymond’s observations prompted Reference ReynaudReynaud (1973) to expand upon Nye’s numerical approach by introducing a basal boundary condition that would allow for lateral variations in basal sliding. In Reynaud’s model, friction between ice and rock is controled by the effective pressure at the bed (which varies laterally), and with this Reynaud was able to produce results that matched the general characteristics of the empirical data from Athabasca Glacier. However, there are two limitations in Reynaud’s approach: first, to match Raymond’s empirical data required that Reynaud assume a piezo-metric surface position significantly different from that observed by Raymond; and, secondly, Reynaud’s approach required that the marginal velocity be specified in advance, which limits its usefulness for applications in which marginal velocities need to be predicted.

As part of a larger study of the long-term development of glacial landforms (Reference HarborHarbor, 1990b), I have investigated the application of a general sliding law to the problem of simulating flow in a glacier cross-section. Some initial results of this application were outlined in Reference Harbor, Hallet and RaymondHarbor and others (1988), although the emphasis there was not so much on the flow pattern as on the implications of this pattern for erosion distributions and the long-term development of glacial valley cross-sections. In this paper, I present in more detail the results of the flow-simulation modeling, in particular as they relate to constraining the form of a general sliding law.

A General Sliding Law

In attempting to explain marked lateral variations in sliding velocities observed at Athabasca Glacier, Reference RaymondRaymond (1971) observed that a simple relationship between sliding and shear stress was unsatisfactory because at the margin of the glacier there were low sliding velocities but high shear stresses. This led Raymond, following Reference LliboutryLliboutry (1968), to suggest that variations in sliding might relate not only to variations in shear stress but also to variations in effective pressure. More formally, Reference RaymondRaymond and Harrison (1987) tested a general law for basal sliding over a wet bed of the form:

(1)

where U _b is the basal sliding velocity, τ_b is the basal shear stress, is the effective normal stress at the bed, m and p are positive constants, and k is a sliding parameter related to bed roughness and structure. It is important to note that for the case of a glacier cross-section, assuming a horizontal piezometric surface (Reference BindschadlerBindschadler, 1983), the distribution of N across the section includes a maximum at the intersection of the piezometric surface with the glacier bed (Fig. 2; Reference HarborHarbor, 1990a). Under Equation (1), with positive values of p, flow is retarded in the region of maximum Ν which tends to produce relatively small marginal sliding velocities (Reference RaymondRaymond, 1971). Equation (1) is consistent with most previous experimental and theoretical analyses of glacial sliding (e.g. Reference WeertmanWeertman, 1964; Reference LliboutryLliboutry, 1968, Reference Lliboutry1979; Reference BuddBudd and others, 1979; Reference FowlerFowler, 1986), although estimates of the values of m, p and k have varied widely in past work (Reference PatersonPaterson, 1981; Reference Raymond and HarrisonRaymond and Harrison, 1987). For example, Reference Raymond and HarrisonRaymond and Harrison (1987) used observed spatial and temporal variations in flow along the length of Variegated Glacier to show that, although the sliding rate varied directly with shear stress and inversely with estimated effective normal stress, there was no single combination of values of m and p that could explain both the spatial pattern and time evolution of flow.

The Numerical Model

To apply the general sliding equation to flow in a glacier cross-section, I use a finite-element scheme that calculates patterns of ice flow through glacier cross-sections, and is set up to model the flow of a power-law fluid in an inclined channel of arbitrary cross-section, assuming the flow is rectilinear (i.e. flow in a straight channel with no convergence with, or divergence from the bed or ice surface). The calculations are based on a variational principle equivalent to the differential equations governing incompressible creeping flow, and the model can handle a variety of basal boundary conditions, including a wide range of possible sliding laws. The details of the finite-element model have been described in Reference RaymondRaymond (unpublished).

For this application, it is assumed that the ice is at its melting point throughout, and that deformation is described by Glen’s flow law:

(2)

where τ and are the shear stress and shear strain rate, respectively, n is a constant, and β depends on ice temperature, crystal size and orientation, impurity content, and possibly other factors (Reference PatersonPaterson, 1981). In the simulations described here, n is taken to be on the order of 3 and β on the order of 2 bar a^1/n, values appropriate for a temperate glacier (Reference Weertman, Whalley, Jones and GoldWeertman, 1973; Reference HookeHooke, 1981). The basal boundary condition is the sliding relationship described in Equation (1), which relates the basal velocity to the basal shear stress and effective normal pressure. One of the aims of the analysis presented here is to determine appropriate values for the exponents m and p in the sliding law.

Initial Results

In an initial attempt to model flow in a glacier cross-section, Reference Harbor, Hallet and RaymondHarbor and others (1988) used a sliding law based on Reference BuddBudd and others’ (1979) empirical work (which in terms of Equation (1) set m = 3 and p = 1) with a horizontal piezometric surface 20% of the maximum ice depth below the ice surface, and set k to yield a situation in which basal sliding accounted for 80% of the surface velocity at the center of the glacier (in an attempt to match the general characteristics of Reference RaymondRaymond’s (1971) empirical data from Athabasca Glacier). Although this form of the sliding law provided an acceptable first attempt at modeling cross-sectional flow patterns, it was unable to reproduce one key element of the empirical data: the model results included marginal velocities on the order of 60–70% of the maximum surface velocity (Reference Harbor, Hallet and RaymondHarbor and others, 1988), whereas Reference RaymondRaymond’s (1971) data (Fig. 1) indicate negligible marginal velocities (less than 10% of the maximum surface velocity) for a case where the side walls are gently inclined. Although there are wide variations between glaciers, borehole experiments indicate that basal sliding may on average account for about 50% of the surface velocity at the center of temperate glaciers (Reference KambKamb, 1964), whereas marginal sliding velocities are generally considerably less: Reference MeierMeier (1960) found marginal sliding velocities of less than 20% of the center-line surface velocity for the Castleguard section of Saskatchewan Glacier (Reference RaymondRaymond, 1971), and Reference Engelhardt, Harrison and KambGlen and Lewis (1961) measured marginal sliding velocities of 10–20% of the maximum surface velocity on Austerdalsbreen in the summer of 1959. However, Reference GlenGlen (1958) reported measurements made at the base of an icefall on Austerdalsbreen which gave marginal sliding velocities up to 65% of the maximum surface velocity, suggesting important spatial variations in marginal sliding velocities as a function of local topographic and glaciologie conditions.

Fig. 1. Velocity distribution for a cross-section of Athabasca Glacier, based on surface-velocity measurements and borehole deformation at locations 1A through 5A. Velocity contours are in m a⁻¹. Modified from Reference RaymondRaymond (1971).

Sensitivity to Variations in the Sliding-Law Exponents

Clearly, it is important to investigate why the initial modeling failed to produce low marginal sliding velocities. Although it is possible that the stress exponent in the flow law (n in Equation (2)) might decrease in low-stress areas near the margin (see Reference BuddBudd (1969) and Reference PatersonPaterson (1981), but also Reference WeertmanWeertman (1969) and Reference ThomasThomas (1973)), this does not explain the elevated marginal velocities in the initial results. In the initial simulations, basal sliding accounted for almost all the motion near the margins, thus any difference between the model predictions and the empirical data in this region presumably relates to the sliding law rather than possible non-Glen flow-law behaviour at low stresses.

To investigate how variations in the sliding law might change the model predictions, a number of simulations were performed for flow in a parabolic channel of half-width 620 m, depth 310 m, and surface slope 3.5° (the approximate characteristics of the Athabasca Glacier sections studied by Reference RaymondRaymond (1971)). Internal deformation was modeled using Glen’s flow law with n = 3 and 1.5 < β < 3.0 bar a^1/n (see below), and Equation (1) was applied as the sliding-law boundary condition. Various combinations of m and p were tested to assess the sensitivity of the flow pattern in a glacier cross-section to variations in the sliding-law parameters.

In all of these simulations, the aim was to produce results with center-line surface velocities close to 50 m a⁻¹and center-line bed velocities close to 40 m a⁻¹ (again matching the general characteristics of the Athabasca Glacier data). Fixing these velocities proved to be straightforward: for all the combinations of m and p investigated, it was possible to find reasonable values of β and k which, when applied to the whole cross-profile, gave the desired center-line velocity values. However, to match exactly the empirical data across the whole section, one would presumably have to experiment with spatial variations of β, n and k for each sliding law investigated and, considering the limitations of the empirical data, it is unlikely that such tedious attempts to match exactly the field data are likely to provide significantly better insights into the controls on the flow of ice in glacier cross-sections. For these reasons, the simulations were initially restricted to spatially uniform values of β, n and k, although β and k were varied between simulations to match target center-line velocities.

By controling center-line conditions, it is possible to observe how changes in the sliding-law exponents affect velocity patterns across the rest of the channel. The results for a situation with a horizontal piezometric surface 0.13 of the maximum ice thickness below the surface (Ps = 0.13, where Ps is the distance from the ice surface to the piezometric surface, scaled relative to the maximum ice thickness) are illustrated in Figure 3, which shows surface and basal velocity variations across half of a symmetric glacier. In these plots, U* is dimensionless velocity (velocity at a point expressed as a percentage of the maximum velocity in the section) and W* is dimensionless distance from the valley center line in the horizontal direction (distance divided by the maximum ice thickness). In dimensionless units (W*), the glacier thickness is 1.0 at the center line and the piezometric surface is 0.13 below the ice surface. For comparison, Figure 3 includes a similar plot of U* for Athabasca Glacier, derived from Raymond’s reconstruction of the velocity field. The piezometric surface in this case (Ps = 0.13), taken from Reference RaymondRaymond (1971), yields very low effective normal pressures at the center of the glacier (0.6 bar), and maximum effective pressures (3.5 bar) where the piezometric surface intersects the bed.

Fig. 3. Sensitivity of bed- and surface-velocity distributions to variations in the sliding-law exponents m and p, and comparison with Reference RaymondRaymond’s (1971) data from Athabasca Glacier. W* is dimensionless distance away from the center of the glacier (distance divided by maximum ice depth) and U* is velocity expressed as a percentage of the maximum velocity for the section. For Athabasca Glacier, sliding velocities are well constrainedfor −1 < W* < 1, and at the margins of the section.

With increasing values of m in the sliding law (a strong dependence on τ_b), marginal velocities increase (Fig. 3), which conflicts with observed, generally low marginal sliding velocities. Because shear stresses are high near the margins of glaciers, a sensitive dependence on τ_b in the sliding law (large values of m) is likely to give high marginal sliding velocities. However, with increasing values of p (a sensitive dependence on effective pressure), progressively lower marginal velocities are predicted (Fig. 3), and with p = 5 essentially negligible sliding velocities. This result arises because effective pressures increase from the center of the glacier to the point where the piezometric surface intersects the bed (Fig. 2). As the sliding velocity is inversely related to effective pressure (Equation (1)), higher effective pressures away from the center of the glacier give rise to lower sliding velocities. With increasing values ofp, this effect is enhanced, so marginal velocities fall. However, the simulations also predicted reversals in the basal velocity field towards the edge of the glacier (above the piezometric surface) which arise because close to the glacier margin Ν decreases with decreasing ice thickness, so with positive values of p this gives increasing sliding velocities towards the edge of the glacier.

Fig. 2. Effective pressure distribution for a glacier with a parabolic cross-section: note that the maximum effective pressure for the section occurs where the piezometric surface intersects the bed. W* is dimensionless distance away from the center of the glacier (distance divided by maximum ice depth), and pressures are scaled relative to the pressure at the base of a column of water of depth equal to the maximum ice thickness for the section.

Thus, with a sliding law sensitively dependent on basal shear stress, the model was unable to produce low marginal sliding velocities, whereas with a sensitive dependence on effective pressure the model predicted low marginal velocities, but also basal-velocity reversals towards the edge of the glacier section. There are as yet no data with which to test the prediction of a reversal in basal velocities towards the margin but, despite the lack of data, it seems reasonable to assess the implications of assuming that this reversal is incorrect. This assumption allows for an assessment of whether relaxing any of the current assumptions used in the flow modeling can lead to results in which basal velocities decrease consistently towards the margin.

Sensitivity to Variations in the form of the Piezometric Surface

In the simulations described so far, it was assumed that the piezometric surface is horizontal and located 0.13 of the maximum ice depth below the surface (Ps = 0.13). If this assumption is relaxed, and the level of the piezometric surface lowered, velocities across the section drop and velocities at the margins progressively increase relative to those at the center (Fig. 4). This arises because, with lower piezometric surfaces overall, effective pressures are increased, reducing sliding velocities, and the point of maximum N on the bed (Fig. 2) is shifted away from the glacier margin and thus has less of a retarding effect on marginal flow. With the current version of the sliding law (Equation (1)), piezometric surfaces lower than a level close to that at which flotation would occur in the center of the glacier are unable to produce marginal sliding velocities that are low compared to those at the center. In the case close to flotation (Ps = 0.1), there is still a reversal in basal velocities close to the margin.

Fig. 4. Sensitivity of bed- and surface-velocity distributions to variations in the level of the piezometric surface (Ps). Sliding-law exponents in each case are m = 1 and p = 2.

In addition to the level of the piezometric surface, the form of the piezometric surface across a glacier appropriate for modeling sliding is also an unknown. Water levels in boreholes fluctuate widely over a range of time-scales, and little information exists on the form of the piezometric surface towards the margins of glaciers. By assuming that subglacial water flows longitudinally in a main channel at the low point of the cross-section, Reference BindschadlerBindschadler (1983) has shown that, if the piezometric surface is controled by the water-pressure distribution in hypothesized transverse secondary channels, then the piezometric surface will be approximately horizontal for cases where the main-channel pressure is relatively high. However, for low values of water pressure in the main channel, Bindschadler’s numerical results suggest that the secondary-channel pressure distribution is equivalent to a piezometric surface that has a minimum at the center of the cross-section, but which rises rapidly with distance away from the center to an approximately horizontal surface. This gives maxima in the distribution of N across the glacier section both in the center of the channel and at the point where the piezometric surface intersects the bed. A simulation with this type of piezometric surface gave low sliding velocities at the center of the glacier (because of the high values Ν), and maximum sliding part-way along the channel walls (Fig. 5). Clearly, this type of velocity distribution is incompatible with the empirical data from Athabasca Glacier, which give both high sliding velocities and low Ν in the center, although it could possibly apply elsewhere.

Fig. 5. Bed- and surface-velocity distributions for a case with a central minimum in the basal water-pressure distribution (Reference BindschadlerBindschadler’s (1983) low main-channel water-pressure case).

Finally, having noted that water levels in boreholes fluctuate widely over a range of time-scales, it is effective roughness at the margins of glaciers is relatively high. Simulations with this assumption will be referred to here as simulations which use an adjusted sliding law.

The tests of possible combinations of m and p were repeated with the adjusted sliding law and showed again that high values of m do not result in low marginal sliding velocities (Fig. 6). However, with m = 1 and p = 1, marginal velocities are low, and with increasing p the basal-velocity distribution becomes more peaked towards the center of the glacier and includes negligible velocities towards the glacier margin. It is important to note that with the adjusted sliding-law velocity reversals no longer occur close to the glacier margin, and the complete flow field looks similar to the empirical data (cf. Figs 1 and 7). With variable piezometric heights (Fig. 8), the model again predicts that relative marginal velocities increase as the piezometric surface is lowered relative to the glacier surface, and only in a case where the piezometric surface is nearly high enough to float the center of the glacier does the model predict low relative marginal velocities.

Fig. 6. Sensitivity of bed- and surface-velocity distributions to variations in the sliding-law exponents m and p, using the adjusted sliding law. W* is dimensionless distance away from the center of the glacier (distance divided by maximum ice depth) and U* is velocity expressed as a percentage of the maximum velocity for the section.

Fig. 7. Velocity solution for a parabolic cross-section using the adjusted sliding law with m = 2 and p = 2. The maximum ice depth in this case was 310 m and the surface width 1240 m, the approximate dimensions of the Athabasca Glacier section shown in Figure 1.

Fig. 8. Sensitivity of bed- and surface-velocity distributions to variations in the level of the piezometric surface (Ps) using the adjusted sliding law. Sliding-law exponents in each case are m = 1 and p = 2.

important to recognize that in the model results presented here it is assumed that the piezometric surface is effectively steady-state. The simulations do not include the effects of temporal variations in basal water pressures. However, field observations have shown that sliding velocities respond quickly to variations in basal water pressures over time-scales ranging from hours to seasons (e.g. Reference Boulton and VivianBoulton and Vivian, 1973; Reference HodgeHodge, 1976), and this lends support to the model assumption that sliding velocities reflect current water-pressure conditions, rather than some complex, integrated response to conditions over a longer time-scale. Thus, it seems reasonable to assume that the piezometric surface is steady-state for this initial sensitivity analysis.

Sensitivity to Variations in Basal Friction

Another major assumption of the simulations presented so far is that k, the friction factor in the sliding law, inversely related to bed roughness and debris concentration, is uniform across the glacier section. However, there may be important cross-sectional variations in the friction factor: one might expect effective roughness at the margins of many glaciers to be anomalously high because of relatively high concentrations of basal and englacial debris from subglacial and subaerial sources (Reference Kamb, Pollard and JohnsonKamb and others, 1976; Reference Engelhardt, Harrison and KambEngelhardt and others, 1978; Reference Harbor, Hallet and RaymondHallet, 1981). Although it is not clear exactly how to treat this problem quantitatively (Reference Harbor, Hallet and RaymondHallet, 1981), it is mathematically convenient to let k scale with N_p towards the edge of the glacier (above the piezometric surface) so that spatial variations in the product of k and N^−p vanish (equivalent to using a Wecrtman-type sliding law in the marginal zone (Reference WeertmanWeertman, 1964)). As k is inversely proportional to effective roughness, and decreases in the marginal zone, this adjustment accords with the expectation that the

Simulation of the Athabasca Glacier Flow Field

To test the ability of the flow model to replicate the empirical data in detail, using the adjusted sliding law, flow was simulated for the cross-section of Athabasca Glacier shown in Figure 1. Subjectively, an acceptable match between the empirical data and the simulations is provided by m = 2, p = 2 (Fig. 9), in particular the results replicate the most important elements of the cross-sectional variation in basal velocities. Where the empirical data are most reliable (solid lines in Figure 9A), the model velocity contours provide a close match to the empirical data. However, it is important to note that in the interior of the flow in this region the shapes of the velocity contours are largely independent of the boundary condition and hence the sliding law. On the channel-side slopes, the agreement between the extrapolated empirical data and predicted velocities is somewhat poorer than in the central zone, and there are some clear differences between the observed and predicted velocity fields. Most noticeably, the model is unable to reproduce the asymmetric velocity field of the empirical data, although this is not surprising as it probably relates to curvature in the long-profile of the glacier which is not accounted for in the rectilinear flow model we are using. In more detail, the velocity contours extrapolated from the empirical data (dotted lines on Figure 9A) suggest a zone of rapidly decreasing basal velocities on the lefthand side of the channel which is not matched in the model results, while in the vicinity of the bedrock shoulder on the righthand side of the channel the model results suggest a more rapid drop off in basal velocities than is indicated by the empirical data. Nonetheless, the most important elements of the cross-sectional variation in basal velocities are predicted. More refined comparisons between theory and observations would require more empirical data for flow in glacier cross-sections, in particular better data for sliding and water-pressure variations towards the margins of glaciers, and for spatial variations in bed roughness.

Fig. 9. Comparison of empirical data (A) for the flow field through a section of Athabasca Glacier with model predictions for a similar cross-section form (B) using m = 2, p = 2 and the adjusted sliding law.

Constraints on the Sliding Law

Although the form of a general sliding law has previously been evaluated with respect to temporal and spatial variations in the controling variables along the central region of a glacier (Reference Raymond and HarrisonRaymond and Harrison, 1987), the attempt to predict sliding velocities across a glacier section provides another way to assess the form of this law. Clearly, on the basis of a comparison with only one data set, it would be inappropriate to argue that m = 2 and p = 2 in the adjusted sliding law represent the most appropriate values for the general case, yet Reference RaymondRaymond’s (1971) data are the only complete empirical analysis of flow in a cross-section that is available. The values of m and p suggested by the results presented here, however, vary markedly from the range given in Reference Raymond and HarrisonRaymond and Harrison (1987) as providing the best fit for spatial and temporal variations in sliding along the center of Variegated Glacier. In particular, Raymond and Harrison’s analysis suggested that m is on the order of 5, a value which in the cross-section application is unable to provide the low marginal velocities observed in the cross-section data. If m is on the order of 5, and does not vary spatially, to produce low marginal velocities would require very significant lateral variations in effective bed roughness, which may not be unreasonable, although there are presently no empirical data with which to assess the existence of such a pattern. Interestingly, even with low values of m, some enhancement of k is required near the margin to produce low velocities, to offset the effect of vanishing effective pressures at the margin. Clearly, there is a pressing need for more complete theoretical and empirical analyses of flow in cross-sections to address the precise nature of controls on sliding close to the margins. As Reference Raymond and HarrisonRaymond and Harrison (1987) concluded, we are not yet able to determine the most reasonable values of the sliding-law parameters from the available data. However, even if the current form of the sliding law is appropriate, considering the range of basal conditions that can exist, it seems unlikely that there exists a single, spatially invariant set of parameters that will fit all cases.

Conclusions

Observations at Athabasca Glacier and elsewhere suggest that basal sliding can account for a significant part of total flow, and that sliding rates vary significantly across a glacier section. The ability to model such spatial variations in basal velocities is important in understanding flow in valley glaciers, as well as in predicting spatial patterns of glacial erosion that drive land-form development models.

With a sliding law in which the basal velocity is dependent on the basal shear stress and inversely dependent on the effective pressure at the bed, it is possible to predict an overall flow pattern that is consistent with the empirical data, if it is assumed that friction increases close to the margin of the glacier. In particular, it is possible to predict high sliding velocities in the central region of the cross-section, which decrease progressively to relatively low values at the glacier margin. To model flow in glacier cross-sections more precisely will require further effort to constrain the form of the sliding law, and in particular it will require further field efforts to establish the nature and scale of cross-sectional variations in surface roughness and basal water pressures, especially close to the margins. In particular, it is surprising that to date there exists only one empirical study of the flow pattern in a glacier cross-section (Reference RaymondRaymond, 1971), whereas several such studies over a range of glaciological conditions is needed to further constrain the form of the sliding law.