Notation

1. Introduction

Drumlins are mounds of sediments produced by the action of ice sheets and glaciers. They range in size from a few metres to tens of metres high and range in length from tens of metres to kilometres. Typically, drumlins have blunt upstream and occasionally downstream faces. The problem of explaining drumlin formation, which has troubled glacial geologists for over a century, was re-invigorated in the 1980s when it was realized that sediment deformation could explain many of the structures found within drumlins (Reference Boulton, Menzies and RoseBoulton, 1987) and it was hypothesized that, if subglacial sediment were essentially deforming as a viscous fluid, then drumlin formation might be a problem explicable using fluid-dynamical principles (Reference Boulton and HindmarshBoulton and Hindmarsh, 1987). This has important consequences, because the same physics which explain drumlin formation can also explain ice-stream lubrication by deforming sediment, as observed beneath Ice. Stream B (Reference Blankenship, Bentley, Rooney and AlleyBlankenship and others, 1986). The ability to explain drumlin Formation can thus be seen as a key lest for any theory of subglacial sediment transport and, if subglacial sediment deformation is an important sediment transport route, the study of drumlins provides information about the large-scale Flow of sediment, information which cannot be measured in the laboratory for practical reasons (Reference HindmarshHindmarsh. 1997).

Recently, the blunt faces of drumlins have been explained as shocks (Reference HindmarshHindmarsh, 1996, in press a) and a nonlinear kinematiccal theory of drumlin formation has been developed. However, while a certain proportion of drumlins are obviously moulded from pre-existent relief, there is a sufficient proportion of drumlins composed of the same material as surrounding till (subglacial sediment) fields to suggest that drumlin formation may be an instability in a viscous flow of a till sheet. This idea has long been debated amongst glacial geologists. In this paper, we examine the stability of a till sheet at wavelengths longer than the depth of the deforming till but shorter than the ice-sheet thickness; this is a typical drumlin wavelength.

In general, the flow of two shearing layers of different viscosities but the same density is stable. In this case, the viscosity of the till layers is affected by the normal traction applied by the ice. Moreover, till discharge does not increase monotonically with till thickness (for a given stress) but decreases, owing to increased effective pressures (Reference HindmarshHindmarsh, in press a). On its own, this does not lead to unstable thickening of the till layer but the interaction of this property with the ice flow might. This is investigated in the present paper. It is found that under certain conditions, but not all conditions, that unstable thickening does occur; sheet flow is not stable. A feature which emerges is that, in general, till bed forms migrate. This is not something which has been extensively discussed in the literature and indeed it is hard to see how it might be inferred from formerly glaciated areas. Rock-cored drumlins might at first sight appear to militate against this idea but, since so many drumlins are not rock-cored, one has to regard rock-cored drumlins as atypical in genesis and therefore in dynamics. It is hard to imagine a viscous theory which does not in general predict migration.

The ice-flow model is a variant of the perturbation approach developed by Reference NyeNye (1969, Reference Nye1970) and Reference KambKamb (1970), called here the NK solution, and, since the ice is in a state of simple shear at large distances from the bed, also bears some relation to the analysis by Reference MorlandMorland (1976a, Reference Morlandb); it is not, however, a gravity-driven flow. Rather than slipping Over a perfectly smooth, fixed bed, it is fixed to a thin deforming fluid, the till, whose stress fields can be computed using lubrication-theory approximations (known in glaciology as the thin-till approximation Reference AlleyAlley, 1989)). The evolution of the till profile can then be computed and it is found that in certain parameter ranges the surface relief grows, although of course the till volume remains constant. This is the principal result of the paper; ice flowing over a thin till sheet can promote unstable amplification of relief. This is a necessary component of a theory of drumlinization if it is believed that not all drumlins are created by moulding of pre-existent relief.

