Stress and velocity fields in glaciers: Part II. Sliding and basal stress distribution

Abstract

Numerical methods are used to examine the interaction between the spatial distribution of the basal shear traction and the corresponding basal velocity for an inclined slab geometry. In our improved treatment, we reject the common assumption that basal velocity is a simple function of local variables in favour of a non-local treatment that includes normal deviatoric stress and takes basal velocity to be an integrated response to spatially varying influences. Computationally, one must either iterate the basal velocity with a friction parameterization that relates basal shear traction to basal velocity or, alternatively, prescribe the basal shear traction that results from bed decoupling and substrate déformation.

The average of basal shear traction over the entire bed of the ice mass is invariant under changes in sliding distribution and thus constitutes a useful reference; any local relative reduction of traction leads to basal movement, either sliding over the bed or moving with a deforming subglacial layer. The local stress réduction is accompanied by a concentration of traction up-and down-glacier of the moving base. Growth, decay and possible migration of basal stress concentrations may be closely related to short-lived sliding events and to surges.

1 Introduction

Thermomechânical modelling of glaciers and ice sheets has become an indispensable tool for the interpretation of field observations. Although thermal conditions, strain rate and the general behaviour of an ice mass can be modelled to acceptable accuracy, the flow contribution of sliding still awaits satisfactory treatment. Most studies that deal with sliding have concentrated on the micro-scale, for example, bed protuberances and thermodynamics (Weertman, 1979; and references therein), bed waviness with or without partial cavitation and water pressure (Iken, 1981; Gudmundsson, 1997a, b; and references therein). Common to most of these studies is the assumption that the representative roughness element is spatially periodic.

The application of friction laws in the shallow-ice approximation obscures the crucial point that in glaciers the local flow and basal velocity is not locally defined. Friction laws can be applied in higher-order flow models, but the higher-order equations must be applied to new and corrected basal boundary conditions that must be iteratively met. The application of a friction law also feeds back to the stress field and thus must be included in the iteration (Hutter, 1983, chapter 4). Such friction laws constitute a functional relation between basal velocity and other basal conditions, such as stresses, water pressure and roughness elements with or without cavitation.

The treatment of sliding in a numerical model faces the problem that sliding is determined by different processes which act on various spatial scale lengths. A numerical model with a given grid size must handle average basal conditions that are not defined by the microscopic sliding law alone. On the micro-scale, a Coulomb-type sliding (Lliboutry, 1968) occurs with a more or less velocity-independent friction coefficient. Sub-grid variability of the friction coefficient and resistance to sliding due to sub-grid topographic variations, or both together, must be parameterized for each gridpoint. However, such sub-grid patterns may lie unlimited in their variety and, thus far, only a few simplified situations have been investigated theoretically in detail (Gudmundsson, 1997a, b; and references therein).

For frictionless sliding acting on the true base, a periodic topographic variation with sinusoidal shape seems to give a power-law relation between average basal velocity and average basal shear traction (Gudmundsson, 1997b). For this type of sliding parameterization, both sliding velocity and basal shear traction are non-locally defined variables in the sense that both depend not only on local basal conditions but also on the conditions in neighbouring and more remote areas.

Sediment beds may be smoother than rock beds and may be flat on a scale length of modeled gridcells. In such cases, the variability of the Coulomb friction defines the sub-grid parameterization of sliding. Friction coefficients are then not dependent on the basal velocity but may vary Spatially and temporally with changing basal roughness and hydraulic conditions. In these cases, basal shear traction is locally defined by local conditions but the basal velocity is still a non-local variable.

Friction and sliding laws themselves are not the topic of the present paper. This study investigates the relation between the spatial pattern of basal shear traction and the spatial pattern of basal and internal shear with emphasis on “spatial pattern” for topographically flat beds having spatial Journal of Glaciology variations in the friction coefficient. The results demonstrate the important influence of deviatoric stress gradients on glacier flow in general and on the sliding pattern in particular. The range of influence of spatial inhomogeneities at the bed is of the same order of magnitude as the horizontal extent of average valley glaciers. This shows that local measurements of sliding must be interpreted with a view to the mechanics of the entire glacier.

This paper makes explicit use of methods presented by Colinge and Blatter (1998).

2 Temporal and spatial variability of ice-bed coupling

Temporal variability in glacier-flow rate occurs at time-scales that range from seconds to centuries (Forbes, 1851; Meier and Post, 1969; Iken, 1981). Temporal variation in sub-glacial water pressure is the usual suspected agent of this variability.

