## Abstract

The Hindmarsh instability theory of drumlin formation is applied to the study of interfacial instabilities, which may arise when ice flows viscously over deformable sediments. Here, the analytic form of this theory is extended to the case where the ice is Newtonian viscous and of finite depth, and where the basal till can be both sheared by the ice and squeezed by basal effective pressure gradients: previous authors assumed infinitely deep ice, based on the assumption that the developing waveforms had wavelength much less than ice depth.

The previous infinite depth theory only allowed transverse instabilities to occur, and these have been associated with the formation of ribbed moraine; one of the purposes of extending the analysis to finite depth is to see whether three-dimensional instabilities, which might be associated with the formation of drumlins or mega-scale glacial lineations, can occur: we find that they do not. A second purpose is to calculate under what circumstances the infinite depth theory provides accurate prediction of bedform development in ice of finite depth *d*_{i}. We find that this is the case if the waveforms have a wavelength less than approximately 1.2*d*_{i}. Finally, the finite depth theory allows us to compute, for the first time, the response of the ice surface to the developing unstable bedforms. We find that this response is rapid, and we give explicit recipes for the surface perturbation transfer functions in terms of the perturbations to the basal stress and the basal topography.

## 1. Introduction

Drumlins are streamlined hills that are usually found in clusters, and there is a scientific consensus that they have formed through the action of formerly existing ice sheets. They have been observed and studied in the literature for well over a hundred years (Kinahan & Close 1872), and a good deal of discussion as to their origin has occurred over that time (e.g. Gravenor 1953; Menzies 1979). More recently, there has been a separate emphasis on ribbed moraine (Dunlop & Clark 2006) and mega-scale glacial lineations (MSGLs), particularly as the latter have recently been observed to be actively forming under the Rutford ice stream in Antarctica (King *et al.* 2009). These three types of bedforms represent a continuum from two-dimensional transverse ridges (ribbed moraine) through drumlins to MSGLs, which are longitudinally oriented rolls aligned with the flow direction (Sugden & John 1976).

In the classical literature, the genesis of drumlins is often discussed in the light of concepts of deposition and erosion. Deposition involves the idea that subglacial till is transported by the ice and then deposited, while erosion refers to the idea that pre-existing sediments are eroded by the ice, or by subglacial water. The latter idea may be more readily consistent with some observations of stratified drumlin interiors (for example, in the Puget Sound), although others may be constructional, as in northern Saskatchewan (Shaw 1996); the former may be consistent with drumlins formed largely of till, though some of these may also be erosional (e.g. Kor *et al.* 1991; Sharpe *et al.* 2004).

Although erosion and deposition are no doubt important agents in drumlin growth, neither concept has been shown to produce the landforms that are observed, based on a quantitative theoretical model. More recently, a formation theory based on the concept of an instability of the ice–sediment interface has been proposed (Hindmarsh 1998*a*), which shows that if the ice is underlain by wet saturated sediments at high water pressure, then instability will almost inevitably occur, under very general assumptions concerning basal till deformation and ice sliding. The principal ingredients of this theory are slow Newtonian viscous ice flow and water pressure-mediated sediment shear flow. While water is essentially present in this instability theory, subglacial stream flow has so far been largely ignored, as has possible sediment removal by such streams. A number of other formation theories have been advanced in recent years, particularly for MSGLs. These include a groove-ploughing theory of Clark *et al.* (2003), and a more recent instability theory based on the behaviour of ice as a non-Newtonian fluid (Schoof & Clarke 2008), although neither theory provides a particularly convincing explanation of these landforms.

However, the principal competing mechanism to describe drumlin formation is the so-called flood theory advanced by Shaw and co-workers (e.g. Shaw & Kvill 1984; Shaw & Sharpe 1987; Shaw 1994), which envisages subglacial flood water eroding sediment and thus carving drumlins. This mechanism has been widely criticized (e.g. Clarke *et al.* 2005), perhaps because the popular version invokes enormous floods, which are less easy to understand, even though it is generally accepted that huge floods arising from ice-dammed pro-glacial lakes have occurred in the past (Bretz 1923, 1969; Clarke *et al.* 2004). Insofar as the flood hypothesis involves the transport of subglacial sediments, some elements of it may have a bearing on the instability hypothesis; specifically, this will involve including the transport of sediment in stream flow as well as its transport by ice. Water is essential in the instability hypothesis, and if there is water, then it can be expected to erode sediments. This is because groundwater flow through till is commonly thought to be insufficient to remove the subglacially produced meltwater, which must therefore form subglacial streams. If the stream flow is sufficiently rapid, then sediment transport will occur. The importance of such stream-derived sediment transport is then largely one of degree, as stream flow and till deformation will compete in transporting the sediments.