This analysis is based upon the viscous flow laws for till where strain rate is proportional to the shear stress to some-low power and inversely proportional to the effective pressure to some low power (Reference Boulton and HindmarshBoulton and Hindmarsh, 1987). The relevant paper was written in the context of six decades of mechanical research which showed that on the small scale sediment was plastic. At Breiбamerkurjökull, the ice overriding deforming sediment does not build up and then slip; en masse, it does not behave plastically. A viscous law, while undoubtedly failing to represent failure, represents the simplest model compatible with large-scale observations. In the same way as observations of dislocations in ice (an event akin to plastic failure) are never held to militate against the validity of viscous descriptions of ice, there appears some to be no contradiction between laboratory experiments which demonstrate plastic behaviour and the posing of a viscous law which applies on a larger scale (Reference Bahr and RundleBahr and Rundle. 1996; Reference HindmarshHindmarsh, 1997). Recent studies, for example, by Reference KambKamb (1991) and Reference Iverson, Hanson, Hooke and JanssonIverson and others (1995), which question the existence of local viscous behaviour, do not explicitly address the question of what an appropriate large-scale law might be. On the other hand, to be of any use, a large-scale viscous law must be shown to be able to reproduce observed behaviour. This paper shows that the viscous flow of till, when coupled with ice flow, can create an instability which appears to be necessary in the format ion of drumlins (Reference HindmarshHindmarsh, in press a). Actually, it also shows that power-law indices as high as 10 (which some might call plastic behaviour) can create drumlins.

2 Stokes Equations and Boundary Conditions For Shear Over a Mobile Bed

The solution domain is an ice-sheet half-space overlying a deforming till layer of finite thickness. The ice and the till are supposed to be fully coupled, so that velocities are continuous across the interface. A shear stress is applied to the ice at a large distance from the bed, which sets up a simple shear in the ice and causes the till to deform. In the base case, the till surface is flat, meaning that a simple shear exists in the ice and that the till deforms with uniform velocity at its upper surface. Small periodic variations in the till surface profile are then introduced, which cause the simple shear in the ice flow to be modified, according to a variant of the NK perturbation of the Stokes equations. Flow in the till is still computed according to lubrication theory principles, implying that the wavelength of the disturbances be greater than the deforming layer thickness. The perturbation to the ice flow implies perturbations to the tangential and normal tractions applied to the till, causing the till flow to be perturbed, and it is the distribution of these which determines whether the till relief grows or decays.

Consider the ice to be occupying a half-space periodic in the horizontal direction in the domain x = [0, 2π/k;] and z ≥ 0, where kis a wave number. We can generate a zeroth-order solution by assuming that the bed is flat and that the till is thus of uniform thickness D =D₀, with the top of thc till at z = 0. Let us suppose that a traction T₀^t is applied to a horizontal surface on the ice at a distance very far from the bed. This sets up a simple shear in the ice and in the till and also creates a horizontal velocity in the ice at the bed z = 0 We find that the zeroth-order velocity is

(1)

where ρ_e0 is the interfacial effective pressure and F_u is a function So be discussed later, This shear implies a flux of till given by

(2)

In this analysis, the effective pressure is a free parameter, which in reality is determined by the subglacial hydraulic system. In this study, horizontal scales are assumed to be sufficiently small that pressure gradients are hydrostatic; there are no horizontal water-pressure gradients. The physics of this has been discussed by Reference HindmarshHindmarsh (1996, in press). Following Reference NyeNye (1969, Reference Nye1970) and Reference KambKamb (1970), we perturb the zeroth-order solution of ice flow by introducing some undulations into the bed. There is a difference here compared with the NK solution, because we have a base solution with simple shear (and thus no possible problems with the Stokes paradox (Reference FowlerFowler, 1981)), the situation is in some ways analogous to the gravity-driven flows computed by Reference MorlandMorland (1976a, Reference Morlandb). Moreover, we follow Reference NyeNye (1970) by working in a perturbed coordinate system where z is measured from above the bed. This transform is covered in more detail in the Appendix.

The perturbation expansion of field values in the ice in this transformed coordinate system is

(3a)

(3b)

(3c)

with basal boundary perturbation expressions

(4)

and it is shown in the Appendix that we obtain the following first-order field equations

with boundary relations

where the normal traction convention is compressive-positive. If we define a new function

(5)

we can retrieve the more familiar set of NK field equations

(6a)

(6b)

(6c)

but now with boundary relations

(7a)

(7b)

The first-order basal velocity is computed using a linearization of the basal velocity relationship (1), so we can write