Spatial variability is more challenging to observe, because it originates at the bed and is diffusively projected on to the upper boundary (Budd, 1971; MacAyeal, 1992; MacAyeal and others, 1995). Recent subglacial measurements of water pressure (Murray and Clarke, 1995), sliding (Blake and others, 1994), ploughing (Fischer and others, submitted to Journal of Glaciology) and subglacial deformation (Iverson and others, 1995; Hooke and others, 1997) confirm that spatial variability is significant over length scales down to or below ~1 m Figure 1, showing subglacial observations from Trapridge Glacier, Canada, illustrates this point. Figure la shows summer 1992 water-pressure measurements from four sensors that have a typical spatial separation of ~5 m. Sensor 92P06 is situated in a hydraulically connected hole and registers a strong diurnal cycle with peak water pressure occurring in late afternoon. The water-flotation level in this region of the glacier is roughly 63 m, so that several of the most pronounced peaks correspond to pressures sufficient to cause artesian outflow from the subglacial water system. Nearby sensor 91P14 displays the same diurnal oscillations as 92P06 but with attenuated amplitude. In contrast, sensor 91P13, which is situated in a region of the bed that is not hydraulically connected to 92P06, shows water-pressure fluctuations of opposite polarity to those for 91P14 and 92P06; sensor 92P09, which seems to represent an intermediate case, alternates between in-phase and out-of-phase response. An explanation of this behaviour, involving the transfer of the ice-overburden support From one point to another, has been proposed by Murray and Clarke (1995). The main lesson, however, is that pronounced spatial and temporal variations in subglacial water pressure can occur.

Because ice-bed coupling is strongly influenced by sub-glacial water pressure, it is no surprise that the sliding and subglacial deformation rates are also spatially heterogeneous. Figure lb shows the summer 1996 measured sliding rate at two contiguous sites at the bed of Trapridge Glacier. Again, the spatial separation is ~5 m and a diurnal cycle can be discerned. The most striking feature of these graphs is that for several cases the peak sliding rate for 96SL01 corresponds to the minimum sliding rate for 96SL02. The obvious interpretation is that reduced ice-bed coupling (and thus a large sliding rate) at one site is compensated by increased ice bed coupling at other sites in order to maintain overall stress balance. Simultaneous measurements of subglacial sediment deformation rates at nearby sites lead to similar conclusions. Figure le shows results from tilt cells 96BG03 and 96B05 which are located within the same array as the sliding sensors. Diurnal variability is present but there is no correspondence between peak deformation rates for the two sensors. In summary, the observations beneath Trapridge Glacier reinforce the central theme of this paper: that the stress field in the vicinity of the ice-bed contact is markedly heterogeneous and that this heteroge-neity adds unwelcome complications to sliding mechanics.

Fig. 1. Evidence for spatial heterogeneity of subglacial water pressure, sliding rate and subglacial deformation rate for Trapridge Glacier, Canada. (a) Measured subglacial water-pressure records at four contiguous sites. (b) Measured sliding rate at two contiguous sites. (c) Measured subglacial strain rate at two contiguous sites.

3 Basal stress and velocity

3.1. Governing equations and boundary conditions

In this section, we develop the relation between basal shear traction, sliding velocity and the spatial extent of soft spots in the simplest case of a parallel-sided slab of isothermal homogeneous ice. The derivations and discussion of the force-balance equations and stress-strain rate relations have been presented in Part I of this study (Colinge and Blatter, 1998).

The geometry of the ice slab is defined by the upper free surface S ≡ H and the basal surface B ≡ 0 in Cartesian coordinates (x,z), with the x axis aligned with the surface slope. The z axis and the direction of gravity include an angle ?. The non-dimensional force-balance equation and stress-strain rate relations have been presented by Colinge and Blatter (1998), but here we only summarize the relevant equations corresponding to the first-order approximation, which are

(1)

(2)

(3)

with the flow law

(4)

where f, σ and ù are the shear stress, deviatoric normal stress and the longitudinal velocity component, respectively, the constant rate factor A ₀,and is proportional to the inverse of the viscosity (= fluidity) at vanishing effective stress. The surface-boundary condition is f_s =0 and at the base a sliding law must be imposed.

The foregoing scaling reduces the number of free parameters for a homogeneous slab to two, namely and A _0. The rate factor A₀ is simply a multiplier for the velocity field and the stress field is independent of A₀, as long as A_0, is homogeneous over the whole domain. This permits normalization of the velocity field for different choices of and makes A₀ a function of ,reducing the number of free parameters to the single parameter . In all examples, we take =0.1 and A₀ = 5/3 such that the asymptotic surface velocity