The instability theory predicts waveforms transverse to ice flow, which can be identified with ribbed (or Rogen) moraine. Furthermore, exploration of the governing parameters has been shown to predict wavelengths of such ridges that match well with extensive observations of the spacing of ribbed moraine ridges (Dunlop *et al.* 2008). While the instability theory does make predictions of amplitude and wavelength, for example (Fowler 2009), it has not yet been shown to produce three-dimensional waveforms such as drumlins. This is due to the fact that the theory does not predict purely longitudinal roll wave instability.^{1} Observations of MSGLs under former and present ice sheets (Clark 1993; King *et al.* 2009) are suggestive of purely longitudinal roll instability, but no such instability was found, assuming the ice is infinitely deep (Schoof 2002*a*). Part of the intention of the present paper is, therefore, to study the instability theory in the case where the ice depth is finite, to see whether that alters the conclusion and indeed to see how much the infinite depth theory is altered quantitatively. A second purpose is to provide a description of the ice surface perturbations produced by the growing bedforms, thus extending earlier work by Gumundsson (2003), which dealt with the surface perturbations produced by prescribed basal topography.

## 2. Mathematical model

A feature of previous theoretical analyses (Fowler 2000; Schoof 2007) is that they assume, for convenience, that the ice is infinitely deep. Formally, this corresponds to the assumption that the corrugation parameter *σ*=*l*/*d*_{i}, where *l* is the drumlin length scale and *d*_{i} is the ice depth, is small. Drumlin lengths are mostly in the range of 200–1000 m (Clark *et al.* 2009), while typical ice depths where they form may range from 500–4000 m, this being the present day range of ice depths in ice caps and ice sheets. It is clearly a possibility that *σ*∼*O*(1) or that *σ*≪1. In fact, the assumption of infinite ice depth allows the drumlin-forming instability to occur, suggesting that the finite depth is not a crucial constituent of the theory. On the other hand, the instability is fundamentally two dimensional, and three-dimensional instabilities do not preferentially occur.

This has been held as a weak point of the theory (Schoof 2007), insofar as drumlins are inherently three dimensional, and MSGLs, which one might suppose represent an extreme drumlinized bedform when the ice velocity is large, are essentially two-dimensional rolls oriented in the flow direction. They cannot appear in the instability theory, at least with infinite ice depth, and this is surprising because there appears to be a perfectly good instability mechanism to form them. This is that where the ice is thicker, the basal stress should be larger, and thus sediment removal will be larger. This provides a positive feedback: where the ice is thicker, sediment removal is more efficient, which further deepens the ice. Such a positive feedback suggests an instability mechanism whereby a small transverse undulation might grow, forming MSGLs.

This instability is essentially analogous to the rill-forming instability studied by Smith & Bretherton (1972), and a clue to its absence in the infinite depth drumlin theory may lie in the fact that the basal stress delivered by the overlying ice flow is not then alterable by the flow over the underlying topography. Since the Smith–Bretherton theory involves essentially the same physical ingredients as we have here (that is to say, flow of a fluid of finite depth over an erodible and transportable substrate), we might suppose that it is the finite depth of the ice flow that will allow three dimensionality in the model. To this end, we pose the model in the general setting of an ice flow of finite depth.

### (a) Effective pressure

We modify the discussion in Fowler (2009) by more specific consideration of the drainage. The instability theory indicates that drumlins can only form where the effective pressure of the drainage system is very low, comparable to the basal shear stress, and we suppose that this is associated with a distributed system of shallow canals, as envisaged by Walder & Fowler (1994). Broadly speaking, we envisage the inter-drumlin space to be swamp-like, with generally low water fluxes. Based on the observations of subglacial floods below Antarctica (Wingham *et al.* 2006), we may expect these low-lying swampy areas to be subject to occasional floods that can flush sediments through the system, somewhat similar to (but much less catastrophic than) the scenario envisaged by Shaw (1983).

In order to describe subglacial drainage, we make the particular assumption that the stream system is isolated from the evolving topography, essentially because the streams are supposed much smaller in width than the landforms, and this allows us to suppose that their drainage leads to a relation between the effective pressure in the stream system, *N*_{c}, and the water flux through the system (which we suppose known). This effective pressure in the drainage system is given, as with other drainage theories (e.g. Röthlisberger 1972; Walder & Fowler 1994), by the difference between far-field (cryostatic) ice pressure, and the water pressure in the streams is determined by the elevation of the subglacial streams, which we may suppose are at the lowest point of the evolving topography; thus, we have
2.1
where *p*_{c} is the water pressure in the local drainage system, *p*_{a} is the atmospheric pressure, *z*_{i} is the ice surface elevation, *ρ*_{i} is the ice density, *g* is the acceleration due to gravity and *s*_{−} is the minimum value of *s*, the ice–till interface elevation. We are specifically assuming that this far-field pressure is ‘far’ from the stream system on the scale of the streams, but ‘near’ the streams on the scale of the bedforms.

From this, it follows that the effective normal stress at the bed is given by
2.2
where
2.3
*ρ*_{w} is the density of water, *τ*_{nn} is the deviatoric normal stress in the ice and *Π* is the reduced pressure defined by
2.4
where *p*_{i} is ice pressure. Equation (2.2) presumes that water at the ice–till interface is connected hydrostatically to the stream system. (Note that we presume it is the far-field effective pressure (not normal stress) that is determined by the stream drainage.)

#### (i) Ice flow

The geometry of the flow is shown in figure 1. We use a coordinate system (*x*,*y*,*z*), where *z* is vertical and *x* points in the direction of the downstream ice flow. The ice surface is *z*=*z*_{i}, and the ice base is at *z*=*s*. In general, if the ice separates from the till surface, we have a cavity or stream of depth *h*_{w}. In the present paper, dealing as it does with linear stability of an undisturbed flow in which there are no cavities, we suppose that *h*_{w} is always zero.

The equations describing the slow three-dimensional flow of an incompressible fluid of viscosity *η*_{i} can be written in the form
2.5
where *z*_{i} is the ice surface, and the velocity is ** u**=(

*u*,

*v*,

*w*). The surface is almost constant on the drumlin length scale, but there is a crucial small slope that provides the driving basal stress, and this must be retained. Because the landforms that we study have a natural length scale which is much smaller than the horizontal extent of an ice sheet, the curvature or divergence of the ice flow streamlines is largely irrelevant to their study, and for this reason, we suppose that the ice flow in the absence of any basal topography is a two-dimensional flow. In this basic state, we therefore take 2.6 where

*d*

_{i}is local ice depth and

*S*is the surface slope. In the unperturbed flow, the

*x*-axis is orthogonal to the elevation contours of the upper ice surface. The slope defines a relevant ice sheet length scale 2.7 The presence of basal undulations will cause perturbations to this basic surface profile.

#### (ii) Boundary conditions

The boundary conditions deviate from those considered by Fowler (2009) and Schoof (2007), insofar as we apply conditions of no stress at the ice surface, together with a kinematic condition:
2.8
where *B* is the accumulation rate, and we define *τ*_{1} and *τ*_{2} as the two (*x* and *y*) components of tangential shear stress. We calculate these in appendix A.

At the bed *z*=*s*, one of the boundary conditions is equation (2.2), i.e.
2.9
which gives *N* in terms of *s*, *Π* and *τ*_{nn}.

The magnitude of the basal shear stress at the bed is
2.10
where again *τ*_{1} and *τ*_{2} are the two tangential components of the shear stress at the bed. The two *x*-tangential and *y*-tangential components of the basal velocity are
2.11
The sliding law at the bed is then taken to be
2.12
where *U*(*τ*,*N*) is the sliding velocity. We will also use the sliding law in the form that is more commonly determined by sliding theories; that is, we invert *U* to find *τ* as a function of *U* and *N*:
2.13
The final condition at the base is the kinematic condition,
2.14
this ignores basal melting, which is negligible in this context. (Basal melting is of the order of mm yr^{−1}, whereas accumulation is of the order of 0.1–1 m yr^{−1}, which, as we show below, is also negligible.)

#### (iii) Till flux

The evolution of the bed *s* is determined by the Exner equation
2.15
in which **q** is the till flux, given by
2.16
where *a* is the depth of deforming till and **V** is the mean till velocity. We suppose that till deformation is caused by a combination of a shearing component and a squeezing component, as described in appendix B:
2.17
where **u**|_{z=s} is the ice velocity evaluated at the bed, *c*≤1, and we suppose that deforming till depth is constrained by the depth at which the increase of effective pressure brings the stress below the yield stress. This leads to the prescription (see Fowler 2009 for details)
2.18
where *μ* is the coefficient of friction, *ϕ* is the porosity of the till and Δ*ρ*_{sw}=*ρ*_{s}−*ρ*_{w} is the difference between sediment and water density.

In reality, following a change in *N* at the ice–till interface, the effective pressure profile with depth will change over a finite time, controlled by the permeability and compressibility of the till. The consolidation equation for adjustment of the till effective pressure follows from the equations
2.19
where *k* is till permeability and *η*_{w} is the viscosity of water; thus (assuming material quantities are constant)
2.20
where *ϕ*(*N*) is determined by the till compressibility (see Clarke 1987 for details). Equation (2.6) determines a time scale for adjustment of effective pressure of the form
2.21
where *d*_{T} (defined below) is an effective measure of deformable till thickness. Since *t*_{T} will be found to be small, the equilibrium prescription (2.18), may be assumed. Fowler (2009) used a modification of equation (2.18), which allowed relaxation towards the equilibrium value, and which included also a flux term **∇**⋅**q**, on the basis of an analogy with erosion and deposition of sediments in rivers.

#### (iv) Non-dimensionalization

We define two depth scales. The first is associated with drumlin elevation, and is given by
2.22
where Δ*ρ*_{iw}=*ρ*_{w}−*ρ*_{i} is the difference between water and ice density. The second is given by
2.23
and is a measure of deforming till depth. We also define a representative length scale *l*, which is determined by balancing the stress and strain rates, thus
2.24
where *u*_{0} is an appropriate ice velocity scale.

We take the basic shear flow without bed perturbations to be
2.25
where
2.26
is the basal shear stress of the underlying ice flow; **i** is the unit vector in the *x*-direction. Note that this is not an exact solution for the boundary conditions on the surface, but is an accurate approximation for small slope *S*. When the aspect ratio *d*_{D}/*l* of the emerging bedforms is small (as is the case in practice), then the perturbation to the velocity is correspondingly small.

We now scale the model by choosing
2.27
the dimensionless velocity is introduced here because the development of the topography causes a change in the mean sliding velocity. Its prescription is given below. From equation (2.18), the dimensionless equilibrium deforming till depth *A* is given by
2.28

With this choice of scaling, the dimensionless model for the ice flow is
2.29
where
2.30
It is convenient to define the perturbed surface profile *H* by putting
2.31
so that the momentum equation becomes
2.32
and note that in the basic unperturbed state, *H*=0. Note also that *h* is constant to *O*(*S*), and this allows simplification of the surface boundary conditions (2.8), adopting the definitions in equations (A2) and (A3), and bearing in mind that the slope terms are of *O*(*S*/*σ*)≪1, to the form
2.33
where
2.34
and the second equation in (2.33), the kinematic condition, can be written in the approximate form
2.35
where
2.36