(8)

(9)

where the derivatives are evaluated around the zeroth-order state. We also need to satisfy the first-order kinematical condition

(10)

where ₁ D₁ (by construction zero in the NK formulation) is given by the linearized till-conservation equation

(11)

(12)

Some readers may feel these perturbations look strange as they appear to lack the expected products of zeroth- and first-order terms. In fact, all the zenith-order terms are in the coefficients R,_t.n.D, Q_t,n,D A quick way of verifying the validity of the method is to remember that a linearization is simply a Taylor expansion; the Taylor coefficients here are simply R,_t.n.D, Q_t,n,D.

The velocities and fluxes are computed following Reference AlleyAlley (1989) and Reference HindmarshHindmarsh (in press a). Till flux can arise either From internal deformation within the till or from fill sliding over the base. When considering internal deformation, the strain rate in the till is given by a double power-law rheology (Reference Boulton and HindmarshBoulton and Hindmarsh, 1987)

(13)

where A_d is a rate factor and a and b are parameters. Then, defining

(14)

whence after integrating with respect to z over the interval -D ₀ ≤ z ≤ D - D₀ we find

(15)

and the flux is given by

(16)

Where

(17)

In these formulae

where g is the acceleration due to gravity, ρ_w, ρ_i and ρ_s are the densities of ice, water and sediment grains, and ø is the porosity of the sediment. We use the values 9.81 m s⁻², 1.0 Mg m ⁻³, 0.917 Mg m ⁻³, 2.7 Mg m ⁻³ and 0.2, respectively, in this paper. These formulae are derived in exactly the same way as the equivalent relationships in Reference AlleyAlley (1989) but account for statically induced vertical interfacial effective pressure gradients arising from the density difference between water and ice Reference HindmarshHindmarsh, 1996). Under static assumptions, an increase in elevation causes an increase in effective pressure because water pressure decreases more rapidly than does ice pressure as one ascends. Alleys formulae may be retrieved by setting α = 0.

We now turn our attention to till—bed sliding. If we use the viscous-type sliding law proposed by Reference HindmarshHindmarsh (1996), we find that

(18)

Sliding and deformation do not occur together in this study. We note that the more complicated derivatives are given by

(19a)

(19b)

(19c)

(19d)

A significant point is that we have assumed that the shear stress is constant within the till. At wavelengths comparable with the deforming till thickness, this approximation breaks down, and gradients in the till thickness and the applied normal stress cause additional flow effects. Our analysis, therefore, is only valid for wavelengths between the thickness of the till layer and the thickness of the ice layer. This is the relevant length scale as most drumlins are of this size.

3 Solution

Let us suppose that the bed profile is given by

(20)

where D₁^α(t), D₁^β(t) (£) are scalar functions of t. The solution to the field Equations (6) is

(21a)

(21b)

(21c)

where the coefficients C^α,P^α,C^β,P^β are determined by the basal boundary conditions. At the base, z = 0, we compute the tangential traction and its gradient as follows

(22a)

(22b)

and the normal traction and its gradient as

(23a)

(23b)

while the interfacial velocities are given by

(24a)

(24b)

Substitution of Equations (20), (22a), (22b), (23a), (23b), (24a) and (24b) into Equation (8) and also into Equation (10), where we have eliminated ∂_t D₁ by using Equation (11) and subsequent multiplication of each of these equations by sin kx and separately by cos kx followed by integration over the horizontal domain yields the following system of four equations for C^α,P^α,C^β,P^β:

(25)

Solutions are discussed below.

Superposability is easy to demonstrate by expanding arbitrary (differentiable) bed shapes in a Fourier series. The above equation holds for each bed frequency, and solutions for the stress and velocity fields may be obtained for each wave number which depend only on the bed geometry at that wave number. This also applies to the linearized till thickness evolution equation discussed below, meaning that evolution equations may be obtained separately for each mode.

The linearized till conservation equation is

Using relations (22b) and (23b), we can show that the phase and anti-phase mode-evolution equations are

(26a)

(26b)

This solution represents a sine wave growing or shrinking in amplitude and translating. The rate of growth of the bed perturbation D₁c| is

(27)

and we can use the idea of the growth-rate constant