(5)

The asymptotic stress and velocity profiles are defined as the solution for the entirely non-sliding slab, i.e. the solution in the absence of or far away from any sliding perturbations. With this assumption, the only freedom rests in the choice of the basal stress pattern over a chosen sliding area. In most of the subsequent examples, the basal shear traction was set to zero within the sliding zone. This is the most extreme case and all other possible distributions of shear traction within the sliding area arc intermediate between this extreme case and the non-sliding case. In this way, the relevant patterns of basal velocity and stress can be summarized with a relatively small number of model computations.

To solve Equations (1)-(3), we used the numerical method presented in Paper I for a parallel-sided slab. An infinitely long plane slab is simulated by connecting the lower end of the finite-model domain of the slab with its upper end.With uniform basal boundary conditions, this produces a homogeneous solution for the entire slab; with non-uniform boundary conditions, this yields a (non-uniform) periodic pattern with wavelength of the finite domain length in the longitudinal direction of the infinite slab. If the influence of some internal inhomogeneity decays to zero towards the end of the domain the model approximates an infinitely long slab. This is the case if the distance of the sliding part from the boundary of the computed domain exceeds roughly five times the ice thickness.

In Colinge and Blatter (1998) we recognized that the closure of the boundary-value problem is related lo the basal sliding problem. The simplest basal boundary condition is the non-slip condition Ũ_b= 0. In the following paragraphs, we investigate the relation between the spatial distribution of basal shear traction and the resulting spatial distribution of basal and interior ice flow. Only flat-bed situations without any sub-grid topographic variations are considered. Such situations may occur in fast-sliding ice streams and under glaciers sliding over sediment beds. In these cases, spatial variations in the friction coefficient may occur due to variations in basal hydraulics and local qualities of the sediment itself.

3.2. Flow over a single sliding spot

Although the englaeial stress field is affected by the distribution of basal movement, there is no possibility from surfieial observations alone of distinguishing whether this movement stems from sliding of the ice sole over the bed, from dragging an underlying layer of deformable material or both together. However, for convenience, the terms “sliding” and “non-sliding” are used in the sense of “Presence of basal movement” and “absence of basal movement”, respectively. The term “basal” refers to the sole of the ice, however sharply this can be defined.

In an initial series of numerical experiments, we calculated the stress and velocity conditions in the vicinity of an isolated sliding spot in an otherwise homogeneous non-sliding slab. Within the sliding region, we prescribe a reduced basal shear traction

(6)

note that ?, ^ = 1 in the homogeneous non-sliding part. The reduction factor q = q(x) of the basal shear traction is a function of position in the sliding area and it represents the decoupling between ice and bed.

Only a small sample of situations is presented and these have been selected to illustrate essential features of the stress and velocity fields across sliding spots. The sliding conditions are specified at seven gridpoints (Fig. 2), for which perfect slip, i.e. vanishing basal shear traction is prescribed. Thus, the total width of the sliding zone is about 7Δx and the associated spatial patterns for stress and velocity are resolved with seven gridpoints.

Fig. 2. Illustration of the vertical gridlines and the corresponding labels in the sliding part of the slab, for which velocity and stress profiles are shown in figures 3 and 4.

The vertical profiles of longitudinal velocity û are shown in Figure 3 for Δx = 0.1, Δx = 0.25 and Δx = 0.5. The profiles along the vertical gridlines from the edge to the centre of the sliding spot are compared with the asymptotic profile. The example with the smallest spot size clearly demonstrates a bridging effect across the sliding spot due to normal stresses. Even with complete Uncoupling of the ice from the bed, when t_b= 0, the sliding velocity remains strongly limited. With growing spot size, the bridging diminishes and the sliding velocity increases. Figure 4 shows the corresponding vertical profiles of shear stress in the sliding area and in the transition zone between sliding and non-sliding in comparison with the asymptotic profile for Δx = 0.5. The sliding velocity with t_b = 0 is the largest possible for a given geometric pattern of sliding. Any larger value of basal velocity would result in negative basal shear traction, i.e. a shear stress acting to oppose the driving stress and produce inverse velocity profiles. Such a situation would be difficult to explain. These examples also demonstrate the extent to which the vertical profiles of deformation velocity and shear stress are changed by spatially variable sliding.

