Introduction
Most of the internal deformation of glaciers takes place in a relatively narrow region close to the base. One must, in general, expect small local bedrock undulations, which protrude into the ice, to affect the basal flow and possibly to cause a flow pattern considerably different from the one predicted by the well-known plane-slab solution. Knowledge of the flow perturbations associated with bedrock undulations is, among other things, important for the interpretation of slope measurements in ice, and because of their possible effect on ice stratigraphy.
This paper develops an analytical solution for a highly viscous medium flowing over a perfectly lubricated sinusoidal bed and analyses its properties. A numerical treatment of the problem for high roughness values, where the roughness r is defined as the ratio of the beds amplitude a to its wavelength λ (r : = a/λ) and using Glen’s flow law, is the subject of a further paper (Reference Gudmundsson,Gudmundsson, 1997).
Notation
Dimensional quantities are usually in lower-case letters and non-dimensional quantities are in capital letters.
Previous Work
Theoretical treatment of flow over undulating bed is difficult and only a few analytical solutions exist (Reference Nye,Nye, 1969, Reference Nye,1970; Reference Kamb,Kamb, 1970; Reference Morland,Morland, 1976a, Reference Morland,b; Reference Fowler,Fowler, 1979, Reference Fowler,1981). These solutions often apply to somewhat idealized conditions at the glacier bed but nevertheless give a valuable insight into the nature of the flow. Numerical work has so far been limited to a few cases (Reference Meyssonnier,Meyssonnier, 1983; Reference Schweizer,Schweizer, 1989; Reference Schweizer, and IkenSchweizer and Iken, 1992; Reference Raymond,Raymond, unpublished).
Reference Nye,Nye (1969) and Reference Kamb,Kamb (1970) found an approximate solution for a highly viscous Newtonian fluid sliding over a perfectly lubricated bed. They used a perturbation approach and calculated the flow field to first order in ε : = ak, where the vertical position of the bed line is given by z 0 = a sin kx.
Ignoring the effect of regulation and assuming no tangential traction, the boundary conditions along the bed line are
and
where tan β(x) : = dz 0(x)/dx. The problem is depicted in Figure 1 and the variables have been defined above under Notation.
The field perturbations are to first order (Reference Nye,Nye, 1969; Reference Kamb,Kamb, 1970):
where τ : = √σʹII, which is sometimes called the effective stress. p ∞ is the pressure applied at the upper boundary of the medium, σij are the components of the stress tensor and σʹij are the components of the deviatoric stress tensor. The driving force of the motion is a constant shear stress applied at the upper boundary. The basal sliding velocity, u b, is given by
where η is the viscosity of the ice, and τ b the driving stress.
Expressions (3a) – (3f) display some interesting features. One of them is the fact that τ, in Equation (3f), shows no dependence on x. This will of course also apply to the second invariant of the strain-rate tensor. Another interesting feature of the linear solutions given above is the occurrence of extrusion flow, which is here defined as an increase of the horizontal flow-velocity component with depth. At the point kx = 3π/2 + 2πl where l ∈ N, where N is the set of integers and z = z min : = 1/k, vx has a local minimum:
From z = z min downwards to the bed the horizontal velocity increases. Note that, since ak ≪ 1, it follows that z min ≫ a.
Extrusion flow, a term introduced by Reference Demorest,Demorest (1941, Reference Demorest,1942) has been a subject of some debate in the glaciological literature. On theoretical grounds, it can easily be shown that a global extrusion flow, that is an increase of the horizontal velocity with depth throughout an entire glacier, is impossible since the overlying mass will then experience a force in the main direction of flow, which is not counter-balanced by any other force, leading to an accelerating velocity (Reference Nye,Nye, 1952). There are, on the other hand, claims of extrusion flow having been directly observed by borehole deformation measurements (Reference Hooke,, Holmlund and IversonHooke and others, 1987) and by observations within subglacial caves close to the bedrock interface (Reference Carol,Carol, 1947). Extrusive flow has also been observed within subglacial sediments (Reference Blake,, Clarke and Gérin.Blake and others, 1992). Arguments supporting (global) extrusion flow based on mass-balance measurements have also been given (Reference Streiff-Becker,Streiff-Becker, 1938; Reference Seligman,Seligman, 1947).
Reference Morland,Morland (1976a) derived second-order solutions for flow over bedrock undulations using boundary conditions (1) and (2), and calculated explicit solutions valid along the bed line of a sinusoidal bed. Solutions based on his work, which are valid for the hall-space above the bed, are given below and discussed.
Reference Meyssonnier,Meyssonnier (1983) and Reference Schweizer,Schweizer (1989) did FE calculations of flow over a sinusoidal bed. Meyssonnier obtained, in some of his numerical calculations, a point of maximum relative horizontal velocity that was situated above the peak of the sine wave, and some of Schweizer’s calculations show a point of relative horizontal velocity minimum situated above the trough of the sine wave.
Second-Order Solutions for a Gravity-Driven Flow
Reference Morland,Morland (1976a) incorporated gravity as the driving force of the motion and calculated terms to second order in ε. For the special case of a sinusoidal bed, he gave expressions valid for the pressure field and for the velocity components along the bed line, i.e. at z = 0. Using Morland’s results, one can calculate the velocity and the stress field as functions of x and z, only somewhat laborious work is involved. All equations in this section follow from Reference Morland,Morland (1976a).
The basal sliding velocity u b is
where the controlling wave number k * is defined by Equation (17), and
is defined as Table 1 compares the notation that is used here with the notation of several other authors.
If k/k, ≪ 1, the effects of regelation are negligible and the basal sliding velocity is given by
On the other hand, if k/k,* ≫ 1, which is the pure regelation limit,
u b has a minimum at k = k *. For a given amplitude-to-wavelength ratio the largest part of the drag is contributed by the Fourier components of the bed with wavelengths around λ*, where λ* = 2π/k * is given by Equation (15) (Reference Nye,Nye, 1969). Note that
where
Equation (9) is a useful relation that can be used to eliminate the sliding velocity from the following equations.
The velocity field is given by
and
where A m is defined through Equation (18), and h is the mean glacier thickness.
For the special case of z = 0, Equation (11a) and (11b) reduce to solutions for the velocity field along the bed line given earlier by Reference Morland,Morland (1976a). Note that for z = 0 the first-order term in Equation (11) vanishes.
Equation (11a) and (11b) lead to the following expressions for the strain rates:
and
The second invariant of the strain-rate tensor is then found to be
Equation (17) shows that, to second order. depends on both x and z and not only on z as is the case in the first-order Nye/Kamb solution.
Finally, the pressure distribution is given by
These expressions can be used to calculate the flow and the sliding velocity for a general bed geometry as long as ε ≪ 1.
The Effect of Regelation on the Flow Field
Regelation is only important at wavelengths comparable to or smaller than the transition wavelength λ* (Reference Weertman,Weertman, 1957, Reference Weertman,1964, Reference Weertman,1979; Reference Nye,Nye, 1969; Reference Kamb,Kamb, 1970). where
and where K 1 and K B are the thermal conductivities of the ice and bed, respectively.
Two parameters and A m) enter the flow solutions (Equation (11a) and (11b) that describe the relative importance of regelation to viscous flow (Reference Morland,Morland, 1976a). is the ratio of the bed wavelength to the transition wavelength, i.e.
where k * is the controlling wave number, given by
In the no-regelation limit β = 1 (β is defined by Equation (10)) and in the pure-regelation limit β 1 = 0. A m is given by (Reference Morland,Morland, 1976a)
The effect of freezing and melting on the flow field is negligible if A m ≪ 1 and that is almost always the case (Reference Morland,Morland, 1976a), which is the reason for ignoring the effect of regelation on the flow field in the following discussion of the properties of Equation (11a) and (11b).
Dimensionless Form of the Flow Solutions
For the following discussion, it is of advantage to rescale the dimensional quantities and to put the equations in a dimensionless form. To this end, dimensionless vertical and horizontal length scales, denoted by capital letters, are defined by
where the wave number k is used as a scaling factor. The velocity field is scaled by the sliding velocity, so that
The dimensionless parameters which enter the problem are the slope parameter ε and the thinness parameter δ : = (kh)−1. Regelation will be ignored so that the following discussion is only valid for λ ≫ λx.
Using the above-defined scalings, the velocity field is
and
where use has been made of Equation (7).
Overall Features of the Horizontal Velocity Field
Before going into a somewhat tedious mathematical discussion of the properties of the flow field, let us look at some contour plots of the horizontal velocity field to get an overall idea of the flow perturbations caused by the sinusoidal bed.
Figure 2 depicts VX (X, Z) as a function of X and Z, for ε = 0.01 and δ = 0, according to Equation (11a). The bed line is flat, since in the mathematical solution the sinusoidal bed profile has been projected on to the line Z = 0. Note that and correspond, respectively, to the peak and the trough of the sinusoid. The most conspicuous features of the figure are the stationary points situated above the peak and the trough of the sinusoidal curve at Z ≈ 1.
In Figure 3, the value of ε has been changed from ε = 0.01 to ε = 0.1 as compared to Figure 2. The effect of this increase in ε is to move the local maximum of the horizontal velocity field upward away from the bed, and the minimum of the velocity field down towards the bed. The amplitudes of the velocity perturbations are also considerably larger. Furthermore, a saddle point, where there is a local maximum in horizontal direction but a local minimum in vertical direction, can be seen above the velocity maximum at .
Increasing the value of ε even further, as has been done in Figure 4, where brings the minimum towards the bed line. The maximum point and the saddle point at have disappeared. As will be shown below, the maximum point and the saddle point cancel each other for δ = 0 at ϵ ≈ 0.138.
The effect of changing the value of δ somewhat on the horizontal velocity held can seen by comparing Figure 5, where ε = 0.1 and δ = 0.1, to Figure 3, where ε = 0.1 but δ = 0.0. The velocity maximum moves slightly towards the bed line and the saddle point further away from the bed as δ is increased from 0 to 0.1. Increase in δ causes the local minimum of the horizontal velocity field to move away from the bed line.
One of the interesting features of Figures 2, 3 and 5 is that above the velocity maximum and below the velocity minimum a region where the horizontal velocity increases with depth (extrusion flow) is found. In the next section, the exact conditions under which extrusion flow develops are determined. Since the following discussion is somewhat tedious, the reader who is not interested in the line details of the matter may find it better to skip the next section and read the summary of the results given in the last section of the paper.
Extrusion Flow
It is of particular interest to know when extrusion flow occurs according to Equation (11a). This question can be answered by investigating when VX (X, Z) has a local maximum or minimum for Z > 0. A necessary criterion for a stationary point of VX , Z) is that ∇VX (X, Z) = 0.
The horizontal positions of the stationary points
Differentiation of Equation (11a) gives
and shows that ∂VX (X, Z)/∂X = 0 has as solutions X = π/2 and X = 3π/2 with no restrictions on Z. (It is to be Understood that because of the periodicity of the sine function, an integer multiple of 2π can always be added to the values of X, although it will not be explicitly so written.)
There is another interesting set of solutions given by
It can be shown that this solution set is a second-order effect which is only important in the immediate vicinity of the bed line. This solution branch will not be discussed here further. The interested reader can find a detailed analysis of this solution set out in Reference Gudmundsson,Gudmundsson (1994a). Ignoring this second-order effect there are, hence, only two stationary points. One is located above the peak of the sinusoidal curve, at X = π/2, and can be shown to be a point of maximum velocity with respect to X. The other stationary point is situated above the trough, at X = 3π/2, and is a point of minimum horizontal velocity with respect to X.
The Vertical Positions of the Stationary Points
Differentiating V X(Х, Ζ) with respect to Z and setting the resulting expression equal to zero gives
The interesting cases to be considered are X = π/2 and X = 3π/2, but there are also solutions at X = 0 and X = π with Z = 0. These two solution points, situated at the bedrock interface, are saddle points, where the horizontal velocity obtains a maximum with respect to X but a minimum with respect to Z. The existence of these points is a second-order effect (Reference Gudmundsson,Gudmundsson, 1994b) and they will not be discussed further.
Stationary Points Situated above the Peak of a Sinusoid
Let us begin with the case X = π/2 in Equation (24) and see if there is a solution to the resulting equation
This is a non-linear equation that does not have a solution in a closed form. By plotting the lefthand side (L(π/2, Ζ, δ)) and the righthand side (R(π/2,Z,ε)) separately, as is done in Figure 6, one sees that there will be a solution to Equation (25) if the righthand side, for at least one value of Ζ, becomes greater than, or equal to, 1 − δZ. This will happen if ε is less than some particular value, that will now be called ε crit π/2. There will, hence, only be stationary points at X = π/2 for some δ if ε < ε crit π/2. To determine ε crit π/2, write ε in Equation (25) as a function of Z and δ:
Then ε crit π/2 can be found by maximizing ε(Ζ, δ) given by Equation (26) subjected to Z > 0 and 0 ≤ δ ≪ 1. By solving ∂ε(Ζ, δ)/δZ = 0 numerically, ε crit π/2 can be calculated as a function of δ. The details of these calculations can be found in Reference Gudmundsson,Gudmundsson (1994b).
It is found that there are in general two stationary points along the line X = π/2. An inspection of the determinant of the Hessian matrix of VX (X, Z) at the point X = π/2 shows that one of the stationary points is a point of relative maximum. This velocity maximum will now be called U max π/2. The other stationary point is a saddle point, where VX (X, Z) has a maximum with respect to X, but a minimum with respect to Z. The horizontal velocity component at this point will be called Usaddle π/2.
ε crit π/2 is depicted in Figure 7 as a function of δ. For δ = 0, U max π/2 and U saddle π/2 exist as long as ε < ε crit π/2 ≈ 0.138. The figure also shows dial the value of ε, above which the horizontal velocity has no stationary points at X = π/2, increases with δ.
The changes in the vertical coordinates of both U max π/2 and U saddle π/2 as ε increases from zero are interesting. For ε = 0, U max π/2 is at Z = 1 and saddle π/2 at Z = 1/δ. As ε increases from zero towards ε = ε crit π/2, U max π/2 moves upwards away from the bed, while U saddle π/2 moves downward towards the bed line. Both stationary points then unite at Z crit ≈ 1.98 for ε = εcrit π/2, and disappear (Reference Gudmundsson,Gudmundsson, 1994b). These changes in the vertical coordinates of U max π/2 and U saddle π/2 are illustrated in Figure 8.
The stationary points U max π/2 and U saddle π/2 are not the only possible stationary points above the peak of the sinusoid which the horizontal velocity field can have. Figure 6 shows that the curvature of R(π/2, Ζ, ε) curves (dotted lines) can change, in which case the L(π/2, Z, δ) curves (solid lines) can cross the dotted ones not only twice but three times. In this case, a third stationary point, U Up π/2, will be found along the line X = π/2 together with U max π/2 and U saddle π/2. It can be shown (Reference Gudmundsson,Gudmundsson, 1994b) that the vertical position of U Up π/2 will always be above Z = 3 and that it will in general be found close to the surface.This stationary point may appear because of the assumption, made in the derivation of the flow solutions, that the surface remains flat at all times. This assumption will only be approximately true when the bedrock is undulated. The properties of this stationary point will, hence, not be discussed further.
Stationary Points Situated above the Trough of a Sinusoid
The other possible X value, beside X = π/2, for a stationary point of VX (X, Ζ) is X = 3π/2. Inserting X = 3π/2 into Equation (24) leads to
The lefthand side will always be greater than or equal to zero, because Z must always be within the range 0 ≤ Z ≤ kh, The 1/ε − e−Z term on the righthand side will always be positive. There is therefore no solution possible for Z > 1, but for Z ≤ 1 there will always be a solution to Equation (27) as long as This limiting value for ε will be called ε crit 3π/2, so that for δ ≪ 1. For ε → 0, the vertical coordinate of this stationary point will be at Z - 1, and for the point will be situated at Z = 0 independently of the value of δ. By looking at the determinant of the Hessian matrix, one finds that this stationary point is a point of relative minimum, and it will be referred to as U min 3π/2.
Equation (27) can be solved for ε, giving
As δ varies, the position of U min 3π/2 changes somewhat. This is depicted in Figure 9, which shows the position of U min 3π/2 as a function of ε for two different δ values. For a given ε value and δ > 0, there will be two solutions for Z of Equation (28). Only one of these solution points will be situated below Z = 1/δ (the other one is above the surface, and is a mathematical artifact) and it is this solution that is depicted in Figure 9.
Note that the perturbation approach is only valid as long as ε ≪ 1. A value of is, in this respect, rather large. It is therefore not clear whether the prediction that the stationary point above the trough of the sinusoid exists as long as саn be trusted. Numerical approach seems to be the only possibility of getting a definitive answer.
The existence of the stationary points U min π/2, U saddle π/2 and U min π/2 shows that there will be two regions of extrusion flow. One is at X = π/2, which extends over the region that lies between the saddle point (U saddle π/2) and the maximum point (U max π/2), and another one at X = 3π/2 that extends from the bed towards the point of local velocity minimum (U min π/2). The region of extrusion flow along the line X = π/2 will only exist if ε < ε crit π/2 (see Fig. 7) and the extrusion flow at X = 3π/2 only if
Discussion
Comparison of the First- and the Second-Order Velocity Solutions
The solution of Nye and of Kamb predicts, as said earlier, regions of extrusion flow. In this solution, the vertical positions of the local maximum of the horizontal velocity field above the peak of the sinusoid. U max π/2, and the local minimum above the trough, U min 3π/2, are always at Z = 1 independent of the value of ε. The solutions of Nye and of Kamb have no saddle point. These findings of the first-order solution are reproduced by the second-order solution in the limit ε → 0, as is to be expected. As ε increases somewhat, there are, however, profound differences between these two solutions. The inclusion of the effects of gravity on the flow field limits the occurrence of extrusion flow to a certain range of ε values, and causes the vertical positions of the stationary points of the horizontal velocity field to depend on the amplitude-to-wavelength ratio of the bedrock undulations.
Implications for Borehole Deformation
The influence of the wavy nature of the bed on the flow is not largest at the bed but at some distance above it. Typically, this distance will be about 1 (or z = l/k). Often the deformation of a borehole is used to get information on the rheological behaviour of the ice. One must interpret carefully the data from the lowest part of a borehole, since the exact form of the bedrock, which is usually not known, can have a large effect on the flow. Figuire 3, for example, shows that the perturbed flow (the second and the third terms of Equation (21a)) dominates the gravity-driven plane flow (the first term of Equation (21a)) in the region kz < 3.
The effect of extrusion flow on the deformation profile with depth will be to reverse it with respect to what one would expect from a simple plane-slab flow. Although the discussion here has been limited to one particular type of boundary condition (free-slip), extrusion flow can be expected to occur under other types of boundary conditions as well. It will, in general, be the lowest section of the bore-hole, within the vertical distance λ/2π. where λ is the wavelength of a typical bedrock undulation, which will be affected by extrusion flow.
A Physical Explanation for Extrusion Flow
Upon reflection, it becomes evident that the vertical contraction and expansion of the ice close to the bed is responsible for the extrusion flow. At some distance sufficiently far above the bed, let us say at z = z 1, the ice moves parallel to the mean bed slope. For kx = π/2, a high-pressure zone develops above the bed, which causes a Poiseuille flow, superimposed on the gravity-driven plane flow (GDPF) solution. The maximum of the Poiseuille velocity profile is at (z 1 − z 0)/2 and, if its decrease above that point is faster than the increase of the GDPF velocity profile, a velocity maximum will be found. Since the influence of the bed profile on the velocity field is (because of the factor e−kz) limited to a zone of height proportional to 1/k, z 1 will be proportional to 1/k and one will expect the position of this maximum also to be proportional to 1/k. As a a matter of fact Equation (3a) has a maximum at z = 1/k for kx = π/2. At kx = 3π/2, the ice is expanded vertically and the Poiseuille flow profile reverses, causing a velocity minimum, again at z = 1/k.
The Nye and Kamb solution ignores GDPF and actually expresses nothing but this contraction and expansion due to the wavy nature of the bed. The velocity minimum and maximum therefore never disappears no matter how ε is varied. If, on the other hand, the GDPF is present, it induces a subtle interplay between the Poiseuille flow and the GDPF. Therefore, only certain ε values give rise to this interesting flow behaviour.
Higher Harmonics
An interesting feature of Equation (11a) and (11b) is the presence of the first harmonic of the fundamental period. A simple physical argument shows that one must expect frequencies other than sin kx to appear in the solution, unless ε ≪ 1, and that they will become more pronounced as ε becomes larger.
Let us suppose that there were no higher harmonics in the expressions for vx , and vz , for all values of ε, so that
where ĉ 0 and ĉ 1 are some unknown constants. It follows that
On the other hand, using the exact boundary condition (1) and z 0 = a sin kx one obtains
at the base. Comparing Equation (30) to Equation (31) shows that ĉ 1 must be equal to u b ak and that ĉ 0 has to be zero if expression (29) is to be true. But, this must be true for all ε values, because Equation (31) is always valid. On physical grounds. ĉ 0 = 0 can, however, be rejected; for high ε values ux will certainly not be a constant on z = z 0, which means that the starting-point of Equation (29) must be incorrect. The velocity will therefore, in general, not be a single harmonic, although this may be approximately true for small roughness values.
The Stress Field
Using results from Reference Morland,Morland (1976a), the two-dimensional stress field can be calculated. Detailed results have been given in Reference Gudmundsson,Gudmundsson (1994b). It turns out that, in contrast to Equation (3f), the τ does depend on x if second-order corrections are considered, as can in fact also been seen from Equation (13). Another interesting fact is that τ attains its largest value at kz = 1 independently of ε.
Conclusions
Basal flow has been analysed using analytical solutions for a two-dimensional flow over bedrock undulations and criteria for extrusion flow have been given. Except for some second-order effects, the horizontal velocity field can have at most three different stationary points close to a sinusoidal bed line:
If , there is a minimum point (U min 3π/2) situated above the trough of the sinusoid (at X = 3π/2). As ε increases from zero towards . the vertical position of U min 3π/2 moves from Z = 1 towards the bed line. As is reached, U min 3π/2 hits the bed line and disappears. Since VX (X, Z) increases with depth below U min 3π/2, a zone of local extrusion flow is found close to the bed line as long as U min 3π/2 exists.
Above the peak of the sinusoid (at X = π/2), there is a point of maximum surface-parallel velocity (U max π/2), at which vertical position Z is always within the bounds of 1 < Z < Z crit. It is found that limε→0 Z crit ≈ 1.98 for δ = 0. U max π/2 moves upwards away from the bed line with increasing ε.
Also, above the peak of the sinusoid (at X = π/2), a saddle point of (U saddle π/2 of VX (X, Z) develops, which exists for exactly the same range of ε values as U max π/2. U saddle π/2 is characterized by a maximum with respect to X and a minimum with respect to Z of VX (X, Z). Between the vertical positions ofU saddle π/2 and U max π/2, VX increases with depth.
It must be stressed that the extrusion flow discussed here, which is caused by bedrock undulations, is of local character. It should therefore possibly be called local extrusion flow, in order to distinguish it from extrusion flow encountered earlier in the glaciological literature, which was of global nature in the sense that it extended over a large area. The local extrusion flow does not add to the ice flux as global extrusion flow was thought to do. Local extrusion flow close to bedrock undulations is understandable in simple physical terms and must be expected to be a general feature of basal flow.
Acknowledgements
This research was supported by Swiss Science Foundation grant 20-29619.90. I thank A. Ike, C. F. Raymond, K. Hutter and D. Bahr for thorough reviews of the manuscript, and B. Hallet for bringing Η. Carol’s paper on the formation of roches moutonnées to my attention.