(28)

This analysis shows that an incipient drumlin will move. The phase ω is given by

and the migration velocity M of a sine wave with wave number k is given -ω/k. It can readily be shown that

and in this case the velocity M of a sinusoidal till wave with wave number k is given by

(29)

We call this the migration velocity. Where there is no ice till coupling, the migration velocity is equal to the kinematic wave velocity Q_d. It can be shown that M = Q_d + (dq/df){_x²f/_x²D), where f is a field variable other than D, such as T_t. This result generalizes to an arbitrary number officiel variables. Thus, we do not expect the migration velocity to be equal to the kinematic wave velocity when Q_t or Q_n are non-zero — where there are gradients in the field variables, the elevation of a kinematic wave is not constant. Owing to the relationship between the migration velocity and the kinematic wave velocity, it is hard to conceive a viscous theory which does not have migrating drumlins.

4 Examples and Applications

Without loss of generality, since we are considering cases of individual Fourier modes, we set D^β = 0, which considerably simplifies the solutions given below. It is straight for-Hindmarsh: Stability of a viscous till sheet coupled with ice flowward to show that the force exerted by tangential traction must sum to zero at this order, and the additional drag, which is a second-order quantity (Reference FowlerFowler, 1981), is obtained from the normal traction

(30)

We find that the solutions to the matrix Equation (25) for all the Q zero, implying a bed of fixed geometry, are

(31)

If all the R are zero (no slip) then C^α, C^β are zero and we find that ρ^α, ^β = 0 since U_b is zero for a fixed bed. Nothing happens at this order of the approximation. If the basal conditions are very slippery (R_t very large) but the other R remain zero then we retrieve the NK solution C^α =0,P^α = 2ŋk²U_bD₁^α. It is not obvious that we should get this but nonetheless we do.

Note that non-zero anti-phase components C^β, ρ ^β only occur when there is a dependence of the basal velocity on the effective pressure. Moreover, when the bed is not mobile, the basic Equation set (25) and expressions (24a) and (24b) ensure that there is no anti-phase structure in the normal velocity at the interface nor in the normal traction.

The full set of equations, with non-zero Q_t,Q_n,Q_D do possess solutions in rationals of polynomials involving no radicals, which are nonetheless too complicated to peruse. Numerical solutions involving direct inversion of Equation (25) are readily obtainable. Our main question concerns the possible existence of unstable Fourier modes which cause bed relief to grow. We now show that they must be associated with the existence of anti-phase structure. In consequence, since anti-phase structure is excited by the dependence of the deforming-bed viscosity on the effective pressure and through this, on the applied normal traction, we see that it is the dependence of deforming sediment rheology on effective pressure which causes unstable growth of relief.

First, we note that the rate of change following the wave is given by dD/dt|_w= _tD+ M _xD. Since, at the maximum, _xD = 0, the rate of change of the maximum point following the maximum is simply d_tD. Consider now the bed-evolution equation

The elevation maximum occurs at x = Π/2k. At this point. relationship (22b) shows that (in the absence of anti-phase structure) that the first term on the right hand size is zero, as is the third term. In general, a zero R_n implies a zero Q_n but even if that were not the case, a more detailed examination shows that for Q_n < 0, the usual ease, the sediment relief will decline. The feature which causes the sediment relief to evolve is a non-zero R_n.

There are a large number of factors which influence the evolution of the bed profile; the zeroth-order stresses T₀and p_e0, the sediment deformation and sliding rate factors, and the various indices. Since it is evident that there is a direct dependence of the perturbation on the basal velocity U_b, instead of varying the rate factors as a parameter, we vary the Velocity U_b, as a parameter and choose the rate factor to give the required velocity.

Approximately 240000 calculations considering internal deformation have been carried out, comprising the direct product of the following seven-dimensional parameter space: A, Such a large exploration of parameter space would not have been possible using numerical solutions of the full equations. The main aim has been to identify the regions of parameter space which encourage amplification of relief. Figure 1 shows the proportion of positive growth rates for each parameter value. All cases for a particular parameter value are considered: for example, if the effective pressure is 0.002, into this bin goes every case where this is true.