Fig. 3. Vertical profiles of longitudinal velocity above a sliding area in an otherwise non-sliding ice slab, as shown in Figure 2, together with the asymptotic velocity profile. The local reduction factors of basal shear traction are 0.5 at the gridimes labelled 1 and 0.0 at all other gridlines within lite sliding zone. The line labelled a shows the asymptotic velocity profile. Three examples are shown corresponding to grid sizes Δx =0.1, Δx = 0.25 and Δx=0.5 measured in units of one slab thickness. The dashed profiles correspond to the last non-sliding position (labelled m in Figure 2) adjacent to the sliding area; the dash-dotted, dash-triple dotted, long-dashed and solid lines correspond to the positions labelled 1, 2, 3 and c, respectively.

Fig. 4. Vertical profiles of shear stress in the sliding zone for Δx = 0.5. The curves correspond to the gridlines labelled in Figure 2: m (dotted line), 2 (dashed line), c (dash-dotted line) and to the asymptotic profile (solid line).

The longitudinal gradients within the transition zone between sliding and non-sliding are limited. These transition zones suffer the largest stresses in the basal laver and are therefore prone to structural failure, such as basal crevassing under high static pressure. Figure 5 shows the longitudinal profiles for σ,τ and the effective stress τ_eff= s/๠+ τ² for the five cases Δx = 0.1, 0.2, 0.3. 0.4 and 0.5. In all five cases, the basal shear traction is set to drop to zero over one mesh size. The weak dependence of the stress components on Δx suggests an approximate similarity rule: if the geometric patterns of basal shear traction in sliding areas are similar, then the resulting stress components and stress gradients in the transition between sliding and non-sliding are only weakly dependent on the spatial extent (represented by the grid size Δx) of the sliding areas. The resulting velocity gradients grow with the stress components accorsding to the assumed flow law; thus, the velocity components themselves grow with increasing Δx.

Fig. 5. Longitudinal profiles of scaled basal normal deviatonc stress σ (dolled lines), basal shear traction τ (dashed lines) and basal effective stress τ_eff (solid lines) for Δx = 0.1, 0.2, 0.3, 0.4 and 0.5. The profiles for τ and τ_eff are symmetrical and the profiles for σ are anti-symmetrical with respect to the centre of the sliding area.

The scaled driving stress(Whillans, 1987) is τ_d =1 everywhere along the base of the slab. This driving stress is a measure of the total longitudinal force of gravity per unit area acting on a column of ice. Driving stress only depends on the geometry of the ice mass and thus is independent of the fields of shear stress and normal stress. The resistive force acting against basal driving stress is the basal shear traction. However, these forces are not balanced locally but must balance over the entire bed. The average of the basal shear traction over the entire bed equals the average of the basal driving stress over the entire bed, and thus is an invariant that only depends on the geometry of the glacier.

This relation was used to test the accuracy of the numerical solutions for which we found that the average shear traction, which should average to unity, generally matched this value to within 2%. Accordingly, the longitudinal average of the normal deviatoric stress must vanish. This clearly shows the interplay between shear stress and normal stress. Within sliding areas, shear stress is reduced; however, the normal stress transfers the driving forces to the transition zone between sliding and non-sliding. Just outside the sliding areas, additional large shear stress is built up. This demonstrates the pulling and pushing strength of the interior layers that enables the sliding velocity in one area to respond to the presence of neighbouring non-sliding areas and to other more remote sliding areas.

Fig. 6. Ratio of the difference between maximum surface velocity and asymptotic surface velocity In the maximum basal velocity as a fonction of the size of the sliding area L for various area fartants ?.

The fraction of sliding velocity that is transferred to the surface also depends on the spot size. We define this fraction by

(7)

where f/|> and « are the sliding and surface velocity in the centre of the sliding area, respectively, and u^ is the asymptotic surface velocity (=1).This fraction takes values f_s= 0.13,0.36 and 0.69 For ΔX = 0.1, 0.25 and 0.5, respectively, and f_s approaches unity for large Δx. This result reminds us to be cautious in interpreting observed seasonal variations in the surface velocity of glaciers. The difference between lowest and peak velocities only indicates changes in basal velocity somewhere nearby. The sliding does not necessarily occur below the position of the observation and the changes in sliding velocity may be as large as or even larger than the observed peak velocity itself.

3.3. Spatially periodic series of sliding spots

In a second series of numerical experiments, we assume a spatially periodic pattern of sliding and non-sliding zones. In the sliding parts, basal shear traction is set to zero at all corresponding gridpoints to yield a totally uncoupled bed. The pattern is defined by the length of the sliding zone L and the areal fraction of sliding ?, i.e. the ratio of the length of the sliding zone to the total length of the periodically repeating region.