Using the definitions in equations (A2) and (A3), the basal conditions take the form
2.37
and these are all applied at *z*=*νs*; *β* is defined by
2.38
where we use the magnitude of *b* from (B6), that is,
2.39

With *l*=300 m, *d*_{i}=1500 m, thus *τ*_{b}=0.15 bar with an assumed ice surface slope of *S*=10^{−3}, *d*_{D}=50 m, *N*_{c}=0.4 bar, *d*_{T}=5 m, *u*∼100 m yr^{−1}, typical values of the parameters are
2.40
Additionally, if we take *η*_{w}=10^{−3} Pa s, |*ϕ*′(*N*)|∼10^{−7} Pa^{−1}, *k*∼10^{−15} m^{2} (Clarke 1987), then the till adjustment time scale in equation (2.21) is *t*_{T}∼0.1 yr, and the ratio of this to the drumlin formation time scale *t*_{D}=*d*_{D}*l*/*d*_{T}*u*_{0} in equation (2.27) is *t*_{T}/*t*_{D}∼0.3×10^{−2}, corroborating our assumption of equilibrium deformable till thickness in equation (2.28).

## 3. A reduced model

We now simplify the model by assuming that the aspect ratio *ν* is small. This has two effects. The first is that the terms in *ν* in equation (2.37) can be ignored. The second is that the basal boundary condition can be applied at *z*=0. Of the other parameters, we put *S*, *ε* and *δ* to zero, but we retain *λ*, *α* and *β*. *α* is retained as it is known to provide a regularizing effect on the two-dimensional instability, and indeed to control its growth time. *λ* represents the time scale for adjustment of the surface, so its incluson allows assessment of the stability of the surface. Finally, *β* represents the small lateral till flux term, which is a singular perturbation to the basic model, and is expected to be important when cavities form, as it describes cavity infilling owing to till squeezing.

The reduced model equations for the ice flow are then
3.1
and the boundary conditions become
3.2
while on the base,
3.3
all now applied at *z*=0. It should be noted that the mean sliding velocity is a function of time only, and is to be determined by ensuring that the spatial average of *u* is zero.

### (a) The Exner evolution equation

We can uncouple the determination of the ice flow from that of the base in the following way. The evolution of the bed is found by solving the Exner set of equations for *s*, *N* and *H*, all functions of *x*, *y* and *t*, given by
3.4
where *Θ* and *Ξ* are determined once the ice flow problem is solved by defining
3.5

It is appropriate to comment at this point on the equation for the effective pressure in equation (3.4), or its dimensional equivalent in equation (2.2). The choice of scales for the model assumes that the channel effective pressure *N*_{c} is constant, which is consistent with simple drainage theories, but this is unlikely to be true in general. If we suppose that *N*_{c} varies, then we may use the effective pressure at the origin,
3.6
as the scale for the stresses. If we assume is constant, for example, we would replace equation (3.4)_{1} by
3.7
and this would have a mild effect on the subsequent stability calculation. Of course, the assumption of constant channel water pressure ignores the hydraulics of water flow through the system. A proper resolution of this issue requires a proper discussion of dynamic drainage, and as we discuss in the conclusions, this will form a focus for future development of the model. For the moment, we retain the simpler assumption inherent in equation (3.4)_{1}.

### (b) Ice flow solution

We can use the Fourier transform of a function
3.8
where **k**=(*k*_{1},*k*_{2}), to solve the ice flow equations. This has been done by Schoof 2002*a*), for example. First we define
3.9
The solution for the transformed flow variables is
3.10
where *k*=|**k**| and *α*^{±} and *β*_{±} are functions of **k**.

The eight coefficients that constitute *α*^{±} and *β*_{±} are found by applying the transformed boundary conditions, which take the form
3.11
together with two conditions that arise from requiring consistency of equation (3.10) with the transformed continuity equation
3.12
The primes here denote derivatives with respect to *z*. These eight equations and their solution are given in electronic supplementary material, appendix C, and from these we can construct expressions for the transforms of *Θ* and *Ξ*, for incorporation into equation (3.4). These take the form
3.13
where the coefficients *G*_{k} and Δ_{5} depend on *k*_{1}, *k*_{2} and *σ* only, and are given explicitly in electronic supplementary material, appendix C. Together with equations (3.4) and (3.9), this provides a closed nonlinear system for the evolution of *s* and *H*.

### (c) Linear stability

Equations (3.4), (3.9) and (3.13) have a steady uniform solution in which *s*=0, *N*=1 and *H*=0, and we suppose also that and *A*=1. The nonlinearity in the model is in the advective term in equation (3.4)_{2}, and the sliding law function *F*. We write
3.14
assume *s*, *P* and *H* are small, linearize the terms *A* and *F* and Fourier transform the resulting system, denoting transforms by overhats. We then write
3.15
etc., where
3.16
and *r* is the growth rate and *c*_{w} is the wave speed. We thus obtain
3.17
assuming that .

Eliminating the redundant variables, we can write this as a homogeneous system of three linear equations for , and , and the condition for a solution is that the determinant should vanish. This leads to a quadratic equation for *Σ* in the form
3.18
where
3.19
and the new coefficients *Q* and *m*_{j} are defined in electronic supplementary material, appendix C.

#### (i) The limit *σ*→0

In the limit of infinitely deep ice, *σ*→0 and the coefficients *Q* and *m*_{j} tend to one. In this limit, the quadratic (3.18) can be factorized, and the two roots are
3.20
The first of these is a surface mode, and describes the relaxation of the surface to equilibrium on a rapid time scale *λ*. Note that at large wavenumbers, this mode apparently becomes neutral; as we shall see in a moment, this is not practically the case.

The second mode is that describing the bed instability. Computing its real part, we find that the growth rate is
3.21
and instability occurs if *A*′>0, as has been previously found (e.g. Fowler 2009). The till squeezing coefficient *β* acts diffusively, and is weakly stabilizing. For small values of *α* and *β*, the growth rate is
3.22
Note that purely longitudinal (MSGL) modes are neutrally stable (*r*=0).

#### (ii) *σ*≪1

For finite *σ*, the quadratic can almost be factorized; if we take *m*_{1}=*m*_{4} and *m*_{2}=*m*_{3}, then the roots are
3.23
and we can see from electronic supplementary material, appendix C, that, since *m*_{1}−*m*_{4} and *m*_{2}−*m*_{3} are relatively small, these are good approximations.

Since *Q*(0)=0, the modification of *Σ*_{1} shows that the surface mode is stable at small wavenumbers. Specifically, *Q*∼(*k*/2*σ*) for small *k*, and thus *Σ*_{1}∼(−1/4*λσ*). For small *α* and *β*, the basal mode growth rate is modified to
3.24
where the coefficients *q*_{j} are defined in electronic supplementary material, appendix C. Since *q*_{2}<0 and *q*_{4}>−1, the system is again stabilized at small wavenumbers, but the instability criterion is still essentially that *A*′>0.