Fig. 1. Proportion of cases where the growth rate is positive, plotted using each parameter as the independent variable. Note that a starts just where p_e finishes. Parameter units are indicated on the main plot.

Thus, the proportion of positive growth rates increases with basal velocity, declines with perturbation wavelength, and shows a non-monotone relationship with till thickness, reaching a maximum at till thicknesses around 50 m. This need not be true for any particular parameter value. The proportion decreases with effective pressure and shear stress, increases with shear-stress index and has a non-monotone relationship on the effective pressure index.

Individual cases are not so clear cut. figure 2a shows how growth rate depends on the till thickness and the basal velocity for the indicated values. For low values of till thickness and basal velocity, growth rates are negative, and there is a monotone increasing dependence on basal velocity but a non-monotone dependence on till thickness, but with growth rates remaining positive. Fur a rather more non-linear flow law, figure 2b shows a simpler dependence on till thickness and a different pattern of negative growth rates. figure 2c has basal velocity and wavelength as independent variables. High wavelength means negative growth rates for low velocities. figure 2d shows that high effective pressures can suppress drumlin formation for small till thicknesses and that there is a non-monotone dependence of growth rate upon effective pressure.

Fig. 2. Plots of growth-rate constant D₀c/D₀^c as a function of indicated independent variables, and for various cases of the indices a and b. In cases where they are not, the independent variables, D₀ = 10 m, U_b = 100 m a^−1, , L = 100 m, p_c. = 10⁵ Pa, and = 10^<5 Pa, Growth rates are shaded where positive and white where negative.

In general, these eases suggest that high velocity and greater till thickness encourage drumlin formation. The result that is encouraged at small wavelengths is dependent on the assumption that the thin till approximation holds good. Drumlins presently found in very thin till sequences may therefore be erosional remnants, if cases exist where longitudinal stress gradients within the till suppress drumlin formation. One curious feature is that positive growth-rate constants are very much larger in magnitude than are negative ones and are high enough to suggest that drumlins can form very fast indeed. Calculations invoking till sliding only and covering the same parameter range found no cases where the growth rate was positive. This suggests that drumlins which have formed from till slipping over its bed may also be erosional remnants or may arise from moulding of existent relief.

A fundamental constraint on any linearization technique is that it makes no prediction of what happens when the perturbation becomes large. The small parameter in the NK perturbation is the bed slope and the perturbation becomes invalid when the slope reaches unity. However, when the wavelength is 100 m and the perturbation amplitude 10 m, the error in the Stokes equations is still only 10%, meaning that qualitative predictions may remain correct provided that the original till thickness was sufficiently large that a 10 m amplitude wave remains small.

By the time slopes reach order unity, we anticipate shock formation to have occurred, which is not a process that can occur in the present linearized theory. Λ further factor likely to be of significance is the appearance of higher harmonics when the perturbation of a linear fluid is considered to higher order (Reference GudmundssonGudmundsson, 1997a) and (perhaps separately) when a fluid with a Glen rheology is considered (Reference GudmundssonGudmundsson, 1997b). These analyses show that higher-order and non-linear effects start to play a significant role well before the aspect ratio reaches order unity. These factors could cause relief to be amplified at different wavelengths from that of the perturbation.

5 Conclusions

This paper has been concerned with establishing a fundamental property of ice shearing over sediments deforming according to a double power-law rheology; does it act so as to amplify relief? The answer is that, under certain conditions, infinitesimal perturbations grow. How large they grow and, whether they will continue to grow once shocks start to form, is not answered by this analysis. This analysis is taken as indicating that drumlins could be formed from a till-sheet flow instability, and is consistent with a body of evidence (Reference Boulton, Menzies and RoseBoulton, 1987; Reference HartHart, 1997) that drumlins are some-limes composed of deformation till and are not necessarily formed by the moulding of pre-existent relief.

The drumlin instability therefore represents another aspect of larger-scale behaviour which can be modelled using viscous laws. Framed as they are at the larger scale, there is no reason why they should represent smaller-scale behaviour and it is a misconception to suppose that they must. Furthermore, the viscous-type law must, on the present state of knowledge, be more appropriate for large-scale ice-sheet modelling than the plastic law.