In these experiments, the extent of the sliding area was always taken as ten gridpoints, with an areal fraction of sliding ? = 10/N, where N is the total number of gridpoints in one period. The length of the sliding zone thus becomes L= 10Δx. The sliding patterns are presented for the domain in the parameter space {(?, L) | 0.1 # ? # 0.67 and 1.0 # L # 5.0}.

The ratio between the difference of maximum surface velocity (u_s.max and minimum surface velocity u _s.min within a given spatially periodic region and the maximum basal velocity u_b,max is denoted by

(8)

it represents the maximum sliding velocity expressed as a fraction of the variation of the surface velocity. Figure 6 depicts the dependence of f _p on ? and L. The strength of the local sliding signal clearly fades as areal fractions ? and the sliding area L become smaller. The bridging becomes stronger with smaller scale length of the basal patchiness, and one can again anticipate that the information on the distribution of basal conditions is substantially diffused and effectively lost at the ice surface.

Fig. 7. Basal velocity in the centre of the sliding area as a function of the uncoupling factor 1 — q for an area fraction ? = 0.2 and length of the sliding area L. = 2.5 (dotted line), ? = 0.4,L = 2.5 (solid line), ? = 0.2 and L =1.0 (dashed line).

Fig. 8. Basal velocity in the centre of the sliding area as a function of the size of the sliding area for the spatially periodic pattern and various area fractions ?.

In one scries of numerical experiments, the coupling factor q (Equation (6)) between ice and bed in the sliding area was varied between zero (total uncoupling) and unity (non-sliding). The consequent sliding-vclocity variation is closely proportional to 1-q for all computed cases (Fig. 7). This result provides an additional justification for limiting the investigation of sliding patterns to the extreme case of total uncoupling, q = 0, since all other cases are intermediate between the total uncoupling and total non-sliding conditions.

Figures 8-10 present results for the case of total uncoupling and give the sliding and stress patterns as functions of L and ?. Figure 8 shows the dependence of the sliding velocity in the center of the sliding area on the size of the sliding spot for a sample of areal fractious ?. As expected, the sliding velocity increases with increasing size of the affected area. With smaller ?, the distance between adjacent sliding areas grows and thus the bridging between sliding spots becomes weaker and sliding velocities are reduced.

figures 9 and 10 explore die basal stress concentration in the transition area. The square of the effective stress (see Fig. 9) is a measure of the strain softening of the ice, and additionally indicates the possibility for fracturing (e.g. by assuming a failure criterion). For n = 3, the fourth power of the basal shear traction (Fig. 10) is indicative of the order of magnitude of basal strain heating. These stress values increase with growing areal fraction ?, This is consistent with the previously discussed invariance of the longitudinal average of the basal shear traction. In a periodic pattern, this average is unity over one period. With increasing ?, the area available for the compensation of the zero traction in the sliding part becomes smaller and therefore must support a larger traction than for small ?. The weak dependence of the stresses on L again reflects the above-mentioned approximate similarity rule.

3.4. Haut Glacier d'Arolla

A numerical experiment was conducted with the two-dimensional (longitudinal section) geometry of Haut Glacier d'Arolla, Switzerland (Fig. 11).The first run assumed no sliding throughout the length of the profile. This yields the basal shear traction for full basal coupling. In the second run, some uncoupling at a given section in the lower part of the profile was introduced. The local non-sliding basal shear traction was set to zero to simulate perfect sliding. The stress and velocity fields were then computed by assuming the mixed basal boundary condition, zero basal velocity in the non-sliding parts and the reduced shear traction in the sliding zone.

Fig. 9. Maximum of the square of effective stress in the transition zone between the sliding and non-sliding area as a function of the size of the sliding area for the spatially periodic pattern and various area fractions ?.

Fig. 10. Maximum of fourth power of basal shear stress in the transition zone between the sliding and non-sliding area as a function of the size of the sliding area for the spatially periodic pattern and various areafractions ?.

The stress and velocity fields in the realistic glacier geometry show the same patterns as in the slab solutions (figs 12 and 13). From the result, it is again clear that a locally observed variation in surface velocity only reflects the fact that the sliding velocity has changed somewhere along the glacier bed. The observed change of the surface velocity does not straightforwardly indicate the position or the change in the velocity of sliding.

It is noteworthy that the average t of the basal shear traction over the entire length of the glacier profile is almost the same for both cases; t=1.2365 for the non-sliding case and t_d=1.2352 for the case with the sliding area. The difference lies in the round-off error of the computation due to the specified stopping criterion for the iteration. This invariant, denoted t_d quantifies the total driving force of gravity on the ice mass and thus must also correspond to the average over the bed of the so-called driving stress t_d = pghS' (Whillans, 1987), where p, g, h and S’ are the density of ice, the acceleration of gravity, the local ice thickness and the local surface inclination, respectively. It must be noted that the above observed invariance of the means basal shear traction (averaged over the entire bed) is not proven rigorously in mathematical terms, although it is physically obvious (personal communication from H. Rothlis-berger, 1997). In contrast, the local driving stress deviates substantially from the local basal shear traction, a fact that reflects the intimate interaction between shear stress and normal stress.

Fig. 11. Longitudinal section of Haut Glacier d'Arolla. X denotes the horizontal distance from the top of the profile in metres; h is the altitude above sea level.

Fig. 12. Horizontal velocity component in a longitudinal section of Haut Glacier d'Arolla from the top (left side) to the tongue of the glacier (right side). The line labelled ns shows the surface-velocity distribution for the non-sliding case, the line labelled sl shows the surface velocity for a situation with a sliding zone where total uncoupling is prescribed (see Fig. 13) and the line labelled sb shows the corresponding basal velocity.

Fig. 13. Basal shear traction in the longitudinal section of Haut Glacier d'Arolla from the top (left side) to the tongue of the glacier (right side). The line labelled ns shows the shear traction for the non-sliding case, the line labelled sl shows the shear traction for a situation with a sliding zone where uncoupling is prescribed as shown by this figure and the dotted line shows the driving stress. The surface inclination corresponds to an average over 50m (two gridcells) in all calculations.

Fig. 14. Effect of grid resolution on the calculated longitudinal average of the basal shear stress (solid line) and of the basal driving stress (dashed line) for the non-sliding case and the profile of Haut Glacier d'Arolla. The grid size Δx is given in metres. The dotted line indicates the exact (i.e. computed for very small Δx) average of the basal driving stress.

The average driving stress can be computed from the given geometry alone, independently of any numerical solution. For this reason, the longitudinal average of the basal shear traction yields an accuracy test for the numerical solutions of the higher-order equations. Figure 14 shows the average of the basal shear traction t for the non-sliding situation as a function of the grid resolution Δx, together with the corresponding average of the driving stress t_d. For the specific geometry of Haut Glacier d'Arolla, t_ddecreases almost linearly with increasing Δx. The model solution, on the other hand, produces a more constant t over the same range of Δx, and this t corresponds within 2% to the t _d for very high resolution, i.e. for very small Δx.

3.5. Three-dimensional effects

The three-dimensional dimensionless force-balance equations for a parallel-sided slab are (Blatter, 1995):

(9)

(10)

With the x axis aigned with the steepest surface sope, and the stress strainrate relations are

(11)

(12)

(13)

(14)

(15)

where

(16)

constitutes a prescribed flow law for the fluid under consideration. The boundary conditions at the free surface are t_xz(S)=t_yz(S)=0; the basal boundary conditions in the non-sliding parts are û(B) = û(B) = 0. In the following numerical results, vanishing basal shear traction is prescribed within the sliding parts and the sliding velocity emerges from the model computation.

To study the influence of side drag, three-dimensional model computations were carried out. The multiple-shooting Newton iteration scheme in three dimensions requires large memory and long computation time. Therefore, with the presently available code only about 20 x 20 gridpoints can be handled with reasonable performance on workstations. Better linear algebra, adjusted to the specific form of the problem, and a faster way to solve the large linear systems of the Newton iteration scheme could substantially reduce hardware requirements and computation time. Such investigations are currently under way (paper in preparation by J. Colingé).

The three-dimensional slab is “closed cyclically” in the x and y directions, i.e. in the direction of steepest slope and transverse to it. This produces a periodicity in two dimensions similar to the one discussed in section 3.3 for the plane-flow case. Here, we chose 15 gridpoints in the x direction and 15 gridpoints in the y direction, with a grid size of Δx = Δy = 1.0. The length (X direction) of the rectangular sliding area was taken as 10 gridpoints, and the transverse size was varied between all 15 gridpoints and three gridpoints. The first case corresponds exactly to the plane-strain situation, with areal fraction ? = 0.67, and the three-dimensional computation produced almost exactly the same result as the plane-slab model. Figure 15 shows the X components of surface and basal velocity along a longitudinal (x direction) profile through the middle of the sliding area for different transverse extents of the sliding rectangle. As expected, with decreasing width of the sliding area, the influence of side drag increases and the resulting sliding velocity decreases.

A few computations for infinitely long channel-like sliding areas were performed to investigate the influence of side drag alone. The length of the sliding area was taken as that of the full 15 point x domain and the sliding width was varied from five to nine gridpoints. For these cases, the sliding velocities grow very large (comparable to surging velocities), if the width of the channel exceeds roughly five times the ice thickness (Fig. 16). For ice streams that are wide compared to ice thickness, lateral drag exerts a negligible restraint on sliding velocity. Other limiting factors are a generally larger basal drag and isolated sticky (i.e. non-sliding) spots within a mostly uncoupled bed.

The areal average of the basal driving stress in the x direction is t_b= 1.0 (see also section 3.3) and vanishes in the y direction. Thus, the areal average of t_xz,av=1.0 must be the same for different sliding patterns and correspondingly t_yz.av = 0. In the numerical experiments with a coarse 15x15 horizontal grid and 12 layers in the vertical, this average deviated by as much as 10% from unity.This inaccuracy stems from the gridpoints just up- and downstream of the sliding areas, where stress peaks are partly truncated in the numerical computations.

4 Discussion

4.1. Spatial scales and sliding

The sliding velocity can exhibit large spatial and temporal variations. The patchiness of soft (i.e. sliding) spots seems to follow the basal drainage system (Harbor and others, 1997) and its scale length varies from much smaller than ice thickness to several times the ice thickness. This natural unit of one ice thickness provides a convenient measure for exploring the bridging effect due to normal deviatoric stress. The results for spatially periodic sliding/non-sliding conditions in the slab indicate that a distance of 5-10 times the ice thicknesses is necessary to uncouple the average sliding velocity over one period from that of adjacent periods. For typical valley glaciers, 5-10 times the average ice thickness is comparable to the transverse extent of the glacier tongue as well as to a considerable fraction of the tolal length. Such distances are significantly larger than the observed scale of basal patchiness; thus, at these scale lengths, a sliding area is responding to conditions beyond the immediate neighborhood of the sliding region.

Bridging affects sliding in several ways. (1) It limits the sliding velocity in soft spots by holding the ice between the non-sliding bridge piers. (2) It makes the sliding velocity strongly dependent on the size of the soft spot. (3) It reduces or eliminates the influence of bridged hard spots, effectively increasing the size and importance of neighboring soft spots. As a result, sliding in valley glaciers is a demonstrably non-local process in which the mechanics of bridging plays a central role. This conclusion also holds in the case of realistic three-dimensional ice masses that may have complex topography, complex basal coupling patterns and spatially variable rheology introduced by varying temperature or water content. Ice streams seem to have large sliding areas that are locally interrupted by isolated sticky spots (MacAyeal, 1992), whereas valley glaciers may have a predominantly rough and resistive bed that is locally decoupled along sub-glacial drainage pathways (personal communication from M. Sharp, 1996).

Fig. 15. Sliding and surface velocities in a three-dimensional parallel-sided slab with maximum inclination in the direction of the x axis. Shown are longitudinal profiles of surface and basal velocity components in the x direction for longitudinal sliding fraction ?_l=0.67 and various transverse widths of the sliding rectangles. Horizontal grid size is Δx = Δy = 1 in units of slab thickness.

The spatial variability of sliding and the corresponding basal stress field has consequences for the interpretation of field observations and for the experimental determination of friction laws, even if the sliding velocity is measured directly. The major difficulty is the impossibility of measuring stress components directly. Indirect determination of the basal shear traction requires accurate measurement of near-basal strain rates and reliable knowledge of the form of the flow law and corresponding rheological parameters. The often-used basal driving stress is a poor approximation to the basal drag and thus a poor surrogate for basal shear traction. This is especially true in areas where sliding conditions yield a shear traction that decreases with increasing sliding velocity, while driving stress remains unaffected.

Fig. 16. Sliding and surface velocities in a three-dimensional parallel-sided slab with maximum inclination in the direction of the X axis. Shown are transverse profiles of surface and basal velocity components in the x direction for a longitudinal infinite sliding channel and various transverse widths of the channel. Horizontal grid size is Δx = Δy — 1 in units of slab thickness.

Furthermore, the validity of proposed flow laws and their usefulness for indirect determination of stress components must be questioned. Basal ice that flows through a region of patchy conditions with small spatial scales undergoes rapid changes in stress; the associated ice-deformation changes on time-scales of days to hours. The stress changes are mainly associated with the large gradient of the effective stress, which leads to large stress rates (rate of change of effective stress) as well as to changes in the directions of the principal stresses which rotate at the transition from non-sliding to sliding and vice versa.

4.2. Sliding in higher-order models

Field experiments that investigate basal sliding and friction laws need to be complemented by higher-order models of glacier mechanics. Information on the spatial distribution of sliding conditions is a necessary model Input. In return, the model results may then supply an approximate stress field and can thus help lo constrain friction-law parameters.

The stress and velocity fields in glaciers are determined by the force equations, mass-continuity equation and the stress strain-rate relations (e.g. Hutter, 1993). The boundary conditions are vanishing shear traction at the free surface and a combination of the velocity and/or shear traction along the glacier bed. For non-sliding conditions, the basal shear traction is unknown and must be determined iteratively. Mixed basal boundary conditions may be prescribed, such as basal (or zero) velocity at some parts of the bed and basal shear traction at the rest of the bed. Iterative solution of the equations then yields the basal shear traction and basal velocity in the respective other parts of the bed. From these considerations, we conclude that for numerical modelling the basal shear traction is the appropriate measure of the basal “drag force”.

Wc now investigate a method for estimating basal velocity from the basal shear traction that can be supported by the bed. Our approach is to introduce a local stress reduction at the bed, either caused by decoupling between ice and bed or by dragging a soft deformable layer of material beneath the ice, or both, thus reducing the local basal shear traction t_b.soft to a fraction q of the traction t_b.hard in the non-sliding case

(17)

where 0 < q < 1. Once the basal boundary conditions are given (zero velocity in the non-sliding parts and reduced basal shear traction in the sliding areas) the stress and velocity fields are well defined. The sliding velocity can then be computed with a mechanical model which, among other things, calculates deviatoric stress gradients.

In a numerical implementation of this concept, the coupling factor q characterizes the mechanical coupling for a gridcell centred at the gridpoint under consideration. This concept is certainly problematic, as all concepts of sliding laws are difficult to constrain, and its applicability is restricted to glacier beds that are smooth on scale lengths of the grid resolution. For complete uncoupling, q= 0 and the upper limit on sliding velocity for given a basal patchiness of soft and hard spots can he estimated. This provides a realistic control on the range of plausible basal velocities for given basal patchiness -likely a better one than that obtained by applying a local friction law with poorly constrained coefficients and exponents.

The nature of basal stress reduction and the parameterization of the relevant physical processes lies beyond the scope of this paper. It can be anticipated that q is a function of water pressure. Other factors that may influence q include basal roughness and stress-enhanced flow around obstacles Weertman, 1957), which may make _q dependent on the sliding velocity itself, variable thermal conditions, e.g. patchy cold and temperate sole conditions, and the rheologieal properties of the substrate. We cannot expect that the coupling factor q is solely determined by local basal conditions.

4.3. Basal stress and transition to surge

The invariant average of the basal shear traction (see sections 3.3 and 3.4) constrains the boundary conditions and indicates that the completely non-sliding case provides a unique reference. The non-sliding shear traction gives the lower limit of the traction that must be supported by the bed to maintain non-sliding. If a given region of the bed cannot support this minimum traction, then sliding must occur at this position, accompanied by a local reduction of the shear traction and a corresponding increase in the traction in its neighbourhood.

Sliding can nevertheless occur in places where the shear traction exceeds the above-defined reference. In an area of enhanced shear stress produced by a nearby sliding zone. the enhanced shear stress acts as new reference for non-sliding; if the bed cannot support this enhanced shear stress, then sliding will occur. This makes it necessary to introduce a threshold which determines whether an enhanced shear traction can be supported by the bed without giving way to sliding. If not, sliding must occur at this position. Whether, and to what extent, the shear traction is then reduced at the onset of sliding at this place depends on the type of bed and the corresponding friction law. We can anticipate thai this slip condition, under enhanced stress in the transition zones between initially sliding and non-sliding areas, is closely related to the temporal and spatial variability of sliding.

If the neighbouring bed is not able to support this enhanced shear traction, the sliding zone must expand into this region accompanied by a stress reduction. The adjacent zone with stress enhancement then shifts further away. If the bed can nowhere support the stress concentrations, this expansion continues and the ice mass must dynamically destabilize. These considerations open questions about the propagation and expansion of such sliding areas and the conditions necessary for surge instability.