#### (iii) *λ*≪1

A more directly accurate approximation uses the fact that the parameter *λ*≪1. The instability mode is then
3.25
(which is indeed equation (3.23)_{2} if *m*_{1}=*m*_{4} and *m*_{2}=*m*_{3}). Ignoring the small terms in *β*, the growth rate is
3.26
with the same conclusion as before.

#### (iv) *k*_{1}=0

Finally, we consider the exact quadratic for purely longitudinal modes with *k*_{1}=0. In this case, one readily finds that the coefficients of the quadratic are all real and positive, and therefore both modes are stable, and approximately,
3.27

However we look at it, we see that the finite depth hardly affects the conclusions of the infinite depth theory. Further, consultation of *Q*, *q*_{3} and *q*_{4} in figures S1–S3 of electonic supplementary material, appendix C, indicates that the infinite depth theory (in the form (3.25)) should be accurate for , i.e. for dimensional wavelengths
3.28
This seems to be well satisfied in practice, suggesting that the infinite depth theory is sufficient. Note that even with the supposition of large depth, it is still possible to calculate the surface perturbation, which is found (see equations (3.4) and (3.13) by equating *Ξ* to zero.

## 4. Discussion

The instability theory of drumlin formation (Hindmarsh 1998*a*; Fowler 2000; Schoof 2007; Dunlop *et al.* 2008) has been successful in predicting two-dimensional transverse waveforms interpretable as ribbed moraine. However, drumlins are three dimensional, and MSGLs, while essentially two dimensional, are aligned with the ice flow direction. There has been a focus on MSGLs recently, particularly as they have been recently found under the Rutford ice stream (Smith *et al.* 2007; King *et al.* 2009; Smith & Murray 2009), and a number of alternative ideas for their formation have been presented (Clark *et al.* 2003; Schoof & Clarke 2008). Of these, the paper by Schoof and Clarke proposes an instability based on non-Newtonian ice rheology, a radically different idea to that presented here, but one which is able to produce longitudinal rolls.

While it must be emphasized that the absence of a three-dimensional linear instability does not at all preclude the possibility that finite amplitude computations will yield three-dimensional wave forms, it would certainly be easier for the instability theory if such three-dimensional instabilities did exist. Particularly, MSGLs are very suggestive of purely two-dimensional longitudinal rolls, as opposed to the transverse ridges of ribbed moraine.

A particular motivation of the present paper was therefore to examine whether the instability theory would allow three-dimensional modes to grow when the ice is assumed to be of finite depth. If the instability theory is to provide a unified theory for the formation of subglacial bedforms, then it must be able to explain not only the formation of ribbed moraine, but also drumlins and MSGLs. A secondary motivation has been to examine the conditions of accuracy of the infinite depth theory, and indeed to provide a framework for future studies, as well as to integrate into the theory the signal which the bedforms produce at the ice surface.

The instability theory in the form presented by Schoof (2007) and Fowler (2009) has been adapted here to allow for the squeezing of till under effective pressure gradients. This is found to be a small effect (the parameter *β*≪1), and one which is weakly stabilizing. If the present theory were to be applied to the computation of finite amplitude bedforms in the manner outlined by Fowler (2009), then this squeezing term will become significant, since the presence of cavities allows the effective pressure to jump discontinuously. We would expect that the squeezing term would replace these sharp jumps with steep boundary layers, and in particular, till will be squeezed into cavities.

Inclusion of a finite depth allows the determination of the perturbation of the ice surface in terms of the bedrock undulations below. This is a subject that has been of quite some interest, insofar as it provides an indirect method for potentially assessing basal conditions, and in particular the variability of the sliding velocity at the base (Hindmarsh 1998*b*; Gumundsson 2003, 2008; Raymond & Gumundsson 2005; Schoof 2002*b*), but most previous work has not incorporated the dependence of the sliding law on the effective hydraulic pressure (an exception is the paper by Hindmarsh (1998*b*)). In the present theory, we find that the surface relaxes rapidly to an equilibrium that is determined by *Ξ*=0 in equation (3.13)_{2}, that is,
4.1
and using the definitions of the coefficient functions in electronic supplementary material, appendix C, this can be written in the form
4.2
where *κ*=*k*/*σ*, and
4.3

The real transfer functions *G*_{SC} and *G*_{SB} represent the effects on the ice surface of perturbations to the basal sliding velocity and to the bed surface, and are plotted in figures 2 and 3. Note in particular that purely longitudinal rolls (MSGLs) have no effect on the surface in this theory. Figures 2 and 3 show that, as we might expect, short wavelength features have relatively little effect on the surface. From equation (2.31), the amplitude of these surface perturbations is of *O*(*d*_{i}*S*), which, for ice of depth 1000 m and slope 10^{−3}, is 1 m.

## 5. Conclusions

In this paper, we have extended the instability theory of drumlin formation to the case of finite ice depth, and we have examined the stability of fully three-dimensional perturbations.

In essence, we repeat calculations of Hindmarsh (1998*b*), with the distinction that our results are analytic. We find that the infinite depth theory is accurate for bedform wavelengths less than the ice thickness, but that even for longer wavelengths, the essence of the theory is not substantially altered, either qualitatively or quantitatively. In particular, the instability remains fundamentally two dimensional, and associated with the formation of transverse ridges, which we identify with ribbed moraine. Purely longitudinal rolls (corresponding to MSGLs) are essentially neutrally stable, although the inclusion of the effects of till squeezing renders them weakly stable. The instability theory in its present form therefore does not allow for the direct formation of MSGLs, although the possibility that they form through a secondary instability of transverse ribbed moraine still exists. Thus, we support Schoof’s (2007) scepticism as to whether the instability theory is a realistic candidate to explain the variety of subglacial bedforms.

An issue for the future is thus whether a more sophisticated version of the instability theory will have the capacity to produce drumlins and MSGLs directly. The key requirement would seem to be the necessity of an unstable purely longitudinal mode. The key physical simplification that has been made in the theory is that the stream flow is assumed to be passive: we suppose that there is an articulated stream system that passes through the evolving bedforms without being affected by them. It seems fairly evident that this assumption is unrealistic, since the development of transverse ridges will act as a barrier to water flow. Therefore, the next step in the development of the theory is to pose a more realistic model for stream flow, which allows it to interact dynamically with the deforming sediment.

## Supplementary material

The five files of supplementary material are an Appendix C, which gives the results of solving (3.10), and the computations of the coefficients *G*_{k}, Δ_{5}, *Q* and *m*_{j} as used in (3.13) and following. The file supplementary pdf gives the iterative equations which are solved in computing these coefficients, and msgl.gnu gives a gnuplot_le (described in Description of msgl.gnu) which was used for numerical checking and computation. There is also a nomenclature.

## Acknowledgements

I acknowledge the financial support of the Mathematics Applications Consortium for Science and Industry (`www.macsi.ul.ie`) through the Science Foundation Ireland mathematics initiative grant 06/MI/005. The support of NERC through grant NE/D013070/1, {*Testing the instability theory of subglacial bedform production*}, is also acknowledged. For continuing fruitful discussions, my thanks to Chris Clark, Chris Stokes, Felix Ng, Heike Gramberg, Matteo Spagnolo, Paul Dunlop and Richard Hindmarsh.

## Appendix A. Interfacial stress components

To calculate the normal and tangential components of stresses and velocities at an interface, for example at the bed *z*=*s*(*x*,*y*,*t*), we use the normal, *x*-tangential and *y*-tangential vectors, which we define as
A1
where **∇** *s*=(*s*_{x},*s*_{y}). From this, we find that the normal deviatoric stress is defined by
A2
The two tangential shear stresses at the ice–till interface (*τ*_{1},*τ*_{2}), where *τ*_{i}=**n**⋅** τ**⋅

**t**

_{i}, are A3 (Analogous expressions hold at

*z*=

*z*

_{i}, with

*s*replaced by

*z*

_{i}.)

## Appendix B. Sediment flux

That subglacial sediment can be transported by ice sheets, we consider to be undeniable. It follows that till may be deformable, and the issue then arises of describing an appropriate and effective rheology for it. Boulton & Hindmarsh’s (1987) early expression of a nonlinear viscous rheology was later refuted by Kamb (1991), who proposed a plastic rheology. The debate led to controversy (Tulaczyk *et al.* 2000; Iverson & Iverson 2001; Fowler 2003), partly fuelled by a lack of appreciation of what plastic flow actually consists of. Few would argue that till, as a granular material, has a yield stress, but how the till deforms when this yield stress is reached still requires a rheological choice. Although this is still a subject of current research, viscoplastic rheologies such as that of Jop *et al.* (2006) suggest that plastic till flux can be described by an effective viscous flow law, in which the effective viscosity is determined by the yield stress.

Let us suppose that this is the case, and that till flows as a Couette–Poiseuille flow, driven by a shear stress *τ* and an effective pressure gradient (for example, in the *x* direction) −*N*′. If the deformable till is of depth *a* and effective viscosity *η*_{T}, then one finds that the till velocity at the ice–till interface is
B1
while the sediment flux is
B2

We anticipate that *N*′*a*≪*τ*; then it follows that
B3
and the sediment flux can be written as
B4
If we suppose that there is no slip between till and ice, so that *u*_{b}=*u* is the basal ice velocity, then we may generalize equation (B4) to three dimensions as
B5
where we write **q**=*a***V**, and here
B6
note that equation (B5) is independent of the effective till viscosity. The term in **∇***N* in equation (B5) represents also a small correction, but one that turns out to be singular, and therefore needs to be retained for that reason.

## Footnotes

↵1 By longitudinal rolls, we mean sinusoidal wave forms whose axes are aligned with the ice flow direction.

- Received January 12, 2010.
- Accepted February 10, 2010.

- © 2010 The Royal Society